Issue 
A&A
Volume 676, August 2023



Article Number  A42  
Number of page(s)  18  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202346155  
Published online  04 August 2023 
Tensortoscalar ratio forecasts for extended LiteBIRD frequency configurations
^{1}
Institute of Theoretical Astrophysics, University of Oslo, Blindern, Oslo, Norway
^{2}
IRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France
^{3}
International School for Advanced Studies (SISSA), Via Bonomea 265, 34136 Trieste, Italy
^{4}
INFN Sezione di Trieste, Via Valerio 2, 34127 Trieste, Italy
^{5}
IFPU, Via Beirut, 2, 34151 Grignano, Trieste, Italy
^{6}
Université de Paris, CNRS, Astroparticule et Cosmologie, 75013 Paris, France
^{7}
Instituto de Astrofísica de Canarias, 38200 La Laguna, Tenerife, Canary Islands, Spain
^{8}
Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{9}
Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 2778583, Japan
^{10}
NIST Quantum Sensors Group, 325 Broadway, Boulder, CO 80305, USA
^{11}
National Astronomical Observatory of Japan, Mitaka, Tokyo 1818588, Japan
^{12}
Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 2, Roma, Italy
^{13}
INFN Sezione di Roma, P.le A. Moro 2, 00185 Roma, Italy
^{14}
University of Milano Bicocca, Physics Department, Piazza della Scienza, 3, 20126 Milano, Italy
^{15}
INFN Sezione Milano Bicocca, Piazza della Scienza, 3, 20126 Milano, Italy
^{16}
Instituto de Fisica de Cantabria (IFCA, CSICUC), Avenida los Castros s/n, 39005 Santander, Spain
^{17}
Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, School of Natural Sciences, The University of Manchester, Oxford Road, Manchester M13 9PL, UK
^{18}
SLAC National Accelerator Laboratory, Kavli Institute for Particle Astrophysics and Cosmology (KIPAC), Menlo Park, CA 94025, USA
^{19}
Stanford University, Department of Physics, Stanford, CA 943054060, USA
^{20}
Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK
^{21}
Dipartimento di Fisica e Scienze della Terra, Universitá di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
^{22}
INFN Sezione di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
^{23}
INAF – OAS Bologna, Via Piero Gobetti, 93/3, 40129 Bologna, Italy
^{24}
Dipartimento di Fisica e Astronomia “G. Galilei”, Universitá degli Studi di Padova, Via Marzolo 8, 35131 Padova, Italy
^{25}
INFN Sezione di Padova, Via Marzolo 8, 35131 Padova, Italy
^{26}
INAF, Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy
^{27}
Dipartimento di Fisica, Universitá degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy
^{28}
INFN Sezione di Milano, Via Celoria 16, 20133 Milano, Italy
^{29}
School of Physics and Astronomy, Cardiff University, Cardiff CF24 3AA, UK
^{30}
Dipartimento di Fisica, Universitá di Roma Tor Vergata, Via della Ricerca Scientifica, 1, 00133 Roma, Italy
^{31}
INFN Sezione di Roma2, Universitá di Roma Tor Vergata, via della Ricerca Scientifica, 1, 00133 Roma, Italy
^{32}
University of California, Berkeley, Department of Physics, Berkeley CA 94720, USA
^{33}
University of California, Berkeley, Space Sciences Laboratory, Berkeley, CA 94720, USA
^{34}
Lawrence Berkeley National Laboratory (LBNL), Computational Cosmology Center, Berkeley, CA 94720, USA
^{35}
Kavli Institute for Particle Astrophysics and Cosmology (KIPAC), Stanford University, Stanford, CA 94305, USA
^{36}
Centre Spatial de Liège (STAR Institute), University of Liège, Avenue du PréAily, Angleur 4031, Belgium
^{37}
Institute of Particle and Nuclear Studies (IPNS), High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 3050801, Japan
^{38}
International Center for Quantumfield Measurement Systems for Studies of the Universe and Particles (QUP), High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 3050801, Japan
^{39}
Dpto. de Física Moderna, Universidad de Cantabria, Avda. los Castros s/n, 39005 Santander, Spain
^{40}
INFN Sezione di Bologna, Viale C. Berti Pichat 6/2, 40127 Bologna, Italy
^{41}
Tohoku University, Graduate School of Science, Astronomical Institute, Sendai 9808578, Japan
^{42}
Japan Aerospace Exploration Agency (JAXA), Institute of Space and Astronautical Science (ISAS), Sagamihara, Kanagawa 2525210, Japan
^{43}
The Graduate University for Advanced Studies (SOKENDAI), Miura District, Kanagawa, 2400115 Hayama, Japan
^{44}
Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver BC V6T1Z1, Canada
^{45}
Institut d’Astrophysique de Paris, CNRS/Sorbonne Université, Paris, France
^{46}
Lawrence Berkeley National Laboratory (LBNL), Physics Division, Berkeley, CA 94720, USA
^{47}
Laboratoire de Physique de l’École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France
^{48}
Space Science Data Center, Italian Space Agency, Via del Politecnico, 00133 Roma, Italy
^{49}
Université ParisSaclay, CNRS, Institut d’Astrophysique Spatiale, 91405 Orsay, France
^{50}
Gran Sasso Science Institute (GSSI), Viale F. Crispi 7, 67100 L’Aquila, Italy
^{51}
David A. Dunlap Department of Astronomy and Astrophysics, 50 St. George Street, Toronto, ON M5S3H4, Canada
^{52}
Institute of Astrophysics, Foundation for Research and TechnologyHellas, Vasilika Vouton, 70013 Heraklion, Greece
^{53}
Department of Physics and ITCP, University of Crete, 70013 Heraklion, Greece
^{54}
INAFOsservatorio Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius, Italy
^{55}
Physics and Astronomy Dept., University College London (UCL), London, UK
^{56}
The University of Tokyo, Department of Astronomy, Tokyo 1130033, Japan
^{57}
Suwa University of Science, Chino, Nagano 3910292, Japan
^{58}
INFN Sezione di Pisa, Largo Bruno Pontecorvo 3, 56127 Pisa, Italy
^{59}
Okayama University, Department of Physics, Okayama 7008530, Japan
^{60}
Université ParisSaclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France
^{61}
National Institute of Technology, Kagawa College, 355 Chokushicho, Takamatsu, Kagawa 7618058, Japan
^{62}
UniversitätsSternwarte, Fakultät für Physik, LudwigMaximilians Universität München, Scheinerstr.1, 81679 München, Germany
^{63}
NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
Received:
15
February
2023
Accepted:
10
May
2023
LiteBIRD is a planned JAXAled cosmic microwave background (CMB) Bmode satellite experiment aiming for launch in the late 2020s, with a primary goal of detecting the imprint of primordial inflationary gravitational waves. Its current baseline focalplane configuration includes 15 frequency bands between 40 and 402 GHz, fulfilling the mission requirements to detect the amplitude of gravitational waves with the total uncertainty on the tensortoscalar ratio, δr, down to δr < 0.001. A key aspect of this performance is accurate astrophysical component separation, and the ability to remove polarized thermal dust emission is particularly important. In this paper we note that the CMB frequency spectrum falls off nearly exponentially above 300 GHz relative to the thermal dust spectral energy distribution, and a relatively minor high frequency extension can therefore result in even lower uncertainties and better model reconstructions. Specifically, we compared the baseline design with five extended configurations, while varying the underlying dust modeling, in each of which the HighFrequency Telescope (HFT) frequency range was shifted logarithmically toward higher frequencies, with an upper cutoff ranging between 400 and 600 GHz. In each case, we measured the tensortoscalar ratio r uncertainty and bias using both parametric and minimumvariance componentseparation algorithms. When the thermal dust sky model includes a spatially varying spectral index and temperature, we find that the statistical uncertainty on r after foreground cleaning may be reduced by as much as 30–50% by extending the upper limit of the frequency range from 400 to 600 GHz, with most of the improvement already gained at 500 GHz. We also note that a broader frequency range leads to higher residuals when fitting an incorrect dust model, but also it is easier to discriminate between models through higher χ^{2} sensitivity. Even in the case in which the fitting procedure does not correspond to the underlying dust model in the sky, and when the highest frequency data cannot be modeled with sufficient fidelity and must be excluded from the analysis, the uncertainty on r increases by only about 5% for a 500 GHz configuration compared to the baseline.
Key words: ISM: general / cosmology: observations / cosmic background radiation / polarization / cosmological parameters / Galaxy: general
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
One of the key predictions of the current cosmological inflationary paradigm is the existence of a stochastic background of primordial gravitational waves created shortly after the Big Bang (Starobinsky 1980; Sato 1981; Guth 1981; Albrecht & Steinhardt 1982; Linde 1982, 1983). If such waves do exist, they should induce a particular and unique polarization signature in the cosmic microwave background (CMB) on large angular scales, corresponding to socalled Bmode or divergencefree polarization (Kamionkowski et al. 1997). Detecting this imprint in the CMB ranks among the top observational priorities in modern cosmology, and huge efforts are currently being made to develop the necessary instrumentation, data analysis techniques, and theoretical modeling required for this task.
Among the leading international Bmode efforts is LiteBIRD (the Lite (Light) satellite for the study of Bmode polarization and Inflation from cosmic background Radiation Detection), a satellite concept chosen by JAXA as one of their top priorities for the coming decade, aiming for a launch from 2028 to 2029 (Hazumi et al. 2020; LiteBIRD Collaboration 2022). The development of LiteBIRD has already been ongoing for more than 13 years, evolving gradually from a relatively simple and lightweight concept that originally included only six frequency channels between 60 and 280 GHz (Hazumi et al. 2012; Matsumura et al. 2014) into the current baseline configuration spanning 15 frequency bands between 40 and 402 GHz (Hazumi et al. 2020). This increased sophistication came as a direct response to observational insights gained regarding the astrophysical sky as measured by CMB experiments such as Planck (Planck Collaboration I 2020) and BICEP2/Keck (Ade et al. 2016). These observations pointed toward a complex reality in which both polarized synchrotron and thermal dust emission must be modeled and subtracted with exquisite precision, while simultaneously accounting for unpolarized complications such as carbon monoxide and timedependent zodiacal light emission.
In the following, we revisit the upper limit of the LiteBIRD frequency range, and examine its implications for the instrument’s ability to distinguish between CMB and astrophysical foregrounds. While LiteBIRD Collaboration (2022) has already demonstrated that the baseline configuration fully achieves the primary mission goal – namely to constrain the tensor to scalar ratio with the total uncertainty, δr, down to δr < 0.001 – it is still desirable to optimize the uncertainties from component separation, which corresponds to onethird of the total budget error on r, in order to increase the overall mission margins. Exploring that is the main motivation for the current work. At the same time, we emphasize that this paper only considers the raw sensitivity and component separation aspects of the instrumental design. We need to perform a comprehensive trade study before the collaboration can decide on a baseline change. The scope of the current paper is simply to determine whether sufficiently significant gains are available to warrant such a study.
To this end, we compared the current baseline configuration with five alternative configurations in which the upper limit was varied between 400 and 600 GHz in steps of 40 GHz. The absolute upper center frequency limit of 600 GHz was dictated by the current bolometer design, which has a hard cutoff at 680 GHz (Novotny & Meincke 1975). To additionally minimize the number of changes required for the overall satellite design, we chose to only modify one of the three LiteBIRD telescopes, namely the HighFrequency Telescope (HFT); the Low and MidFrequency Telescopes (LFT and MFT) have been left unchanged. Frequencies up to 600 GHz can be considered only because LiteBIRD is a space mission, while groundbased experiments are restricted to frequencies below 300 GHz due to the atmospheric windows (Liebe 1981; Pardo et al. 2001).
The main question considered in this paper is the quantitative relationship between the frequency range and the tensortoscalar ratio uncertainty marginalized over polarized thermal dust emission. Previous analyses (e.g., Remazeilles et al. 2016) have suggested that a low cutoff on the high frequency side only supports very limited ability to constrain the spectral energy distribution (SED) of thermal dust emission. At the same time, it is well known that the CMB SED falls nearly exponentially above 300 GHz relative to the thermal dust SED (see e.g., Fig. 35 of Planck Collaboration IV 2020), and therefore even relatively small changes in the cutoff frequency may have a dramatic impact on the confusion between CMB and thermal dust. Quantifying the importance of including frequencies higher than where the CMB SED drops off in a realistic setting is the main goal of the current paper. We note that while we subsequently address LiteBIRD in particular, the main arguments are instrumentagnostic, and the primary conclusions are therefore generally applicable to any future space mission or balloonborne experiment. Indeed, similar analyses have recently been published for both the PICO (Aurlien et al. 2023) and ECHO (Sen et al. 2022) CMB satellite concepts with consistent results.
This challenge of optimizing the frequency range may be separated into two important and complementary aspects, namely optimizing our sensitivity to the tensortoscalar ratio overall and maximizing our ability to reject wrong sky models. The following analysis is organized accordingly. To address the first question, we analyze ideal foreground simulations in which the fitting model matches the simulated sky, and estimate the tensortoscalar ratio uncertainty as a function of frequency coverage. In this case, we adopt a standard onecomponent modified blackbody (MBB) thermal dust model with spatially varying spectral parameters, following closely in the footsteps of LiteBIRD Collaboration (2022). In order to ensure robustness in terms of analysisdependent details, we employed five independent componentseparation methods for this task, three of which implement different variations of parametric fitting (Commander – Eriksen et al. 2008, FGBuster – Errard & Stompor 2019, and Moment Expansion – Vacher et al. 2022), while the last two implement (semi)blind internal linear combination (ILC) fitting (Needlet ILC, NILC and constrained moment ILC, cMILC – Remazeilles et al. 2011, 2021).
To assess the capability of rejecting incorrect data models, we analyzed simulations based on an extended thermal dust model, namely the physically motivated silicateironcarbon model (SiFeC; Hensley & Draine 2017). We evaluated the residual χ^{2}, and compared this with the derived tensortoscalar estimates. For this study, we note that only Commander currently provides χ^{2}based goodnessoffit estimates, and we therefore only reported results for that code in this part of the paper. However, we note that work is currently ongoing on implementing similar statistics for other methods, and those results will be reported elsewhere.
It is also important to note that the foreground sky in this paper was approximated as a “single layer model,” and no attempt was made to take into account the full 3D structure of the Milky Way. Full lineofsight integration of 3D effects would introduce additional frequencydecorrelation effects (Tassis & Pavlidou 2015; Chluba et al. 2017; Pelgrims et al. 2021; Ritacco et al. 2023; Vacher et al. 2023a), and the following results therefore represent an optimistic view of the true sky in terms of complexity. In reality, there will be additional structures, both spatially and spectrally, and fully resolving these would in principle require even more data than considered here, for instance in terms of additional frequency bands, wider frequency range, or dedicated 3D constraints on the Galactic magnetic field. The following optimization results should therefore be considered as a lower bound on what is actually needed for a future production analysis.
The paper is organized as follows. In Sect. 2 we provide a brief overview of the LiteBIRD concept and current baseline, and we define a set of alternative configurations to be considered. In Sect. 3 we describe the simulations used for the present analysis, and in Sect. 4 we briefly survey the componentseparation methods used in this paper; algorithmic details are deferred to the Appendices. In Sect. 5, we present the results derived from an ideal model analysis, while in Sect. 6 we consider modeling errors. Final conclusions are drawn in Sect. 7.
2. Baseline and extended LiteBIRD instrument configurations
In this section, we discuss the LiteBIRD mission characteristics that are relevant to the present analysis. A comprehensive review of the instrument concept can be found in LiteBIRD Collaboration (2022). First, the LiteBIRD detector arrays are populated with a total of 4508 lenslet or horncoupled transition edge sensor (TES) bolometers fabricated on silicon wafers, covering 15 frequency bands between 40 and 402 GHz. The LiteBIRD bolometers provide exquisite sensitivity across this full frequency range, but the technology itself is limited to frequencies < 678 GHz, corresponding to the superconducting gap of niobium – the material used for mmwave transmission lines on the detectors (Novotny & Meincke 1975). This value therefore naturally defines a clear upper limit for the feasible frequency range of the LiteBIRD instrument in any modification to the baseline.
As a second consideration, to ensure robust separation between the actual polarization signal and a wide range of intensityspecific contaminants, LiteBIRD will employ halfwave plates (HWPs) spinning at a rate of 0.5 to 0.8 Hz. Important strengths of mesh filterbased HWPs are low mass requirements and high transmission. However, an important drawback is finite effective bandwidth, typically corresponding to a ratio between its maximum and minimum effective frequency of about 2.3 (LiteBIRD Collaboration 2022). As such, the full relative LiteBIRD frequency range of 402 GHz/40 GHz ≈ 10 cannot be supported within a single telescope, but rather three separate telescopes are required. For LiteBIRD these are called the Low, Mid, and HighFrequency Telescopes, respectively, or LFT, MFT, and HFT as a shorthand. While this organization comes at a significant cost in terms of system complexity, it also allows significant optimization of each system, and different technologies may be used as appropriate.
In total, the three telescopes span 22 independent frequency channels, covering 15 frequency bands; there is an overlap of three channels inside the LFT, three channels between the LFT and MFT, and an overlap of one channel between the MFT and HFT. These overlaps represent important crosschecks against systematic errors in any of the three telescopes. Each frequency band overlaps the adjacent bands by approximately half the bandwidth. The radiation is detected by seven different wafer types, each containing one or two types of multichromatic pixels. Each pixel, for each polarization state, distributes power through one, two, or three different bandpass filters to independent TES detectors. The HFT detectors are implemented in terms of three different detector modules; one with 3color pixels, one with 2color pixels, and the highest frequency band on singlecolor pixels (Montier et al. 2020).
The main goal of the current paper is to investigate the optimal LiteBIRD frequency range with respect to the removal of the thermal dust from the CMB polarization signal. Given the mature state of the current baseline, modifications must be realistic and the scientific improvements for doing this must be substantial. One of the parameters that may still be adjusted, however, is the effective HFT frequency range. Specifically, it is possible to shift the entire frequency range by a constant factor, which ensures that the relative factor between the highest and lowest frequency bands is unchanged, as required by the meshfilter HWP. In this paper, we therefore defined a set of extended HFT configurations parametrized by a scaling factor, f, that varied between 1.0 and 1.5 in steps of 0.1, and shifted the entire HFT frequency range. These models were denoted Mi, such that M0 corresponded to the baseline, while M5 corresponded to an extended model with a highest frequency of 1.5 × 402 GHz ≈ 603 GHz. For each modified frequency channel, we calculated new instrumental parameters, such as sensitivity or beam size, using the same methodology as used for the baseline in LiteBIRD Collaboration (2022). This included evaluating the system temperature through ray tracing from the feed to terminals, either inside the telescope or the sky. The beam full width at half maximum (FWHM) values were estimated via the edge taper at the telescope aperture, as calculated from the feed farfield pattern. The resulting values are summarized numerically in Table 1, while Fig. 1 shows the sensitivity as a function of frequency for the extended bands.
Fig. 1. Sensitivity versus frequency for the various LiteBIRD configurations listed in Table 1. The black crosses indicate LFT or MFT frequencies that are not modified in this paper, while colored symbols indicate the various modified HFT frequencies. 
Summary of LiteBIRD HFT instrument configurations considered in this paper in terms of frequency, beam size, and sensitivity.
It is interesting to note that the sensitivity of two channels in different configurations can vary significantly, despite the fact that both their center frequencies and detector counts are nearly identical; a typical example of this is channel 1 in M2 and channel 2 in M0. The reason for this is that the optics for the HFT are optimized for the telescope’s full frequency bandwidth, not any single channel. Frequencies near the edge of the telescope bandwidth typically have poorer optical absorption and reflection properties than frequencies nearer the center of the band.
We also note that the LFT and MFT are left unchanged in this exercise, and their instrumental parameters are not listed in Table 1. The primary motivation for this choice is programmatic; changing the baseline at this late stage will carry a nonnegligible cost, and it is therefore highly desirable to make as few changes as possible. At the same time, it is important to note that this choice does come at a cost in the form of weaker systematics crosscheck between the MFT and HFT, since the two telescopes no longer have overlapping frequency ranges. Thus, these modified configurations may be considered as trading off robustness toward “unknown unknowns” (for instance, subtle HWP frequency effects) against robustness to a known unknown (polarized thermal dust emission). If this issue is considered sufficiently important at a programmatic level, it is of course also possible to modify the lowerfrequency telescopes, albeit at a somewhat higher cost.
3. Simulations
The LiteBIRD simulations considered in this paper are generated directly in pixel space, and account only for sky signal and white noise. The sky emission components for all 22 frequency channels are generated using the PySM (Python Sky Model; Thorne et al. 2017) package. The LFT and MFT channels are defined according to the baseline instrument configuration in LiteBIRD Collaboration (2022), while the HFT channels are generated in six versions, called M0–M5, as tabulated in Table 1. Here, M0 is the LiteBIRD baseline, and M1–M5 are the frequency extended data sets. The instrumental summary parameters listed in this table are derived from first principles using the same methodology as the baseline configuration in LiteBIRD Collaboration (2022), properly accounting for realistic optics, bandpass, and bolometer effects.
We simulate the polarized microwave emission considering three different components, namely a cosmological CMB signal and both diffuse synchrotron and thermal dust radiation from our own Galaxy. The CMB component is simulated as a Gaussian random realization of the Planck 2018 ΛCDM bestfit power spectrum (Planck Collaboration VI 2020). In the following, we set the tensortoscalar ratio, r, to zero. Obviously, a main goal of LiteBIRD is to actually measure a nonzero value of r, and it could therefore be of interest to also consider nonzero values of this parameter. However, this paper is primarily concerned with the relative performance of different instrument configurations. This comparison is not sensitive to the exact value of r, and then for the purpose of our analyses we decided to use r = 0 as reference, which also provides the benefit of not having to deal with the cosmic variance noise term. Future work focusing on absolute detection level performance should obviously revisit this with a range of alternative cosmological parameter values, including r, the optical depth of reionization τ, and others.
Galactic synchrotron radiation is physically generated by cosmicray electrons that are accelerated by the Galactic magnetic field. This signal is highly polarized with a polarization fraction that can reach 20% at intermediate and high Galactic latitudes (Planck Collaboration XXV 2016; Kogut et al. 2007). Synchrotron emission represents the dominant foreground at frequencies below about 70 GHz (Krachmalnicoff et al. 2016; Planck Collaboration IX 2016). As a first approximation, the synchrotron SED follows a power law with a spatially varying spectral index (Lawson et al. 1987; Platania et al. 1998; Bennett et al. 2003). Recent analyses support a spectral index of β_{s} ≈ −3 (Fuskeland et al. 2014, 2021), with variations on the order of 10% across the sky on degree angular scales (Krachmalnicoff et al. 2018). In the following, we adopt the s1 PySM model for synchrotron emission, corresponding to a strict powerlaw spectral behavior with spatially varying spectral index. We note that synchrotron emission is largely irrelevant for the highest LiteBIRD frequencies, and any conclusions regarding the optimal HFT frequency range will be largely independent of the synchrotron model.
The same does, however, not hold true for thermal dust emission, which is the dominant foreground component at frequencies above 70 GHz (Krachmalnicoff et al. 2016; Planck Collaboration IX 2016). Thermal dust emission is generated by vibrating nonspherical dust grains in the interstellar medium (ISM), and is polarized because the small axes of the grains are preferentially aligned parallel to the local magnetic field (Hoang & Lazarian 2016; Hensley & Draine 2023). The detailed physical processes involved in this are, however, more complicated than for synchrotron emission, and the predicted emission from any region of the ISM will depend on a wide range of local properties, including dust grain composition, the local radiation field and magnetic field structure. For robustness, we therefore consider two different thermal dust SED models in the following, namely; (1) a phenomenologically motivated single modified blackbody parametrized by a spatially varying spectral index and temperature (1 MBB; PySM model d1); and (2) a physically motivated silicateironcarbon model (SiFeC; PySM model d8; Hensley & Draine 2017). The 1 MBB model is very often used in the literature. The second model is more exotic, and is physically motivated rather than datadriven. They also differ in complexity in terms of spatial variation. The d1 model has both a spatially varying spectral index and temperature, while d8 has spatially constant spectral parameters. Future work should also consider models that include physically motivated frequency decorrelation, for instance by implementing proper multilayer 3D models. This is however beyond the scope of the current paper.
Figure 2 shows the dust SEDs for the two different models computed on 70% of the sky. Sky masks are obtained by retaining the cleanest fraction of the sky in polarized intensity at 100 GHz. As seen in the bottom panel, the relative difference between the two models can reach tens of percent across the LiteBIRD frequency range, and properly accounting for these variations will be essential for making a robust CMB extraction. For comparison, the gray and black curves show the synchrotron and CMB SEDs. The near exponential dropoff of the CMB spectrum is clearly seen above about 300 GHz, leading to a rapidly increasing ratio between thermal dust and CMB emission between 400 and 600 GHz.
Fig. 2. SEDs of the two dust models considered in this work, the 1 MBB, (red line) and the physical dust model SiFeC, (green line). The black line shows the SED of the CMB, while the gray refers to synchrotron emission. All the SEDs are in brightness temperature units and are computed considering the rms of Stokes Q and U maps on 70% of the sky. Sky masks are obtained by retaining the cleanest fraction of the sky in polarized intensity at 100 GHz. Vertical gray lines show the central frequency of the baseline LiteBIRD frequency channels, and in blue dashed lines we report the central frequency of the HFT M5 extension. The bottom panel shows the percentage difference of SiFeC with respect to MBB. 
Instrumental effects are modeled under highly idealized assumptions, since these are considered to be subdominant regarding the central question of frequency range versus thermal dust reconstruction, and rather only add computational complexity and cost. In particular, we model all instrumental bandpasses in terms of Dirac δ functions, all instrumental beams as azimuthally symmetric Gaussians, and the noise is assumed to be uncorrelated, Gaussian, and spatially isotropic. Center frequencies, beam FWHMs, and array noise equivalent temperatures (NETs) are all listed in Table 1. Under these assumptions, the simulation procedure is very straightforward, and does not require any lowlevel timedomain processing.
The data model for a given frequency band ν may be described as
where d represents a given (simulated) sky map, s denotes the sky signal, and n is instrumental noise. The simulated sky maps analyzed in this paper are based on two different thermal dust emission models. For the baseline model, we only fit the onecomponent MBB model, while the SiFeC model is used to assess modeling errors and is therefore only present in the simulations, and not in the data model in the methods. Including also CMB and synchrotron emission, Eq. (1) may then be written in the following explicit form, adopting thermodynamic CMB temperature units,
where {a_{CMB}, a_{s}, a_{d}} are the signal amplitudes relative to the reference frequency ν_{0} for each signal, γ(ν) is the conversion factor between the RayleighJeans brightness temperature and the CMB thermodynamic temperature, {β_{s}, β_{d}} are the spectral indices, T_{d} is the temperature for thermal dust, and h and k are the Planck and Boltzmann constants, respectively. In principle, this sky model is valid for both intensity and polarization, with individual parameters fitted for each Stokes parameter. In practice, however, we only include CMB, synchrotron and thermal dust emission in polarization, and we additionally assume β_{Q} = β_{U}.
All simulations are generated both at full angular resolution (adopting a HEALPix pixelization with a resolution parameter of N_{side} = 512; Górski et al. 2005) and at low angular resolution (10° FWHM, N_{side} = 16). The former set is used by the Moment expansion, FGBuster and NILC/cMILC analyses, while the latter is used by the Commander analysis. We note that the lowresolution maps do not include subpixel (or subbeam) structure, but are generated natively at the target resolution. The motivation for this choice is that subpixel structure may excite spurious Bmode power, both from beam averaging and from parallel transport inaccuracies during HEALPix downgrading. At the same time, a future LiteBIRD Commander analysis will be performed in the timedomain with full angular resolution data, following a methodology similar to that described by BeyondPlanck Collaboration (2023), and this complication will therefore not be relevant for the final LiteBIRD analysis.
4. Componentseparation algorithms
In order to ensure that the general results are robust with respect to algorithmic details, we employ a total of five different componentseparation algorithms in the following. These may be divided into two main groups, namely parametric (Commander, FGBuster, and Moment Expansion) and (semi)blind methods (NILC and cMILC).
Commander (Eriksen et al. 2008) is a Bayesian Gibbs sampler that has been successfully applied to Planck, WMAP, LiteBIRD, and many other data sets (e.g., Planck Collaboration IX 2016). The defining feature of this method is an explicit joint parametric model that simultaneously accounts for cosmological and astrophysical parameters, and this global model is explored by textbook Markov chain Monte Carlo methods. A major strength of this method is the availability of welldefined goodnessoffit χ^{2} statistics and realizationspecific uncertainty estimates, while a major weakness is a high computational cost. In this paper we only apply this method to lowresolution simulations; for further details, see Appendix A.1.
FGBuster (Puglisi et al. 2022) is based on a similar statistical foundation as Commander, but uses nonlinear optimization to explore the likelihood function rather than Monte Carlo sampling. It further speeds up the estimation process by first estimating all nonlinear spectral foreground parameters marginally with respect to linear amplitude parameters; then it estimates amplitudes conditionally with respect to the spectral parameters. Finally, angular power spectra are computed from the amplitude maps. This method has been used extensively for LiteBIRD forecasting and optimization (e.g., LiteBIRD Collaboration 2022). For further details, see Appendix A.2.
The momentexpansion method (Chluba et al. 2017; Vacher et al. 2023b) also has a similar statistical foundation as Commander and FGBuster but it allows for Taylorexpansionbased distortions in the SED models for each foreground to account for nonlinear averaging effects, for instance from lineofsight integration, beam convolution or harmonic expansion. This method has already been applied to baseline LiteBIRD simulations by Vacher et al. (2022). For further details, see Appendix A.3.
While each of the above methods is based on explicit and nonlinear parametric foreground models, the Needlet Internal Linear Combination (NILC; Delabrouille et al. 2009) method takes a fundamentally different approach, and assumes simply that the total foreground signal may be written as a linear combination of the true CMB signal and the total foreground signal at each frequency channel, and then form the linear combination of frequency maps that minimizes the variance of the final product, see Appendix B.1. To account for possible spatial variations in the foreground SED, the linear combination weights are computed separately in needlet space, allowing for localized optimization both in terms of sky position and angular scale. This method has been applied to a wide range of data sets, including Planck (e.g., Planck Collaboration IX 2016) and LiteBIRD. For further details, see Appendix B.2.
Finally, we consider the more recently developed constrained moments ILC method (cMILC; Remazeilles et al. 2021), which combines the momentexpansion and NILC methods. Specifically, additional constraints are imposed on the ILC weights that explicitly cancels out individual foreground components through their Taylorexpanded SEDs. Each moment corresponds to one additional linear constraint in the ILC solution, and the combination of foregroundspecific and the variance constraints is solved through a single Lagrange multiplier system. For further information, see Appendix B.3.
5. Ideal model analysis: Estimation of statistical uncertainties
When contemplating a major modification of a satellite’s baseline design, there are at least two aspects that must be considered carefully. The first aspect is sensitivity – we must investigate how much stronger (or weaker) constraints on the tensortoscalar ratio the proposed modification will lead to. The second aspect is goodnessoffit – we must also investigate whether the proposed modification will affect our sensitivity to modeling errors or our ability to detect such errors. These two main aspects are addressed individually in this and the next section.
Starting with the sensitivity aspect, we approach this by analyzing the simulations summarized in Sect. 3 for each of the instrument configurations described in Sect. 2. In this section, we primarily focus on the simplest foreground model, namely the onecomponent MBB thermal dust model, for which modeling errors are minor (it is worth noting that they are by no means nonexistent, since both informative priors and spatial variations in the fitted parameters can lead to biases).
We also note that each of the five componentseparation methods discussed in Sect. 4 have their particular strengths and weaknesses, and our goal in this paper is not to perform a headtohead componentseparation algorithm comparison, but it is rather to assess the capabilities and limitations of the LiteBIRD frequency selection itself. In the following, we therefore present selected results from among the five methods, depending on which product is most convenient for a particular application. In general, however, we note that all five methods provide qualitatively similar results.
5.1. Map residuals
First, to build intuition regarding the overall impact of the frequency range, we show in Fig. 3 the reconstructed Bmode NILC^{1} map for one arbitrarily selected simulation and each of the two most extreme instrument configurations, M0 (top row) and M5 (bottom row). While residual foreground contamination is clearly visible in the Galactic disk region of the CMB Bmode map for M0, we see that these residuals are significantly reduced for the M5 configuration, with a wider frequency range.
Fig. 3. Comparison of the reconstructed CMB Bmode maps from NILC for the M0 (top row) and M5 (bottom row) instrument configurations. The two panels in each row show opposite hemispheres aligned with the Galactic plane along the equator, and the north and south Galactic poles at the top and bottom, respectively. The left and right panels are centered on the Galactic center and anticenter. 
A complementary view from the parametric Commander code is shown in Fig. 4 in the form of residual (outputminusinput) maps for four key quantities, namely the Stokes Q and U CMB maps and the thermal dust spectral index and temperature for the M0, M2, and M5 instrument configurations^{2}. In all cases, we see that the quality of the fit improves with increasing maximum frequency. In the CMB maps, this is most clearly seen in the form of reduced scatter around the Galactic plane, while for T_{d} the dark blue area in the high Galactic latitudes systematically fade from M0 to M5; correspondingly, the dark red area in β_{d} also gradually fade away. These observations may be understood intuitively by noting that an increased frequency range improves the ability to break the wellknown degeneracy between β_{d} and T_{d} (e.g., Juvela & Ysard 2012) that arises when fitting for them jointly.
Fig. 4. Commander residual (output minus input) maps for selected parameters and instrument configurations for one arbitrarily selected Monte Carlo sample. Columns show, from left to right: (1) CMB Stokes Q; (2) CMB Stokes U; (3) thermal dust spectral index; and (4) thermal dust temperature. Rows show M0, M2, and M5. The maps show one typical realization. Gray pixels indicate the masked pixels with f_{sky} = 90%. 
5.2. Power spectrum uncertainties
Next, we consider the impact of changing the instrument configuration from M0 to M5 in terms of residual foreground and noise angular power spectra. This is shown in Fig. 5 for both NILC (solid lines) and cMILC (dashed lines). The left panel shows the projected foreground power spectrum, while the right panel shows the noise sample variance. First, we see that cMILC provides a significantly lower residual foreground contamination than NILC across a wide range of multipoles for both instrument configurations. This is due to the additional deprojection of foreground moments, which is more efficient at suppressing true foreground residuals. However, this additional foreground reduction performance does come at a significant price in the form of higher noise uncertainties caused by the same additional degrees of freedom. This noise penalty is greatly reduced when extending the LiteBIRD frequency range from 402 to 603 GHz, to the point that it is comparable to the CMB cosmic variance. In addition, the overall level of residual foreground contamination is also greatly reduced for the extended frequency range, both for NILC and cMILC, across all multipoles.
Fig. 5. Power spectrum of residual foreground contamination (left panel) and noise (right panel) with NILC (solid lines) and cMILC (dashed lines) on f_{sky} = 50% of the sky for the M0 (blue) and M5 (orange) configurations. For reference, we show the input CMB Bmode signal (r = 0, black dashdotted line), the primordial Bmode signal expected from theory for values ranging from r = 10^{−2} down to r = 10^{−3} (grayshaded area), and the residual lensing cosmic variance (yellowshaded area) for no delensing (A_{L} = 1) up to 90% delensing (A_{L} = 0.1). All the spectra are binned with a multipole bin size of Δℓ = 16. 
5.3. Tensortoscalar ratio constraints
We now turn our attention to tensortoscalar ratio constraints, and we start with probability distributions as derived using the NILC and cMILC methods. These are summarized in Fig. 6^{3}. These distributions are derived by averaging over an ensemble of simulations, and therefore correspond inherently to ensemble averages. The top and bottom panels show the results for the M0 and M5 configurations, respectively. Solid and dotted lines show recovered tensortoscalar ratio distributions for NILC and cMILC, while red and black lines show results with and without delensing.
Fig. 6. Comparison of recovered tensortoscalar ratio distributions from NILC (solid lines) and cMILC (dotted lines) for the M0 (top panel) and M5 (bottom panel) instrument configurations. Results without delensing are shown as black curves, while results with full delensing are shown as red lines. 
Since the fiducial value of the tensortoscalar ratio in these simulations is r = 0, any apparent bias (r ≠ 0, defined by the peak of the likelihood) is due to the power spectrum of the residual foreground contamination that projects into the CMB Bmode maps shown in Fig. 3. In contrast, the width of the distribution, that is σ_{r}, receives contributions from the cosmic variance of the lensed CMB signal, the sample variance of projected foregrounds, and noise.
Starting with the unlensed NILC results (solid black lines), we first note that both M0 and M5 have a welldefined nonzero peak, which is indicative of a nonnegligible foreground residual, consistent with the nonnegligible foreground power seen in the left panel of Fig. 5. For the delensed case (red curve), this spurious detection is statistically significant at the 5σ level for M0, while for the lensed case it is nearly consistent with zero at the 2σ level; however, this difference is simply due to the additional uncertainty added by lensing, and not lower foreground residuals as such. For M5, the overall NILC bias is reduced by about 30% in both cases, while the uncertainty remains unchanged, since this is dominated by cosmic variance.
Being a more constrained version of the ILC, the semiblind cMILC method allows further foreground deprojection, and thereby suppresses some of the residual biases that NILC suffers from. Specifically, by assuming an MBB with pivots and a powerlaw with pivot as the baseline zerothorder SEDs for dust and synchrotron, cMILC imposes four nulling constraints to deproject the zerothorder moments of dust and synchrotron and the firstorder moments of dust, whose respective SEDs are f_{sync}(ν), f_{dust}(ν), , and .
The recovered distributions of the tensortoscalar ratio from cMILC are shown as dotted lines in Fig. 6. In contrast to NILC, cMILC shows unbiased recovery of r = 0 thanks to deprojection of foreground moments. However, for M0 the extra cMILC constraints also significantly increase the noise contribution to σ_{r}. On the other hand, for M5 cMILC provides unbiased recovery of r = 0, with a negligible increase in noise. Thus, while M0 does not provide enough constraining power to allow deprojection of thermal dust temperature moments at a useful level, M5 does support this.
All the above results correspond to ensembleaveraged distributions derived with blind internal linear combination methods. To understand the behavior for individual realizations, Fig. 7 shows corresponding posterior distributions as derived with Commander for 20 realizations, plotted as thin lines for each of the six instrument configurations. No delensing is applied in any of these cases, and a uniform prior on r > 0 is assumed. The red vertical lines show the mean of the upper 68% confidence limits. Firstly, we note that these distributions appear (at least visually) statistically consistent with the input value of r = 0. Secondly, we also see that the upper limit decreases monotonically with instrument configuration, with M5 providing the tightest overall constraints. Thirdly, the internal scatter between individual realizations is notably smaller with a higher maximum frequency, resulting in an overall lower sample variance.
Fig. 7. Tensortoscalar ratio posterior distributions for the single MBB component model (for both simulating and modeling) as derived with Commander using multipoles ℓ = [2, 12] and a sky fraction of f_{sky} = 73% for 20 independent simulations (black lines), each drawn from a ΛCDM spectrum with r = 0. The red vertical lines correspond to the mean upper 68% confidence limit, σ_{r}. 
These general observations are summarized more quantitatively in Fig. 8, which shows σ_{r} as a function of instrument configuration, as derived using the different componentseparation methods. We note that σ_{r} for NILC is dominated by the cosmic variance of the bias, so this method is omitted from the figure. We note again that each of the methods involve different masks, different multipole ranges, and different foreground models, and the relative absolute shift between the curves is therefore entirely expected. Furthermore, we also note that the baseline configuration does reach the target sensitivity of σ_{r} < 0.001 for both the FGBuster and Commander methods, as already reported by LiteBIRD Collaboration (2022); even lower absolute uncertainties can be achieved by exploiting additional sky coverage (which is particularly relevant for FGBuster) or multipole range (which is particularly relevant for Commander).
Fig. 8. Tensortoscalar ratio uncertainty, σ_{r}, as a function of instrument configuration for the component separation methods considered in this paper. 
More importantly for the main topic of this paper, however, is the relative behavior between different configurations when we shift the entire HFT toward higher frequencies, and this is qualitatively similar for all methods: the uncertainty on r decreases monotonically for all methods over the frequency range considered in this paper.
In Fig. 9 we plot the ratio for different sky fractions and analysis methods. Overall, we see that M5 typically has 30–50% smaller uncertainties than M0, and the gains are larger when more sky is included in the analysis. This makes intuitive sense, since accurate componentseparation is relatively more important near the Galactic plane.
Fig. 9. Ratio between the tensortoscalar ratio uncertainties for M5 and M0 for different analysis choices and componentseparation methods. The black points are Commander analyses using different masks and sky fractions. 
5.4. Excluding the highest frequency channel
As seen in the above calculations, extending the LiteBIRD frequency range from 402 to 603 GHz leads to a substantially smaller σ_{r}. However, as discussed in greater detail in the next section, an extended frequency range is also associated with a potentially larger bias from modeling errors, since the overall component separation process becomes more dependent on assumptions regarding the thermal dust SED, both in terms of its specific parametric shape and its general spatial coherence between different frequency channels.
Recognizing these challenges, it is interesting to investigate how much the uncertainty on r degrades if we exclude the highest frequency channel. This represents the scenario in which we are unable to model thermal dust at frequencies above 500 GHz. In this exercise, we adopt the M2 configuration, with the highest frequency channel at 482 GHz, as a worked example. When excluding the highest frequency channel in this particular configuration, we recover an instrument model that is very similar in frequency coverage to the baseline, M0, except with slightly lower sensitivity in the remaining HFT channels.
Analyzing this truncated data set with Commander in the same manner as those above, we find σ_{r} = 9.8 × 10^{−4}, which is to be compared with σ_{r} = 9.2 × 10^{−4} for the nominal M0 configuration with otherwise identical settings, or an increase of about 5%. The potential loss in sensitivity is thus relatively modest for this scenario. However, it is important to note that an additional effect of this modification is that there is no frequency overlap between the MFT and HFT in the extended M2 configuration, while there is such overlap in the baseline M0 configuration. Effectively, the extended instrument configurations thus trade off control of unknown systematics in the overlap between the MFT and HFT at the 195 GHz band for additional sensitivity to thermal dust emission at the highest frequencies.
6. Nonideal model analysis: Sensitivity to modeling errors
We now turn our attention to the issue of modeling errors, and we consider two different types. The first type consists of SED modeling errors, in which the fitted SED model is significantly different from the SED of the true sky signal. Being able to reject a spurious foregroundinduced detection is key for a robust instrument design. The second type arises from oversmoothing of spatially varying SED parameters, which also can introduce biases.
6.1. Detectability of SED model errors
To assess the ability to detect SED modeling errors as a function of instrument configuration, we employ the parametric Commander method, which provides various goodnessoffit quantities in terms of residual maps and χ^{2} statistics among its output products. Specifically, we analyze simulations based on the complex SiFeC model summarized in Sect. 3, but using a simple onecomponent MBB dust model for fitting. For each instrument configuration and simulation, we compute both the tensortoscalar ratio posterior distributions, as in Sect. 5.3, and the maplevel reduced normalized χ^{2} statistic as defined by
In this expression, the sum runs over all unmasked pixels, p, and all frequencies, ν, and n_{d.o.f.} is the total number of degrees of freedom per pixel. Since the χ^{2} distribution converges to a Gaussian with mean n_{d.o.f.} and variance 2n_{d.o.f.} for large n_{d.o.f.}, this quantity measures statistical outliers in units of standard deviations. In this section, both the posterior distributions and the χ^{2} statistics are evaluated on the unmasked pixels encompassing f_{sky} = 60%. However, we note that the choice of sky fraction is arbitrary, as the main topic of this section is the impact of modeling errors; a smaller sky fraction will necessarily lead to smaller biases, but larger uncertainties, and vice versa. The choice adopted here is simply chosen as a representative compromise between minimizing both bias and uncertainties.
To build intuition regarding the following results, we show in the top section of Fig. 10 residual maps (r_{ν} = d_{ν} − s_{ν}) for one arbitrarily selected realization and four frequency channels for the ideal onecomponent MBB model. In this case we see no significant residuals in any frequency channels or instrument configurations, simply because the fitted SED model matches the actually simulated sky. All these maps are thus consistent with white noise, and the corresponding χ^{2} statistic is also consistent with the Gaussian noise expectation.
Fig. 10. Commander residual maps (r_{ν} = d_{ν} − s_{ν}) for four different frequency channels and three different instrument configurations. The top section show residual maps for the ideal onecomponent MBB model, while the bottom section show results for the nonideal SiFeC model, where a model mismatch was included. 
In the bottom section of Fig. 10, we show similar residual maps for the complex SiFeC simulations. Here we first see that the baseline M0 configuration appears visually consistent with the ideal MBB simulation in the top panel, and the parameter estimation algorithm will clearly not be able to identify the intrinsic model mismatch. Extending the frequency range to 482 GHz (M2), however, leaves noticeable residuals tracing high dust emission regions in the Galactic plane, and at 603 GHz (M5) the model mismatch is obvious.
The effect of these mismatches in terms of the tensortoscalar ratio posterior distributions is quantified in Fig. 11. Each panel shows 20 realizations as black lines, while the red line shows the median of the ideal onecomponent MBB simulations as a reference; since the overall foreground levels of the two sky models are comparable, and both the fitted model and the instrumental noise parameters are identical between the SiFeC and MBB models, the corresponding posterior distributions should also be roughly equal.
Fig. 11. Tensortoscalar ratio posterior distributions with a fiducial tensortoscalar ratio of r = 0, as derived with Commander when fitting SiFeC simulations using the onecomponent MBB model for six different instrument configurations. A χ^{2} statistic is reported for each configuration, indicating the maplevel goodness of fit. All results are derived from a sky fraction of f_{sky} = 60%. The red dashed lines correspond to the median distribution from the ideal onecomponent MBB fit to an MBB simulation shown in Fig. 7 and serves as a reference. 
Starting with the baseline configuration in the top panel, we see that the tensortoscalar ratio is biased high by a factor of about two (as seen in the shift between the black and red lines), even though each distribution is individually consistent with zero. The χ^{2} statistic is completely consistent with the Gaussian hypothesis, with a deviation of only −1.2σ. The M1 configuration behaves very similarly.
For the case M2, we see that the bias in r is greatly increased, to the point that this configuration would have reported more than a 3σ detection of nonzero Bmode power. At this point, the χ^{2} value has also increased somewhat, but is still only 2σ high compared to the Gaussian expectation. It is only with frequency ranges equal to or higher than 523 GHz that χ^{2} is sufficiently deviant that we are able to conclusively detect the model mismatch, with a statistical significance of 9σ or more.
While the combination of a significant bias of r and the statistical acceptance of M2 appear disconcerting, it is important to recall that these results are obtained with a relatively conservative sky fraction of 60%. For an analysis of real data, one would also exploit the information in the high dust emission regions to identify a statistically acceptable parametric model, and then fit that also at high latitudes. As seen in the residual maps in Fig. 10, such an analysis would indicate that the onecomponent MBB is sufficient for M0, while for M2 and higher configurations the residuals are an obvious red flag. From a model selection point of view, frequencies higher than ≈500 GHz are critically important for identifying poor thermal dust models.
6.2. Resolution of spatially varying SED parameters
Next, we consider the impact of different smoothing scales for SED spectral parameters. In this case, we analyze the single MBB simulations using the FGBuster algorithm, which supports tuned pixel sizes for each spectral parameter as part of its basic implementation, quantified in terms of the HEALPix resolution parameter, N_{side}. To disentangle the SED error effects from the spatial resolution effect, we once again consider the ideal onecomponent MBB simulations, but fit coarsely pixelized spectral index maps, while the true spectral index maps are smoothly varying on the sky. In particular, we focus on the thermal dust temperature, which is key for the highest LiteBIRD frequencies, and we allow this to take on different resolutions at low, intermediate, and high Galactic latitudes, denoted ; see Appendix A.2 for details.
We are now interested in quantifying the recovery accuracy of r as a function of , and the results from this analysis are summarized in Fig. 12 in terms of the recovered mean tensortoscalar ratio r (black crosses) and its uncertainty (bars). Groups along the horizontal axis indicate different combinations, while colors indicate different LiteBIRD configurations.
Fig. 12. Fisher uncertainty (bars) and bias (black crosses) on r for different instrumental (colors) and FGBuster componentseparation configuration (bar groups). The values on the horizontal axis define the N_{side} resolution at which the temperature of the thermal dust component is fitted for in the best 20%, 20–40% and 40–60% sky fractions. 
The first thing to note is that when the HFT channels are shifted to higher frequencies, the uncertainty on the spectral indices decreases, and this leads to a decrease of the statistical residuals, and therefore also to an improvement in σ_{r}; this is the same effect as was illustrated in Fig. 8, but now we also see that this holds true almost independently of .
At the same time, we see that low values of lead to a strong bias in r, and this effect increases rapidly with increasing frequency range. On the other hand, this bias can be mitigated by increasing at the expenses of an increase in the statistical uncertainty. In particular, for = [32,32,32], all LiteBIRD configurations lead to a bias 𝒪(3 × 10^{−4}), with an uncertainty σ_{r} ≈ 𝒪(6 × 10^{−4}). In general, a maximum central frequency above 500 GHz requires at least = 16, even for the cleanest 20% of the sky in order to keep the bias under control.
7. Discussion and conclusions
One of the most important challenges for nextgeneration CMB Bmode experiments is astrophysical foregrounds, and in particular polarized thermal dust emission from our own Milky Way. Properly choosing the frequency range is thus one of the key decisions to be made for any given Bmode space mission or balloonborne experiment. In this paper, we have revisited this issue for LiteBIRD. Although the established baseline configuration has already been shown to achieve the main mission goals (LiteBIRD Collaboration 2022), it is still desirable to optimize the overall uncertainties to increase the overall system margins.
In this paper, we have addressed this issue by analyzing simulations to measure the recovered uncertainty on the tensortoscalar ratio using several different componentseparation algorithms for six different possible instrument configurations, where the entire HFT frequency range is shifted toward higher frequencies. We find that the total statistical uncertainty on the tensortoscalar ratio, r, after foreground cleaning may be reduced by 30–50% by extending the current maximum frequency from the current baseline of 402 GHz to 500 GHz, or higher for the single MBB dust model, also depending on the sky fraction.
A wider frequency range also increases the instrument’s sensitivity to the shape of the foreground SED, and this has two important consequences. On the one side, a wider frequency range imposes stronger requirements on the foreground model itself, since modeling errors may more easily contaminate the CMB results. For instance, if a given componentseparation algorithm (explicitly or implicitly) assumes that the highfrequency channels correlate perfectly with the thermal dust emission at lower frequencies, while the true sky exhibits frequency decorrelation due to the 3D structure of the Milky Way, some fraction of the difference between the assumed and real sky model will leak into the CMB (Tassis & Pavlidou 2015). It is therefore essential that the method of choice allows for sufficient flexibility to capture such uncertainties and model mismatches.
At the same time, a wider frequency range also provides the tools for actually identifying and constraining that model in the first place. Indeed, one of the most important concerns regarding a limited frequency range is the possibility of having lowlevel foreground residuals that mildly bias cosmological parameters, but still result in an acceptable goodnessoffit, (for instance as measured by a χ^{2} statistic). Higher frequency channels provide important safeguards against this type of error. Additionally, having access to higher frequency information increases the number of basis functions that can be used to model the thermal dust emission as a function of frequencies. A particularly important example of such a mode is the spatial distribution of the thermal dust temperature, which is useful for both blind and nonblind methods.
An important intuition that underlies all the results discussed in this paper is simply the fact that the CMB SED falls nearly exponentially (relative to dust) above 300 GHz. This implies that even relatively small modifications of the frequency range can lead to significant improvements in a given instrument’s ability to separate CMB from thermal dust. In particular, at the sensitivity level of LiteBIRD, a 402 GHz channel has a nonnegligible signaltonoise ratio for CMB fluctuations, while at 500 GHz this is essentially zero. For an extended configuration, the highest frequency channel therefore represents essentially a pure thermal dust map, while for the baseline configuration it represents a weighted sum of CMB and thermal dust emission. As shown quantitatively in this paper, a clean foreground tracer substantially improves tensortoscalar constraints for both parametric and nonparametric (ILC) componentseparation methods.
In this paper we have primarily addressed the impact of Galactic foregrounds from the point of view of CMB confusion. An important lesson learned from both WMAP and Planck, however, is that Galactic foregrounds also play a key role in understanding and mitigating lowlevel instrumental systematic effects. For instance, confusion between transmission imbalance factors and foregrounds turned out to be an important uncertainty for WMAP (Page et al. 2007; Watts et al. 2023), while confusion between gain calibration and foregrounds was a key limitation for Planck LFI (Planck Collaboration I 2020; BeyondPlanck Collaboration 2023; Gjerløw et al. 2023). One may therefore argue that experience shows that the most robust approach to mitigating the full effect of astrophysical foregrounds is not by avoiding them, but rather by measuring them. In addition to improving CMB constraints, a wide frequency range will also increase the amount of Galactic science that may be derived from these observations, and it will thereby also increase the overall legacy value of the mission for other fields of astronomy.
While the various instrumental setups studied in this paper appear as possible options, we would like to stress that their technical feasibility has not been performed in depth yet. This would have to include both aspects of optical design and detection chains, but also the potential impacts on the needs for optical modeling and calibration facilities. This tradeoff study would also have to deal with the impact of such modifications of the HFT design on the overall systematics effects, which should be propagated into systematics uncertainties on r. As noted in the introduction, performing such a full tradeoff study lies far beyond the scope of this paper. Rather, the goal of this paper is to understand whether sufficiently significant gains are achievable under ideal conditions to motivate a proper detailed study – and with typical gains at the 30–50% level in the statistical uncertainty on r after foreground cleaning, the results do appear interesting enough to be investigated further.
This large number of samples is the main motivation for using Commander1 rather than Commander3; producing 𝒪(10^{6}) posterior samples at full angular resolution would require tens of millions of CPU hours. While feasible for a final production analysis of real data, such a cost is not justified for the current exploratory analysis.
Acknowledgments
This work is supported in Japan by ISAS/JAXA for PrePhase A2 studies, by the acceleration program of JAXA research and development directorate, by the World Premier International Research Center Initiative (WPI) of MEXT, by the JSPS CoretoCore Program of A. Advanced Research Networks, and by JSPS KAKENHI Grant Numbers JP15H05891, JP17H01115, and JP17H01125. The Canadian contribution is supported by the Canadian Space Agency. The French LiteBIRD phase A contribution is supported by the Centre National d’Études Spatiales (CNES), by the Centre National de la Recherche Scientifique (CNRS), and by the Commissariat à l’Énergie Atomique (CEA). The German participation in LiteBIRD is supported in part by the Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy (Grant No. EXC2094 – 390783311). The Italian LiteBIRD phase A contribution is supported by the Italian Space Agency (ASI Grants No. 20209HH.0 and 201624H.12018), the National Institute for Nuclear Physics (INFN) and the National Institute for Astrophysics (INAF). Norwegian participation in LiteBIRD is supported by the Research Council of Norway (Grant No. 263011) and has received funding from the European Research Council (ERC) under the Horizon 2020 Research and Innovation Programme (Grant agreement No. 772253 and 819478). The Spanish LiteBIRD phase A contribution is supported by the Spanish Agencia Estatal de Investigación (AEI), project refs. PID2019110610RBC21, PID2020120514GBI00, ProID2020010108 and ICTP20210008. Funds that support contributions from Sweden come from the Swedish National Space Agency (SNSA/Rymdstyrelsen) and the Swedish Research Council (Reg. no. 201903959). The US contribution is supported by NASA grant no. 80NSSC18K0132. We also acknowledge funding from the European Research Council (ERC) under the Horizon 2020 Research and Innovation Programme (Grant agreement No. 725456 and 849169) and The Royal Society (Grant No. URF/R/191023)
References
 Ade, P. A. R., Ahmed, Z., Aikin, R. W., et al. 2016, Phys. Rev. Lett., 116, 031302 [NASA ADS] [CrossRef] [Google Scholar]
 Albrecht, A., & Steinhardt, P. J. 1982, Phys. Rev. Lett., 48, 1220 [Google Scholar]
 Aurlien, R., Remazeilles, M., Belkner, S., et al. 2023, J. Cosmol. Astropart. Phys., 2023, 034 [CrossRef] [Google Scholar]
 Azzoni, S., Abitbol, M. H., Alonso, D., et al. 2021, J. Cosmol. Astropart. Phys., 2021, 047 [CrossRef] [Google Scholar]
 Bennett, C. L., Hill, R. S., Hinshaw, G., et al. 2003, ApJS, 148, 97 [Google Scholar]
 BeyondPlanck Collaboration (Andersen, K. J., et al.) 2023, A&A, 675, A1 [Google Scholar]
 Chluba, J., Hill, J. C., & Abitbol, M. H. 2017, MNRAS, 472, 1195 [NASA ADS] [CrossRef] [Google Scholar]
 Delabrouille, J., Cardoso, J.F., Le Jeune, M., et al. 2009, A&A, 493, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eriksen, H. K., O’Dwyer, I. J., Jewell, J. B., et al. 2004, ApJS, 155, 227 [Google Scholar]
 Eriksen, H. K., Jewell, J. B., Dickinson, C., et al. 2008, ApJ, 676, 10 [Google Scholar]
 Errard, J., & Stompor, R. 2019, Phys. Rev. D, 99, 069907 [NASA ADS] [CrossRef] [Google Scholar]
 Errard, J., Stivoli, F., & Stompor, R. 2011, Phys. Rev. D, 84, 043529 [NASA ADS] [CrossRef] [Google Scholar]
 Fuskeland, U., Wehus, I. K., Eriksen, H. K., & Næss, S. K. 2014, ApJ, 790, 104 [Google Scholar]
 Fuskeland, U., Andersen, K. J., Aurlien, R., et al. 2021, A&A, 646, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gjerløw, E., Ihle, H. T., Galeotta, S., et al. 2023, A&A, 675, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [Google Scholar]
 Guilloux, F., Faÿ, G., & Cardoso, J.F. 2009, Appl. Comput. Harmon. Anal., 26, 143 [CrossRef] [Google Scholar]
 Guth, A. H. 1981, Phys. Rev. D, 23, 347 [Google Scholar]
 Hazumi, M., Borrill, J., Chinone, Y., et al. 2012, Proc. SPIE Int. Soc. Opt. Eng., 8442, 844219 [NASA ADS] [Google Scholar]
 Hazumi, M., Ade, P. A. R., Adler, A., et al. 2020, SPIE Conf. Ser., 11443, 114432F [Google Scholar]
 Hensley, B. S., & Draine, B. T. 2017, ApJ, 834, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Hensley, B. S., & Draine, B. T. 2023, ApJ, 948, 55 [Google Scholar]
 Hoang, T., & Lazarian, A. 2016, ApJ, 831, 159 [Google Scholar]
 Juvela, M., & Ysard, N. 2012, A&A, 539, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. D, 55, 7368 [NASA ADS] [CrossRef] [Google Scholar]
 Kogut, A., Dunkley, J., Bennett, C. L., et al. 2007, ApJ, 665, 355 [Google Scholar]
 Krachmalnicoff, N., Baccigalupi, C., Aumont, J., Bersanelli, M., & Mennella, A. 2016, A&A, 588, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krachmalnicoff, N., Carretti, E., Baccigalupi, C., et al. 2018, A&A, 618, A166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lawson, K. D., Mayer, C. J., Osborne, J. L., & Parkinson, M. L. 1987, MNRAS, 225, 307 [Google Scholar]
 Liebe, H. J. 1981, Radio Science, 16, 1183 [Google Scholar]
 Linde, A. D. 1982, Phys. Lett. B, 108, 389 [Google Scholar]
 Linde, A. D. 1983, Phys. Lett. B, 129, 177 [Google Scholar]
 LiteBIRD Collaboration 2022, Prog. Theor. Exp. Phys., 2023, 042F01 [Google Scholar]
 Mangilli, A., Aumont, J., Rotti, A., et al. 2021, A&A, 647, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matsumura, T., Akiba, Y., Borrill, J., et al. 2014, J. Low Temp. Phys., 176, 733 [Google Scholar]
 Montier, L., Mot, B., de Bernardis, P., et al. 2020, SPIE Conf. Ser., 11443, 114432G [NASA ADS] [Google Scholar]
 Narcowich, F., Petrushev, P., & Ward, J. 2006, SIAM J. Math. Anal., 38, 574 [CrossRef] [Google Scholar]
 Novotny, V., & Meincke, P. 1975, J. Low Temp. Phys., 18, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJS, 170, 335 [Google Scholar]
 Pardo, J., Cernicharo, J., & Serabyn, E. 2001, IEEE Trans. Antennas Propag., 49, 1683 [Google Scholar]
 Pelgrims, V., Clark, S. E., Hensley, B. S., et al. 2021, A&A, 647, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2014, A&A, 571, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2016, A&A, 594, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXV. 2016, A&A, 594, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2020, A&A, 641, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2020, A&A, 641, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. LVII. 2020, A&A, 643, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Platania, P., Bensadoun, M., Bersanelli, M., et al. 1998, ApJ, 505, 473 [NASA ADS] [CrossRef] [Google Scholar]
 Puglisi, G., Mihaylov, G., Panopoulou, G. V., et al. 2022, MNRAS, 511, 2052 [NASA ADS] [CrossRef] [Google Scholar]
 Remazeilles, M., Delabrouille, J., & Cardoso, J.F. 2011, MNRAS, 410, 2481 [Google Scholar]
 Remazeilles, M., Dickinson, C., Eriksen, H. K. K., & Wehus, I. K. 2016, MNRAS, 458, 2032 [Google Scholar]
 Remazeilles, M., Rotti, A., & Chluba, J. 2021, MNRAS, 503, 2478 [NASA ADS] [CrossRef] [Google Scholar]
 Ritacco, A., Boulanger, F., Guillet, V., et al. 2023, A&A, 670, A163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sato, K. 1981, MNRAS, 195, 467 [Google Scholar]
 Sen, A., Basak, S., Ghosh, T., Adak, D., & Sinha, S. 2022, arXiv eprints [arXiv:2212.02869] [Google Scholar]
 Starobinsky, A. A. 1980, Phys. Lett. B, 91, 99 [Google Scholar]
 Stompor, R., Leach, S., Stivoli, F., & Baccigalupi, C. 2009, MNRAS, 392, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Tassis, K., & Pavlidou, V. 2015, MNRAS, 451, L90 [NASA ADS] [CrossRef] [Google Scholar]
 Thorne, B., Dunkley, J., Alonso, D., & Næss, S. 2017, MNRAS, 469, 2821 [Google Scholar]
 Vacher, L., Aumont, J., Montier, L., et al. 2022, A&A, 660, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vacher, L., Aumont, J., Boulanger, F., et al. 2023a, A&A, 672, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vacher, L., Chluba, J., Aumont, J., Rotti, A., & Montier, L. 2023b, A&A, 669, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Watts, D. J., Galloway, M., Ihle, H. T., et al. 2023, A&A, 675, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Parametric componentseparation methods
A.1. Commander
Commander is a standard parametric Bayesian componentseparation framework for CMB observations (Eriksen et al. 2004, 2008) that has been used extensively by Planck (e.g., Planck Collaboration XII 2014; Planck Collaboration IX 2016; Planck Collaboration IV 2020; Planck Collaboration Int. LVII 2020). In the current analysis we employ the first implementation of this framework, typically referred to as Commander1, which was used in the Planck 2013 data release. The main advantage of this framework is a very high computational efficiency, which is important when analyzing an ensemble of simulations; at the same time, the main drawback is its requirement of identical angular resolution of all frequency bands, which limits its use to lowresolution data sets. In our current analysis, all input maps are at N_{side} = 16, with a common resolution of 10° FWHM. Another option would be to use Commander3 (BeyondPlanck Collaboration 2023), which supports both multiresolution component separation and lowlevel timeordered data analysis, but this would increase the computational expense of the analysis by orders of magnitude, without actually providing valuable new information regarding the central question of frequency range considered in this paper.
The starting point of any parametric Bayesian analysis is to define an explicit parametric data model, described in Eqs. (1)– (5), and in Commander1, every parameter is fitted per pixel. The main goal is now to fit the parameters θ = {a_{CMB}, a_{s}, a_{d}, β_{s}, β_{d}, T_{d}}, which within a Bayesian framework means computing the posterior distribution P(θd), and then using this to estimate the cosmological parameter r. Using Bayes’ theorem, the posterior distribution may be written as
where P(d) is a normalization factor, P(dθ) is the likelihood, and P(θ) represents some set of priors. In this work, we apply loose Gaussian priors on the spectral parameters, but no priors on the amplitude parameters. Specifically, we adopt a synchrotron prior of P(β_{s}) = N(−3.0, 0.3), where N(μ, σ) denotes a standard Gaussian distribution with mean μ and standard deviation σ. For thermal dust emission, we adopt P(β_{d}) = N(1.54, 0.20) and P(T_{d}) = N(23 K, 7 K) where we fit for the one MBB component.
To map out the distribution in Eq. (A.2), we use Gibbs sampling as implemented in Commander. That is, rather than sampling directly from P(θd), which is computationally unfeasible, we iteratively sample from each of the corresponding conditional probability distributions. For the purposes of this paper, this can be written as
Here ← indicates that the parameters to the left are sampled from the distribution on the righthand side. All amplitude parameters are sampled jointly together with the CMB sky signal to avoid excessive Monte Carlo correlation lengths. The nonlinear parameters, the spectral indices and dust temperature, are sampled using a standard inversion sampler.
We run this iterative Gibbs sampler until convergence, which in the current analysis requires about 40 000 samples after rejecting a burnin period of about 500 samples. We repeat this analysis for 20 different CMB and noise realizations.^{4} Based on the samples from each simulation, we compute a posterior mean CMB map, , and the corresponding covariance matrix, N, and these are subsequently fed into a standard Gaussian likelihood;
Here S(r) is the CMB signal covariance matrix for a given value of r and the bestfit ΛCDM parameters described in Sect. 3. The low pixel resolution of N_{side} = 16 used for Commander in this paper is dictated by the CPU time and RAM required to evaluate this equation, as well as the number of samples needed to establish a robust noise covariance matrix.
We adopt a uniform prior for all positive values of r, and the posterior distribution is therefore numerically identical to the likelihood in Eq. (A.5) for r > 0. This function is mapped out by gridding Eq. (A.5) over a precomputed library of angular power spectra, C_{ℓ}(r), which are computed into covariance matrices in pixel space, S(r), including only multipoles in the range ℓ = [2, 12], as dictated by the low pixel resolution of N_{side} = 16. The Commander results thus only probe the lowℓ reionization peak of the Bmode power spectrum, not the recombination peak around ℓ ≈ 100. However, this is also where diffuse component separation is most challenging, and the question of frequency optimization is most important.
The main analysis mask used for the Commander analysis is the f_{sky} = 73% polarization mask discussed in Planck Collaboration IV (2020). To visualize the impact of sky fraction, other masks have been used in Fig. 9, such as masks based on high latitude cuts, χ^{2}, and by thresholding the 402 GHz polarization amplitude () map, which primarily suppresses thermal dust emission. This latter is used for the 60% mask in Sect. 6.
A.2. FGBuster
A second implementation of parametric component separation is called FGBuster.^{5} This method is based on the same parametric modeling as defined by Eqs. (1)–(5), but uses a computationally efficient nonlinear optimization procedure to perform the fit, as opposed to bruteforce Monte Carlo sampling as implemented in Commander. Following Stompor et al. (2009), Errard et al. (2011), Errard & Stompor (2019), and Puglisi et al. (2022), FGBuster proceeds in three steps. These are: first the estimation of the spectral parameters {β_{d}, T_{d}, β_{s}} through the optimization of a socalled spectral likelihood; second, the estimation of the component maps {a_{CMB}, a_{d}, a_{s}} introduced in Appendix A.1; and third, the estimation of the angular power spectra for the residual foregrounds (i.e., the difference between the input and recovered CMB map), as well as the estimation of the tensortoscalar ratio r with a Gaussian likelihood defined by the previous quantities, the theoretical BB spectrum, and the noise after component separation.
Regarding the first step, we note that FGBuster supports different resolutions (as defined by the HEALPix resolution parameter, N_{side}) for each spectral parameter, while the corresponding amplitudes are fitted at full angular resolution. Choosing a low N_{side} implies that few foreground parameters are fitted for, resulting in low postcomponentseparation noise; however, this may also result in a high bias due to sky complexity that the model is not allowed to capture. Optimizing N_{side} for each parameter is an important point of the algorithm, and we quantify these two competing effects by splitting the foreground residuals into two categories. First, the systematic residuals, as obtained from a noiseless simulation, are driven by the mismatch between input sky and the assumed parametric model, and its power spectrum, B_{ℓ}, defines the bias on r. Second, the statistical residuals are computed as the average power spectrum from 10 noisy simulations, S_{ℓ}, which are driven by the statistical error while estimating spectral parameters, and by coaddition of the noise at different frequencies.
We consider the following approximate likelihood for the power spectrum:
where f_{sky} = 0.49, D_{ℓ} = L_{ℓ} + S_{ℓ} + B_{ℓ}, C_{ℓ} = rT_{ℓ} + L_{ℓ} + S_{ℓ}, L_{ℓ} is the BB power spectrum from gravitational lensing, and T_{ℓ} is the tensor power spectrum for r = 1. The bias on r is estimated by maximizing this likelihood, while we define σ_{r} to be the Fisher estimate obtained by assuming B_{ℓ} = 0. We prefer this choice over computing σ_{r} from the width of Eq. (A.6), because the latter procedure produces values of σ_{r} that are strongly dependent on the bias, while here we want it to reflect only the statistical constraining power of a given configuration. The various spectra involved in the FGBuster analysis are illustrated in Fig. A.1.
Fig. A.1. Example of the results from the FGBuster componentseparation runs. Two instrumental configurations are considered: the baseline M0 (blue) and the high frequency extension M3 (red). The solid lines represent the postcomponentseparation noise and foreground residuals S_{ℓ}, computed by averaging 1000 simulations for the baseline, and averaging 10 simulations for the extension. Together with the power from gravitational lensing and primordial Bmodes, they define the error bars on the CMB power spectrum (boxes). 
The free parameters of our setup are the N_{side} values of the spectral parameters, and the LiteBIRD instrumental configuration. Regarding the former, we follow LiteBIRD Collaboration (2022), and split the sky into three (almost) isoGalacticlatitude regions, each covering approximately 20% of the sky. In each region, every parameter can have different values of N_{side}, tuned according to its local signaltonoise ratio. We define the three regions using Planck postprocessing masks adopted for temperature component separation (Planck Collaboration IV 2020), the socalled GAL20, GAL40, and GAL60. We find this (arbitrary) choice to significantly help coping with the varying signaltonoise, but the use of more of these separate regions is certainly possible and could further improve the postcomponent separation performance – especially if specialized for dust polarization. Regarding the choice of the N_{side} values of the spectral parameter, we find that an accurate estimation of β_{d} is particularly essential for the LiteBIRD configuration, due to a small statistical uncertainty but potentially large biases coming from spatial variability; we use in all the regions. Less accurate estimation of β_{s} and T_{d} are usually tolerable as was experienced with foregroundonly simulations, and following the preliminary results from Errard & Stompor (2019). For synchrotron emission, we adopt for the three regions (counting from lower latitudes), and this is kept constant for all instrumental configurations; we do not expect changes at the high frequencies to significantly affect the fit of the lowfrequency foreground.^{6} In contrast, increasing the HFT frequencies does increase the constraining power for the dust SED, and therefore requires a more precise modeling of the spatial variability of T_{d} to keep the systematic residuals under control. The analysis sky mask covers f_{sky} = 49% of the sky as introduced in LiteBIRD Collaboration (2022). In addition to the Galactic plane, this mask removes regions with high residuals.
A.3. Moment expansion
A third parametric componentseparation method is referred to as the “momentexpansion method,” which models the complexity of the polarized dust SED using parametric fitting with moment expansion coefficients in harmonic space. This formalism was introduced by Chluba et al. (2017) and generalized to polarization in Vacher et al. (2023b,a) and at the crossfrequency power spectra level in Mangilli et al. (2021). An application of this formalism for component separation with the LiteBIRD mission at high frequencies can be found in Vacher et al. (2022). The formalism and the pipeline developed in this last study will be applied identically in the present work.
In true experimental conditions, averaging over different SEDs emitted from distinct regions with different physical properties is unavoidable. They appear along the line of sight in the 3D sky, between different lines of sight, inside the beam of the instrument, or when performing a spherical harmonic decomposition to calculate the angular power spectra over large regions of the sky. Since SEDs depend of the frequency in a nonlinear fashion, summing over them when spectral parameters are varying will lead to distortions that should be understood and modeled. The momentexpansion method aims at accounting for SED distortions from a canonical model by introducing new terms coming from a Taylorinspired expansion of the SED around a pivot value and with respect to its spectral parameters (e.g., β_{d} and T_{d} for an MBB). One can generalize this expansion at the crossfrequency power spectra level 𝒟_{ℓ}(ν_{i} × ν_{j}), providing an analytical model that can be fitted over the foreground signal.
As discussed in Vacher et al. (2022), providing a correct model for the distortions due to spatial variations of the dust temperature is critical for a mission like LiteBIRD, while it could safely be ignored for missions like Simons Observatory (Azzoni et al. 2021), with different frequency coverage and lower sensitivity. If the dust signal in every Galactic pixel is indeed a modified blackbody, reaching higher frequencies will ease component separation by breaking degeneracies between the moments in β_{d} and T_{d}, providing a better characterization of the impact of the dust temperature variations on the shape of the power spectra SED.
For the present application, we extract the crossfrequency power spectra of N_{sim} = 500 simulations, including different realizations of Gaussian noise and CMB lensing, with all the highfrequency configurations from M0 to M5. The only foreground component considered is thermal dust with a single MBB. The simulation maps are produced at N_{side} = 512 and masked in order to keep a sky fraction of f_{sky} = 70%, using a raw mask made from Planck intensity data at 353 GHz as in Vacher et al. (2022). Only the nine highest frequencies above 100 GHz are considered, leading to 13 bands (counting the ones sharing the same frequency) and thus 91 crossspectra. The spectra are binned up to ℓ = 200 with bins of Δℓ = 10. Since we expect it to be the only fitting procedure allowing us to recover an unshifted r posterior for the standard LiteBIRD configuration, we shall only consider the rβ–T fitting scheme, using the moment coefficients up to order 1 in both β_{d} and T_{d} (Vacher et al. 2022). For each simulation, a χ^{2} minimization is done at the crossfrequency power spectra level
with the residual and the covariance matrix: , allowing us to extract one bestfit value for the tensortoscalar ratio per simulation. A histogram is built out of the N_{sim} best fit values of , which gives us the final posteriors, with associated standard deviation found by fitting a Gaussian curve over it.
Appendix B: Minimum variance componentseparation methods
B.1. Internal linear combinations
Generally speaking, internal linear combination (ILC) methods perform a weighted linear combination of the data d_{ν} across frequencies:
which is designed to minimize the overall variance due to foregrounds and noise without altering the CMB signal. More explicitly, the ILC weights w ≡ {w_{ν}} are the solution of a constrained minimization problem that can be formulated by the Lagrangian,
in which the matrix C collects the elements C_{νν′} = ⟨d_{ν} d_{ν′}⟩ of the data covariance matrix for all pairs of frequencies (ν, ν′) and the vector a ≡ {a_{ν}} collects the frequency spectrum of the CMB component across the channels. The Lagrange multiplier λ is used to ensure that the overall variance w^{T}C w of the ILC estimate is minimized under the constraint ∑_{ν}w_{ν} a_{ν} = 1 for the preservation of the CMB signal. The ILC weights are thus given by the saddle point of the Lagrangian defined in Eq. (B.2):
and as such they do not rely on any explicit parametric foreground model.
B.2. NILC
The Needlet ILC (NILC; Delabrouille et al. 2009) method is a specific version of the ILC approach that performs this optimization using a needlet (i.e., spherical wavelet) vector basis. In the NILC method, the ILC weights, Eq. (B.3), are computed independently in different regions of the sky and different ranges of angular scales using a wavelet decomposition of the data. Needlets have excellent localization properties, both on the sphere (pixel space) and in harmonic space (Narcowich et al. 2006; Guilloux et al. 2009), which allows the NILC algorithm to adjust the ILC weights depending on the local variations of the foreground and noise contaminations across the sky and across different scales.
B.3. cMILC
The Constrained Moment ILC (cMILC; Remazeilles et al. 2021) method is an extension of the NILC method that includes additional nulling constraints to deproject spectral moments of the foreground emission arising from lineofsight integration and beam convolution. It combines the socalled constrained ILC approach (Remazeilles et al. 2011) and the moment expansion technique for dealing with foreground emission (Chluba et al. 2017). Because of spectral variations of foreground emission along the lineofsight and/or inside the beam, the lineofsight integration and beam convolution effects in sky observations cause spectral distortions to the baseline SEDs of thermal dust and synchrotron, and thus induce partial decorrelation across frequencies. Such higherorder corrections to the baseline foreground emission can be significant compared to the CMB Bmode signal at low tensortoscalar ratio values.
These effects are addressed by cMILC through moment expansion of the foreground emission beyond the leadingorder SED:
where is the baseline SED evaluated at fixed pivot values for a set of spectral parameters (e.g., dust spectral index and temperature) that vary with the direction on the sky and along the line of sight. The expansion defined by Eq. (B.4) highlights the foreground moment components,
which have homogeneous SEDs given by the derivatives of the baseline SED with respect to spectral parameters:
where α = (α_{1}, ⋯, α_{n})∈ℕ^{n}.
cMILC enables us to deproject several moments of the foreground emission from the recovered CMB Bmode map by generalizing the Lagrangian in Eq. (B.2) as
in which extra Lagrange multipliers μ_{α} are introduced to impose several nulling constraints on the ILC weights, ∑_{ν}w_{ν} b_{α}(ν) = 0, with respect to the moment SEDs in Eq. (B.6). The cMILC weights are thus given by the saddle point of the extended Lagrangian in Eq. (B.7),
where the mixing matrix A = (a b_{1} ⋯ b_{m}) contains the CMB SED vector a in the first column, the moment SED vectors, b_{1}, ⋯, b_{m}, as defined by Eq. (B.6) in the other columns, and e^{T} = (1 0 ⋯ 0). The unconstrained foregrounds (i.e., highorder moments), which are not deprojected by cMILC, are simply mitigated through blind variance minimization, as in the NILC algorithm.
For both NILC and cMILC, we transform the Stokes Q and U fullsky simulations defined in Sect. 3 into fullsky Bmode maps before computing the fits described above, resulting in pure Bmode constraints. We then apply these methods to the two most extreme instrument configurations defined in Table 1, namely M0 and M5. Foreground cleaning is performed on the full sky for both NILC and cMILC, while the power spectra and r estimates are computed after masking f_{sky} = 50% of the NILC and cMILCBmode maps, where the mask is obtained by thresholding the observed 402 GHz Bmode map until 50% of the sky is masked.
The r estimates are derived from the angular power spectrum of the NILC/cMILCBmode map using the Gaussian likelihood
where is the Bmode power spectrum of the reconstructed NILC/cMILC map, is the postcomponentseparation noise, and C_{ℓ}(r,A_{L}) is the theoretical Bmode power spectrum as a function of r that we fit to the data, with lensing amplitude fixed to A_{L} = 1 or 0, depending on either no delensing or full delensing assumptions. The covariance matrix is given by
which implicitly includes the sample variance of the residual foregrounds and noise that are left in the reconstructed NILC/cMILC map.
All Tables
Summary of LiteBIRD HFT instrument configurations considered in this paper in terms of frequency, beam size, and sensitivity.
All Figures
Fig. 1. Sensitivity versus frequency for the various LiteBIRD configurations listed in Table 1. The black crosses indicate LFT or MFT frequencies that are not modified in this paper, while colored symbols indicate the various modified HFT frequencies. 

In the text 
Fig. 2. SEDs of the two dust models considered in this work, the 1 MBB, (red line) and the physical dust model SiFeC, (green line). The black line shows the SED of the CMB, while the gray refers to synchrotron emission. All the SEDs are in brightness temperature units and are computed considering the rms of Stokes Q and U maps on 70% of the sky. Sky masks are obtained by retaining the cleanest fraction of the sky in polarized intensity at 100 GHz. Vertical gray lines show the central frequency of the baseline LiteBIRD frequency channels, and in blue dashed lines we report the central frequency of the HFT M5 extension. The bottom panel shows the percentage difference of SiFeC with respect to MBB. 

In the text 
Fig. 3. Comparison of the reconstructed CMB Bmode maps from NILC for the M0 (top row) and M5 (bottom row) instrument configurations. The two panels in each row show opposite hemispheres aligned with the Galactic plane along the equator, and the north and south Galactic poles at the top and bottom, respectively. The left and right panels are centered on the Galactic center and anticenter. 

In the text 
Fig. 4. Commander residual (output minus input) maps for selected parameters and instrument configurations for one arbitrarily selected Monte Carlo sample. Columns show, from left to right: (1) CMB Stokes Q; (2) CMB Stokes U; (3) thermal dust spectral index; and (4) thermal dust temperature. Rows show M0, M2, and M5. The maps show one typical realization. Gray pixels indicate the masked pixels with f_{sky} = 90%. 

In the text 
Fig. 5. Power spectrum of residual foreground contamination (left panel) and noise (right panel) with NILC (solid lines) and cMILC (dashed lines) on f_{sky} = 50% of the sky for the M0 (blue) and M5 (orange) configurations. For reference, we show the input CMB Bmode signal (r = 0, black dashdotted line), the primordial Bmode signal expected from theory for values ranging from r = 10^{−2} down to r = 10^{−3} (grayshaded area), and the residual lensing cosmic variance (yellowshaded area) for no delensing (A_{L} = 1) up to 90% delensing (A_{L} = 0.1). All the spectra are binned with a multipole bin size of Δℓ = 16. 

In the text 
Fig. 6. Comparison of recovered tensortoscalar ratio distributions from NILC (solid lines) and cMILC (dotted lines) for the M0 (top panel) and M5 (bottom panel) instrument configurations. Results without delensing are shown as black curves, while results with full delensing are shown as red lines. 

In the text 
Fig. 7. Tensortoscalar ratio posterior distributions for the single MBB component model (for both simulating and modeling) as derived with Commander using multipoles ℓ = [2, 12] and a sky fraction of f_{sky} = 73% for 20 independent simulations (black lines), each drawn from a ΛCDM spectrum with r = 0. The red vertical lines correspond to the mean upper 68% confidence limit, σ_{r}. 

In the text 
Fig. 8. Tensortoscalar ratio uncertainty, σ_{r}, as a function of instrument configuration for the component separation methods considered in this paper. 

In the text 
Fig. 9. Ratio between the tensortoscalar ratio uncertainties for M5 and M0 for different analysis choices and componentseparation methods. The black points are Commander analyses using different masks and sky fractions. 

In the text 
Fig. 10. Commander residual maps (r_{ν} = d_{ν} − s_{ν}) for four different frequency channels and three different instrument configurations. The top section show residual maps for the ideal onecomponent MBB model, while the bottom section show results for the nonideal SiFeC model, where a model mismatch was included. 

In the text 
Fig. 11. Tensortoscalar ratio posterior distributions with a fiducial tensortoscalar ratio of r = 0, as derived with Commander when fitting SiFeC simulations using the onecomponent MBB model for six different instrument configurations. A χ^{2} statistic is reported for each configuration, indicating the maplevel goodness of fit. All results are derived from a sky fraction of f_{sky} = 60%. The red dashed lines correspond to the median distribution from the ideal onecomponent MBB fit to an MBB simulation shown in Fig. 7 and serves as a reference. 

In the text 
Fig. 12. Fisher uncertainty (bars) and bias (black crosses) on r for different instrumental (colors) and FGBuster componentseparation configuration (bar groups). The values on the horizontal axis define the N_{side} resolution at which the temperature of the thermal dust component is fitted for in the best 20%, 20–40% and 40–60% sky fractions. 

In the text 
Fig. A.1. Example of the results from the FGBuster componentseparation runs. Two instrumental configurations are considered: the baseline M0 (blue) and the high frequency extension M3 (red). The solid lines represent the postcomponentseparation noise and foreground residuals S_{ℓ}, computed by averaging 1000 simulations for the baseline, and averaging 10 simulations for the extension. Together with the power from gravitational lensing and primordial Bmodes, they define the error bars on the CMB power spectrum (boxes). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.