<figure><img src="image_38.png" /><figcaption>**Fig. 36.** Foreground emission surrounding the Crab Nebula. Left: 353-GHz *Planck* map (Stokes I) centered on the Crab Nebula, displaying foreground emission in its surroundings, compared to the typical beam size of *LiteBIRD*. Upper right: Polarization angle computed by convolving the *Planck* Stokes parameters maps at 337 GHz with the *LiteBIRD*-HFT beam, as a function of the integration radius. Lower right: the same at 100 GHz (MFT).</figcaption></figure>

surements at different frequencies, and assuming no variation of the polarization angle with frequency, they derived a value for the Crab's polarization angle of $-87.7^\circ \pm 0.3^\circ$. This error of 18 arcmin is larger than the *LiteBIRD* requirement. Ref. [264] discussed the prospects for improving the global uncertainty to a level of $\simeq 6$ arcmin through the addition of new measurements with individual error bars of around $0.2^\circ$. This could be attempted through independent ground-based observations, using facilities like NIKA2, SCUBA-2, or SRT [264–266].

Such observations should also allow us to study other possible errors derived from incorrect modeling hypotheses, in particular: (i) the independence of the Crab's polarization angle with frequency; (ii) the time dependence (it is well known that Crab's total-intensity flux density fades with time, but little is known about polarization); and (iii) the impact of the background emission that will affect the *LiteBIRD* measurement due to its coarse beam. Therefore, ideally these observations must: (i) cover a large frequency range; (ii) be spread in time; and (iii) cover a sufficiently large region around the Crab Nebula ($\simeq 2^\circ$ across).

The effect of the background contamination is a relevant one, in particular at high frequencies, and we have studied this in detail [267]. This is illustrated in Fig. 36 (left panel), where dust emission in the neighborhood of the Crab Nebula is evident using 353-GHz data from *Planck*. This emission is partially polarized, and we have analyzed its effect using the Stokes parameter maps measured by *Planck* at the relevant frequencies [268,269], convolved with the expected *LiteBIRD* beams. In order to assess the contribution of surrounding diffuse emission we have estimated the measured polarization angle, after convolving the maps with the beam over increasing radii. The resulting measured angle is significantly different from the polar-
ization angle of the nebula itself, and is stabilized for integration radii of the order of 1° (see Fig. 36, center and right panels). The depth of the required surveys of the region has been investigated adding white noise to each 1.7 arcmin pixel in the Stokes Q and U maps and computing the measured angle obtained by convolution of the Q and U maps with the LiteBIRD beams. Repeating the process 10 000 times, we find that the standard deviation of the distribution of measured angles is smaller than the required 3 arcmin if the noise per pixel in the reference Stokes parameter maps is below about 30 µK rms at 100 GHz, and 10 µK rms at 337 GHz. These stringent requirements motivate us to explore further detailed measurements of the Crab using forthcoming surveys in order to have the reference maps ready at the time of the LiteBIRD mission.

#### 4.8.4.4 Prospects for ground calibration

As shown in Fig. 30, the general plan is to characterize the relative polarization angles per focal plane during the end-to-end (E2E) cryogenic operations at the telescope level. Special care will be taken to define a common reference frame or reference transposition from subsystem characteristics to the telescopes, up to the PLM-SVM-integrated instrument.

For reference, the end-to-end ground-calibration test campaign for Planck-HFI ultimately achieved an average of 1.5° accuracy while the subsystem level (i.e., single cooled pixel) testing had achieved an average accuracy of 0.1° [270,271]. For LiteBIRD, since we are using both an HWP and sinuous antenna [272], the bandpasses and the polarization angle are therefore highly coupled. It is possible to test for the overall spectral efficiency of a channel with an unpolarized source, and it is possible to accurately measure the maximum efficiency angle of a given detector, but for an end-to-end test of the whole instrument the two effects are combined in a non-trivial way. We therefore plan to cleanly disentangle them via polarized spectroscopy to obtain the two spectro-polarimetric datacubes [273]. The plan is to test the system with two orthogonal inputs (Fourier-transform spectrometer (FTS), VNA, or other coherent source) and to measure spectra at each small variation of HWP angle position in order to be able to reconstruct the expected modulation curve for any given spectral index source.

The end-to-end calibration operations are complementary to the RF test campaign (cf. Sect. 4.8.2), which will specifically address the impact of the optical distortion and coating inefficiencies. A measurement campaign has already been performed on the LFT small-scale model and is reported in Ref. [174].

In addition, various instrumental parameters may also impact the polarization angle and will be characterized at the component level to check their accuracy with respect to the requirements flow down, in particular on the knowledge of the encoding angle HWP PMU (cf. Sects. 4.2.4 and 4.3.3). In parallel, we are still working on the design and especially on the HWP, since any HWP frequency-dependent fast axis within the band can result in a mixing of the Q and U polarization states. Nevertheless, for a sapphire Pancharatnam HWP, there is a recipe that nulls the frequency response, and the LFT PMU employs such a recipe to minimize this effect [196], considering the fact that the residual variation may still produce E-to-B leakage [274]; the study explicitly for LiteBIRD is in progress.

#### 4.8.5. *Calibration development plan.*
The previous sections addressed the most challenging operations, but other instrumental parameters will need to be measured in the end-to-end cold integrated campaign. These are summarized in Table 11, together with the experimental setup
<table><caption>**Table 11.** Instrumental parameters and calibration strategy to assess them in order to meet the accuracy required to meet the systematics budget beyond subsystem testing.</caption><thead><tr><th></th><th>Instrumental param.</th><th>Instrument-level plans in cryogenic facility</th><th>In-flight plan</th></tr></thead><tbody><tr><td rowspan="3">Beam (co- and cross-polar)</td><td>Far sidelobes</td><td>PLM cross-check</td><td>Absolute gains</td></tr><tr><td>Near sidelobes</td><td>RF</td><td>Planet</td></tr><tr><td>Main beam</td><td>RF</td><td>Planet</td></tr><tr><td rowspan="4">Polarization</td><td>Instrumental polarn.</td><td>E2E/DataCube</td><td>C<sub>l</sub><sup>EB</sup> and Crab</td></tr><tr><td>Absolute angle</td><td>E2E/DataCube</td><td>C<sub>l</sub><sup>EB</sup> and Crab</td></tr><tr><td>Relative angle</td><td>E2E/DataCube</td><td>C<sub>l</sub><sup>EB</sup> and Crab</td></tr><tr><td>Efficiency</td><td>E2E/DataCube</td><td>C<sub>l</sub><sup>EB</sup> and Crab</td></tr><tr><td rowspan="2">Gain</td><td>Relative gain in time</td><td>E2E/ΔT of the facility</td><td>Correlation with HK</td></tr><tr><td>Relative and absolute gains</td><td>E2E/polar. source</td><td>C<sub>l</sub> norm. + Galaxy + planet + dipole</td></tr><tr><td>Cosmic ray</td><td>Cosmic-ray glitches</td><td></td><td>Glitches</td></tr><tr><td>Spectral response</td><td></td><td>E2E/FTS/DataCube</td><td></td></tr><tr><td rowspan="5">Transfer function</td><td>Noise PSD</td><td>E2E</td><td>All data</td></tr><tr><td>Detector time const.</td><td>E2E/chopped source</td><td>Glitches, planets</td></tr><tr><td>Crosstalk</td><td>E2E/polar. source</td><td>Glitches, planets</td></tr><tr><td>Linearity</td><td>E2E/polar. source</td><td>Glitches, planets</td></tr><tr><td>Sensitivity to bkg</td><td>E2E/ΔT of the load</td><td>Galactic crossings</td></tr></tbody></table>

planned to be used to characterize them and the proposed in-flight reconstruction. A dedicated characterization phase between the launch and arrival at L2 is also planned as a performance-verification phase prior to nominal observations; this phase is not detailed in this paper.

Finally, a viability study (not in the baseline plans of the calibration strategy) has been performed for a calibration satellite (L2-CalSat), flying in formation at L2 and able to emit a reference signal [275]. The calibration satellite is based on the CubeSat standard and has been conceived to travel as a piggy-back integrated on the service platform of the main satellite, *LiteBIRD* itself. A 6U CubeSat is expected to meet the volume needs of the ancillary satellite. An alternative implementation option would be a micro-satellite capable of reaching L2 autonomously and allowing a major flexibility at the system-design level. However, a major drawback would be the much larger amount of required propellant that would imply a significantly higher cost. The calibration satellite total mass would be 7.2 kg and the required electrical power would be provided by two steerable solar panels and accumulated in a battery module. The solar panels would continuously generate 30 W, enough to feed all the subsystems of L2-CalSat and charge the battery to provide the peak power required during calibration. L2-CalSat is required to be at a sufficiently far distance, of about 6 km, in order to have thermal power levels at least four times lower than the saturation power of the detectors (thermal power levels are maximum at the highest-frequency bands). The calibration satellite could help to provide strong control over many systematic errors, including those related to gain, non-linearity, spectral response, pointing, beam patterns, or absolute polarization angle. In particular, our analysis shows that the polarization angle error can be reduced below 1 arcmin for each detector. In relation to beam-pattern characterization, the calibration sources could provide noise-floor levels lower
<table><caption>**Table 12.** Baseline observation strategy and parameters related to the sampling rate.</caption><thead><tr><th rowspan="2">α [deg.]</th><th rowspan="2">β [deg.]</th><th rowspan="2">Precession rate [min.]</th><th rowspan="2">Spin rate [rpm]</th><th colspan="3">HWP revolution rate [rpm]</th><th rowspan="2">Sampling rate [Hz]</th></tr><tr><th>LFT</th><th>MFT</th><th>HFT</th></tr></thead><tbody><tr><td>45</td><td>50</td><td>192.348</td><td>0.05</td><td>46</td><td>39</td><td>61</td><td>19</td></tr></tbody></table>

than –70 dB in the overall band of interest (from 40 to 400 GHz) with only one calibration session of about one day’s duration. This would assure a complete and accurate characterization of the beam and in particular of the far sidelobes, placed at more than 10° from the beam center. We will assess the feasibility as the conceptual study progresses.

## 5. Statistical and systematic uncertainties in the tensor-to-scalar ratio measurement

This section describes the way that the forecasts on *r* are derived. Section 5.1 gives an overview of the model of the instrument that is used for the simulations, with a particular focus on the expected statistical noise derivation. Section 5.2 summarizes the assumptions made for the foreground modeling, as well as the strategy concerning component separation. In Sect. 5.3, a review of the main systematic studies and their impact on the derivation of the systematic error on *r* is detailed. Finally, in Sects. 5.4 and 5.5, a summary of the total expected error on *r* is given, without and with the use of external data sets.

### 5.1. Instrument model

The LiteBIRD instrument model (IMo) is a quantitative description of the entire LiteBIRD experiment: the three telescopes (LFT, MFT, and HFT); the payload; and the observational strategy. The IMo is used to store and track changes of the instrument design and is used in and interfaced to the forecasts/simulations/data analysis tools. Its use is required to ensure consistency between all the instrumental parameters, and it gives the best current description of the instrumental design. This section gives an overview of the main parameters entering the IMo.

Table 12 gives the values of the main baseline parameters common to the three telescopes, LFT, MFT, and HFT. The parameters related to the observational scan strategy are further defined by the schematics shown in Fig. 37. The duration of the mission is assumed to be three years throughout this section, except where explicitly stated otherwise.

We optimize the scan strategy parameters including α, β, and the spin rate so as to have the distributions as uniform as possible for the hits over the entire sky area and for the scanning directions in individual sky pixels. A spin rate less than 0.05 rpm is found to give non-uniform distributions for both quantities. The precession rate is determined to avoid moiré patterns using the procedure described in Ref. [276]. The HWP spin rates are set so as not to overlap the HWP harmonics with the science observation frequency band determined by the beam FWHM. The sampling rate is twice the Nyquist frequency of the HWP modulation.

Table 13 summarizes the expected noise and sensitivities of the different frequency bands of LiteBIRD, together with related key parameters such as beam size (calculated using physical optics simulations similar to those described in Ref. [277]) and number of detectors. The details of the sensitivity calculations are given in Sect. 5.1. By default the bandpasses are assumed to be top-hat functions whose bandwidths are given in the same table.
<figure><img src="image_39.png" /><figcaption>Fig. 37. Schematic of the observational parameters. The telescope boresight is at an angle $\beta = 50^{\circ}$ to the spin axis and rotates at a rate of 0.05 rpm. The spin axis is rotated around the anti-Sun direction through precession with an angle $\alpha = 45^{\circ}$ in 3.2058 h. The anti-Sun axis rotates around the Sun in one year. With a combination of the three motions, the boresight can cover the entire sky in half a year.</figcaption></figure>

The noise-equivalent power (NEP) of each detector is predicted from a quadrature sum of the expected photon noise (NEPph), thermal carrier noise of the bolometer (NEPth), readout noise (NEPread), and extra noise (NEPext) sourced from the environment and unknowns. The NEP of a single detector (NEPdet) may therefore be expressed as

$$NEP_{det} = \sqrt{NEP_{ph}^2 + NEP_{th}^2 + NEP_{read}^2 + NEP_{ext}^2}, \quad (27)$$

where all the components are assumed to be mutually uncorrelated. The quantity NEPph is deduced from the expected loading on the detectors (taking account of the various optical elements, their intrinsic characteristics, and the temperature to which they are cooled), while NEPth is derived from the temperature of the bath and the expected critical temperature of the bolometers. The readout noise and the external noise are required to increase the internal detector noise (encompassing the detector and the thermal carrier noise) by no more than 10% and the total detector noise by 15%, respectively.

We define the noise-equivalent temperature (NET) of a detector as

$$NET_{det} = \frac{NEP_{det}}{\sqrt{2}(dP/dT_{CMB})}, \quad (28)$$

where we define the conversion factor from power to CMB thermodynamic temperature through

$$\frac{dP}{dT_{CMB}} = \int_{v_1}^{v_2} \left[ k_B \eta(v) \left( \frac{h\nu}{k_B T_{CMB} (e^{h\nu/k_B T_{CMB}} - 1)} \right)^2 e^{h\nu/k_B T_{CMB}} \right] dv. \quad (29)$$

The parameters $v_1$ and $v_2$ represent the lower and upper edges of the frequency band, respectively, while $\eta$ is an overall optical efficiency factor. Detector NETs for each band are shown in Table 3. The NET for an array is calculated as

$$NET_{arr} = \frac{NET_{det}}{\sqrt{N_{det} \times 0.8}}, \quad (30)$$

where $N_{det}$ is the number of detectors in the given frequency band and the factor of 0.8 represents a degradation factor for the detector yield. The NETs of a frequency band covered by
<table><caption>**Table 13.** LiteBIRD sensitivities. We show values related to the sensitivity in 15 observational bands. From left to right the columns are: the telescope covering the band (or “comb.” for combined values); the band identification number; the band center frequency in GHz; the bandwidth in GHz and its ratio to the center frequency; the main beam size FWHM in arcmin; the total number of detectors; the array NET in $\mu K\sqrt{s}$; and the sensitivity in $\mu K$-arcmin. The two values of the beam size and the sensitivity for the LFT 68-, 78-, and 89-GHz bands in parentheses are for the two detector pixel sizes of 16 and 32 mm. See Table 3 for the detector pixel sizes and detector NETs.</caption><thead><tr><th></th><th>ID</th><th>ν<br>[GHz]</th><th>δν [GHz]<br>(δν/ν)</th><th>Beam size<br>[arcmin]</th><th>No. of<br>detectors</th><th>NET<sub>arr</sub><br>[μK√s]</th><th>Sensitivity<br>[μK-arcmin]</th></tr></thead><tbody><tr><td>LFT</td><td>1</td><td>40</td><td>12 (0.30)</td><td>70.5</td><td>48</td><td>18.50</td><td>37.42</td></tr><tr><td>LFT</td><td>2</td><td>50</td><td>15 (0.30)</td><td>58.5</td><td>24</td><td>16.54</td><td>33.46</td></tr><tr><td>LFT</td><td>3</td><td>60</td><td>14 (0.23)</td><td>51.1</td><td>48</td><td>10.54</td><td>21.31</td></tr><tr><td>LFT<br>comb.</td><td>4</td><td>68</td><td>16 (0.23)</td><td>(41.6, 47.1)</td><td>(144, 24)</td><td>(9.84, 15.70)<br>8.34</td><td>(19.91, 31.77)<br>16.87</td></tr><tr><td>LFT<br>comb.</td><td>5</td><td>78</td><td>18 (0.23)</td><td>(36.9, 43.8)</td><td>(144, 48)</td><td>(7.69, 9.46)<br>5.97</td><td>(15.55, 19.13)<br>12.07</td></tr><tr><td>LFT<br>comb.</td><td>6</td><td>89</td><td>20 (0.23)</td><td>(33.0, 41.5)</td><td>(144, 24)</td><td>(6.07, 14.22)<br>5.58</td><td>(12.28, 28.77)<br>11.30</td></tr><tr><td>LFT/<br>MFT<br>comb.</td><td>7</td><td>100</td><td>23 (0.23)</td><td>30.2/<br>37.8</td><td>144/<br>366</td><td>5.11/<br>4.19<br>3.24</td><td>10.34<br>8.48<br>6.56</td></tr><tr><td>LFT/<br>MFT<br>comb.</td><td>8</td><td>119</td><td>36 (0.30)</td><td>26.3/<br>33.6</td><td>144/<br>488</td><td>3.8/<br>2.82<br>2.26</td><td>7.69<br>5.70<br>4.58</td></tr><tr><td>LFT/<br>MFT<br>comb.</td><td>9</td><td>140</td><td>42 (0.30)</td><td>23.7/<br>30.8</td><td>144/<br>366</td><td>3.58/<br>3.16<br>2.37</td><td>7.25<br>6.38<br>4.79</td></tr><tr><td>MFT</td><td>10</td><td>166</td><td>50 (0.30)</td><td>28.9</td><td>488</td><td>2.75</td><td>5.57</td></tr><tr><td>MFT/<br>HFT<br>comb.</td><td>11</td><td>195</td><td>59 (0.30)</td><td>28.0/<br>28.6</td><td>366/<br>254</td><td>3.48/<br>5.19<br>2.89</td><td>7.05<br>10.50<br>5.85</td></tr><tr><td>HFT</td><td>12</td><td>235</td><td>71 (0.30)</td><td>24.7</td><td>254</td><td>5.34</td><td>10.79</td></tr><tr><td>HFT</td><td>13</td><td>280</td><td>84 (0.30)</td><td>22.5</td><td>254</td><td>6.82</td><td>13.80</td></tr><tr><td>HFT</td><td>14</td><td>337</td><td>101 (0.30)</td><td>20.9</td><td>254</td><td>10.85</td><td>21.95</td></tr><tr><td>HFT</td><td>15</td><td>402</td><td>92 (0.23)</td><td>17.9</td><td>338</td><td>23.45</td><td>47.45</td></tr><tr><td>Total</td><td></td><td></td><td></td><td></td><td>4508</td><td></td><td>2.16</td></tr></tbody></table>

multiple telescopes are combined into the average NET$_{comb}$ as

$$ \text{NET}_{\text{comb}} = \sqrt{\frac{1}{\sum_i (1/\text{NET}_i^2)}}, \qquad (31) $$

where the index *i* runs over the telescopes of LFT/MFT/HFT.

Finally, the statistical sensitivity ($\sigma_S$) in $\mu K$-arcmin is given as

$$ \sigma_S = \sqrt{\frac{4\pi f_{\text{sky}} 2 \text{NET}_{\text{comb}}^2}{t_{\text{obs}}} \left(\frac{10\,800}{\pi}\right)}, \qquad (32) $$

where we assume the sky coverage fraction $f_{\text{sky}} = 1.0$ and $t_{\text{obs}} = 3$ yr (94 672 800 s) $\times \eta$. Here $\eta$ is the observation efficiency, derived by taking into account degradation factors on the observation duty cycle ($\eta_{\text{duty}} = 0.85$), the data loss due to cosmic rays ($\eta_{\text{cr}} = 0.95$), and a margin $\eta_{\text{margin}} = 0.95$. Hence we have $\eta = \eta_{\text{duty}}\eta_{\text{cr}}\eta_{\text{margin}} = 0.77$. We note that the factor of $\sqrt{2}$ here originates from the fact that the polarization is measured using two independent pieces of data.
The total sensitivity is obtained by using the NET values combined with all the NET$_{comb}$ values of individual frequency channels using Eq. (31).

## 5.2. Foreground cleaning

In this section, we describe the reference sky that we adopt, and which elements from the IMo in Sect. 5.1 we exploit in order to generate the LiteBIRD-simulated multi-frequency input maps for treatment through component separation. We then describe various methodologies, as well as the products of these analyses that are exploited in the following sections.

### 5.2.1. Reference skies.
We adopt a native resolution of all maps corresponding to HEALPix $N_{side} = 512$ [278]. The diffuse emission, corresponding to CMB and Galactic foregrounds, is in thermodynamical CMB temperature units, $\mu K_{CMB}$. The CMB angular power spectra correspond to the Planck 2018 best-fitting parameters [279], with the six $\Lambda$CDM parameters corresponding to {$H_0$, $\Omega_b h^2$, $\Omega_c h^2$, $\tau$, $n_s$, $A_s$} = {67.36 km s$^{-1}$ Mpc$^{-1}$, 0.022 37, 0.12, 0.0544, 0.965, 2.099 $\times$ 10$^{-9}$}, with lensing and with no tensors, i.e., $r = 0$.

Diffuse foregrounds are produced with the public package Python Sky Model (PySM) with default implementation of the amplitude of foregrounds, and spatially varying spectral indices [280,281], corresponding to the d1sl parametrization in PySM. The Galactic polarized emission is composed of thermal dust and synchrotron models, dominating the polarization emission at frequencies higher and lower than about 70 GHz, respectively [282]. In this work, the polarized synchrotron emission is modeled as

$$[Q_s, U_s](\hat{n}, \nu) = [Q_s, U_s](\hat{n}, \nu_*) \left(\frac{\nu}{\nu_*}\right)^{\beta_s(\hat{n})}, \quad (33)$$

where $[Q_s, U_s](\hat{n}, \nu)$ are the Stokes parameters of synchrotron emission in Rayleigh-Jeans (RJ) temperature units on a given line of sight $\hat{n}$ at a frequency $\nu$, $\nu_*$ is a reference frequency, and $\beta_s(\hat{n})$ is the spatially varying synchrotron spectral index. On the 50% sky mask shown in Fig. 38, the PySM template gives $\beta_s = -2.993 \pm 0.046$. The polarized thermal dust emission is modeled as a modified blackbody, corresponding to the expressions

$$[Q_d, U_d](\hat{n}, \nu) = [Q_d, U_d](\hat{n}, \nu_*) \left(\frac{\nu}{\nu_*}\right)^{\beta_d(\hat{n})-2} \frac{B[\nu, T_d(\hat{n})]}{B[\nu_*, T_d(\hat{n})]}, \quad (34)$$

where $[Q_d, U_d](\hat{n}, \nu)$ are the Stokes parameters of thermal dust emission in RJ temperature units, $\beta_d(\hat{n})$ is the dust spectral index (power-law index of dust emissivity), $T_d(\hat{n})$ is the effective temperature of dust, and $B(\nu, T_d)$ is a blackbody spectrum with a temperature $T_d$. Similarly to the synchrotron spectral index, PySM templates show $\beta_d = 1.534 \pm 0.028$ and $T_d = 22.23 \pm 1.56$ K variations across the $f_{sky} = 49.5\%$ mask shown in Fig. 38. The total intensity of dust and synchrotron scale in frequency just like polarization. According to predictions [283,284] obtained using the publicly available Point Source ForeCast package,<fnref n="13" /> polarized compact sources are expected not to play a major role in LiteBIRD measurements; bright radio sources will need to be removed with dedicated filters, and for the infrared populations only upper limits exist in polarization (see e.g., Ref. [262]). Thus, we do not include their treatment in the present work, focusing on the capabilities of LiteBIRD to deal with diffuse foregrounds.

<footnote n="13">¹³giuspugl.github.io/ps4c/index.html.</footnote>
<figure><img src="image_40.png" /><figcaption>**Fig. 38.** Top panel: Masks used for foreground analyses. The gray region shows the sky area excluded from the foreground analysis, while the colors show areas with homogeneous resolution for the foreground spectral parameters. Middle (bottom) panel: map of the mean (standard deviation) of the foreground residuals.</figcaption></figure>

5.2.2. *Input maps.* The templates of the reference sky described above are scaled in frequency using the *PYSM* prescriptions and combined with instrumental features from the *IMo*. These are specifically:
* bandpass integration;

* beam convolution;

* addition of white noise.

The IMo objects used for the above operations are described in Sect. 5.1 and we report here their exploitation for achieving the series of multi-frequency maps to be processed by component separation. All templates are convolved with top-hat bandpasses (see Table 13 for the band centers and widths) and co-added to form input sky maps. Concerning beam convolution, two series of maps are produced, for comparison purposes. The first one is a direct convolution with a Gaussian circular beam shape with FWHMs from the IMo. The second one is a timeline-based series of simulations, generated with the public package Time-Ordered Astrophysics Scalable Tools (TOAST).<fnref n="14" /> The convolution of the signal with the beam simulated in GRASP (encoding for the whole $4\pi$ beam pattern) is performed with the libconvipt<fnref n="15" /> library, implementing beam convolution on a sphere based on Ref. [285].

For the first set of simulations, the noise was added to the co-added multi-frequency and multi-component maps for the Gaussian circular beam convolutions. A total of 1000 noise and CMB simulations were made available using the procedures described in the next subsection. We use the second set for a study of the beam far-sidelobe systematic effects in Sect. 5.3.3.

### 5.2.3. Component separation.
CMB extraction and diffuse foreground reconstruction may be achieved through broad classes of techniques (see Ref. [40] and references therein). Parametric fitting consists of constructing a sky model on the basis of a suitable parametrization of foreground and CMB unknowns and estimating the latter by exploiting a multi-frequency data set. Internal linear combination and template-fitting methods consist of the minimization of a linear mixture of the multi-frequency data and the subtraction of foreground templates that come from external data or are derived internally, respectively. Baseline results for the present paper were obtained with such a parametric fitting method, as discussed below.

In the parametric map-based component-separation tool “ForeGround Buster” (FGBuster<fnref n="16" /> [286]), our approach is to assume a foreground model with three spectral parameters determining their spectral energy distributions (SEDs): a spectral index $\beta_d$ and a temperature $T_d$ of a modified blackbody for dust; and the spectral index $\beta_s$ of a power law for synchrotron. These spectral parameters are assumed to be constant over Healpix pixels at a given $N_{side}$. Of course, this represents an approximation of the model, whereas the foreground spectral parameters will in general vary smoothly across the sky. In order to balance the statistical uncertainties sourced by the number of free parameters and the need for a flexible enough model capturing the input sky complexity (especially the spatial variability of the SEDs), we fit the different spectral parameters at different resolutions depending on Galactic latitude. Since the optimal balance depends on the sky region, we divide the sky into (almost) isolatitude areas covering 20% of the sky each; we rank them by their expected level of foreground contamination and consider only the three least contaminated. The $N_{side}$ parameters for these regions are given in Table 14. Here the three columns correspond to the three 20% parts of the sky from lowest to highest Galactic latitude, defined from the Planck post-processing masks adopted for component separation

<footnote n="14">¹⁴hpc4cmb.github.io/toast.</footnote>

<footnote n="15">¹⁵https://github.com/hpc4cmb/libconvipt.</footnote>

<footnote n="16">¹⁶github.com/fgbuster/fgbuster.</footnote>
<table><caption>**Table 14.** Values of $N_{\text{side}}$ for foreground parameter variations in three different 20% large regions of the sky. The $N_{\text{side}}$ values are chosen to capture most of the spatial variability of the spectral indices and to keep the resulting foreground bias under control while not increasing the statistical uncertainty. The angular resolution for $\beta_d$ is kept to the maximum allowed ($N_{\text{side}} = 64$), since *LiteBIRD* turns out to be particularly sensitive to this parameter.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Galactic latitude</th></tr><tr><th>Low</th><th>Medium</th><th>High</th></tr></thead><tbody><tr><th>β<sub>d</sub></th><td>64</td><td>64</td><td>64</td></tr><tr><th>T<sub>d</sub></th><td>8</td><td>4</td><td>0</td></tr><tr><th>β<sub>s</sub></th><td>4</td><td>2</td><td>2</td></tr></tbody></table>

[40], referred to as GAL20, GAL40, and GAL60, available on the Planck Legacy Archive,<fnref n="17" /> and shown in Fig. 38. The frequency coverage and sensitivity of *LiteBIRD* allow us to estimate foreground parameters and CMB amplitudes in suitably small sky areas or resolution elements. Thus, the number of non-linear parameters, {$\beta_d$, $T_d$, $\beta_s$}, as driven by the resolutions in Table 14, computed in the $f_{\text{sky}} = 49\%$ observed sky, is precisely 24 545. The number of $Q/U$ amplitudes contained in the recovered CMB, dust, and synchrotron maps is 144 534. Therefore, in total, the component-separation process deals with 169 079 free parameters. This can be compared to the $10^6$-long data vector as formed by the observed $Q/U$ amplitudes in the 22 frequency channels.

The unmodeled foreground component leaks into the CMB reconstruction and sources the foreground systematic residuals reported in Fig. 39 and can eventually give a non-zero bias on $r$, as introduced in Ref. [287]. The noisy estimation of the spectral parameters, instead, sources the so-called “statistical foreground residuals” and eventually the uncertainty $\sigma(r)$. We estimate the latter with a Fisher-matrix approach, where the total covariance $C$ contains CMB lensing, noise, and statistical residuals. Namely, we calculate the second derivative of the likelihood as $F = -\langle \partial^2 \log(L) / \partial r^2 \rangle$ and obtain

$$ \sigma(r=0) = 1/\sqrt{F}. \qquad (35) $$

Here, $-2 \log L = f_{\text{sky}} (\text{tr}C^{-1}\mathbf{D} + \log \det C)$, where $\mathbf{D}$ is the simulated data vector [288], calculated using the noise and residuals spectra averaged over 1000 foregrounds+noise realizations. For each of them, any non-zero recovered CMB $B$-modes are interpreted as foreground leakage.

5.2.4. Products. Once the residuals are estimated on 60% of the sky, their power spectra are evaluated in a sky region corresponding to $f_{\text{sky}} = 49.5\%$, obtained by intersecting the 60% mask with foreground residuals estimated through component separation (see Fig. 38). The mask is derived by: (1) removing regions with the highest foreground residuals; (2) smoothing it with a 2° beam; and (3) converting it to a binary mask using a threshold equal to 0.5. The threshold for step (1) is chosen so that the bias on $r$ is low enough while keeping the sky coverage large enough to reduce as much as possible the associated $\sigma(r)$. Although this exact methodology will not be applicable to real data sets, a better understanding of foreground morphology, as well as additional observational/theoretical priors on their SEDs, will help us to better define the Galactic masks by the time of the *LiteBIRD* launch. The angular power spectra of

<footnote n="17">¹⁷pla.esac.esa.int.</footnote>
<figure><img src="image_41.png" /><figcaption>**Fig. 39.** LiteBIRD error bars on the primordial *BB* power spectrum for *r* = 0.00461 and *τ* = 0.0544, where *D*<sub>ℓ</sub> = ℓ(ℓ + 1)C<sub>ℓ</sub>/2π. There are two types of error bars: cosmic variance (including primordial and lensing *B*-mode variance), shown as gray boxes, and total error bars (including cosmic variance and averaged, noisy foreground residuals), shown in red. The input foreground maps are described in Sect. 5.2.1. The solid orange curve represents the average total residuals (including statistical and systematic foreground residuals, as well as noise after component separation, accounted for by *N*<sub>ℓ</sub> in Eq. (38)), with the light band corresponding to the scatter across noise simulations. The dashed orange curve represents the noiseless foreground residuals, interpreted as the bias on *B*-mode power, eventually leading to a bias on the tensor-to-scalar ratio. Although the contribution of the statistical foreground residuals, noise, and lensing *B*-modes are included in the error bars, the latter are centered on the theoretical primordial *BB* curve.</figcaption></figure>

the foreground residuals are estimated using pseudo-$C_\ell$ corrected for the mask, evaluated using the Healpy Anafast<sup>18</sup> public code. We checked that the power spectra of foreground residuals do not change if they are estimated by methods such as NaMaster,<sup>19</sup> mainly because foreground *E*- and *B*-modes have similar amplitudes.

Figure 39 shows the result of our foreground-cleaning method applied to the multi-frequency data set specified in Sect. 5.2.2 in the case of Gaussian beam convolution. We display the expected error bars on the primordial *BB* spectra for *r* = 0.00461, corresponding to the reference Starobinsky model introduced in Sects. 2.4, 2.5 [14], and *τ* = 0.0544 [289]. The error bars are shown for pure cosmic variance, including the primordial and lensing *B*-mode variance, as well as for the total uncertainty coming from cosmic variance plus noisy foreground residuals. The total residuals from foregrounds and noise after component separation, averaged over the 1000 simulations, are shown in Fig. 39 (orange area, with average given by the dark orange line), as well as the noiseless foreground residuals (dashed orange line), otherwise known as “system-

<footnote n="18"><sup>18</sup>https://healpy.readthedocs.io/en/latest/generated/healpy.sphtfunc.anafast.html.</footnote>

<footnote n="19"><sup>19</sup>github.com/LSSTDESC/NaMaster.</footnote>
atic foreground residuals”. The first element determines the accuracy on *r* via a power spectrum estimation, namely the statistical error including foreground cleaning and noise, while the second is responsible for the bias on the recovered *BB* spectrum and thus on *r*. Due to their large amplitudes, the statistical foreground residuals, sourced by statistical errors on the estimated spectral indices, dominate the overall error budget.

From these averaged spectra, the corresponding bias on $r_{FG}$, as well as the accuracy $\sigma(r=0)$ after foreground cleaning, is given by $r_{FG} = (3.3 \pm 6.2) \times 10^{-4}$. This result is obtained for an input value $r=0.00$, while the same analysis with an input $r=0.01$ leads to $\sigma(r=0.01)=0.0013$. We find $\sigma(r=0.01)=0.0036$ from the reionization bump only ($\ell < 20$) and $\sigma(r=0.01)=0.0014$ from the recombination bump ($\ell > 20$). Note that $\sigma(r=0.01)$ is the measurement uncertainty including cosmic variance. Since the posterior probability density of $r$ is asymmetric and falls rapidly towards $r=0$, we expect detections of both bumps at a significance exceeding $5\sigma$ (see Sect. 5.4).

Relative to the recombination bump, the reionization feature becomes more relevant as one approaches the Starobinsky limit, since the recombination bump would require substantial delensing in order to be measured. We will come back to delensing effects in Sect. 5.5, where the averaged total and noiseless foreground residuals mentioned above are the inputs to the estimation of the total uncertainty on the tensor-to-scalar ratio $r$.

Along with the recovered CMB map, our component-separation approach also provides estimates of the Galactic foregrounds, and in particular of the dust map. This information is used in Sect. 5.4 to marginalize our estimate of $r$ over the foreground bias.

### 5.3. *Instrumental systematic uncertainties*

We describe here our analysis of instrumental systematic effects. Foreground cleaning and systematic residuals will constitute the total error budget, as described in Sect. 5.4. Throughout this section, we describe specific definitions and objects for each of the systematic effects that we consider. Common to the whole analysis, simulations of the sky include Galactic foreground components as an input, following the model described in Sect. 5.2.1, combined as described in Sect. 5.2.2 to produce simulated frequency-channel data.

#### 5.3.1. *Introduction to systematic uncertainties.*
In this section, we describe potential systematic effects for the *LiteBIRD* mission and evaluate their impact on the final results. Systematic effects originate from imperfect knowledge of the foregrounds, combined with incomplete correction of instrumental or environment effects arising from either inaccurate modeling or limited knowledge of instrumental parameters and environmental contributions.

The main objective of *LiteBIRD* is to measure the tensor-to-scalar ratio $r$. We therefore focus our study on the impact of systematic effects on the measurement of this parameter. We particularly estimate the systematic contribution on the measurement bias $\Delta r$ under the condition of the true value being $r=0$. We give the detailed definition of $\Delta r$, as well as the total error $\delta r$, in the next section.

We divide the systematic effects into several categories. The systematic bias produced by the component-separation residual was described in Sect. 5.2. We describe the remaining systematic sources in this section. Table 15 shows a summary of the most relevant systematic effects as the top-level category. Beam systematic effects are mainly caused by imperfections in the calibrated
<table><caption>**Table 15.** Summary of sources of systematic effects as the top-level category. We categorize the systematic effects in the first column. The details of individual categories are shown in the second column. The third column shows the systematic biases on *r*, except for the cosmic rays, which are supposed to yield an extra noise source. The fourth column describes the sources of systematic effects. The last column shows the types. Types “E” and “R” state the expectation and the requirement, respectively, as described in the text. Note that we assign an error budget of 1% of the total systematic bias $5.7 \times 10^{-6}$ to individual systematic effects as a nominal error budget in order to account for dozens of potential systematic sources. Details are described in the text.</caption><thead><tr><th>Category</th><th>Systematic effect</th><th>Δr</th><th>Source</th><th>Type</th></tr></thead><tbody><tr><td rowspan="5">Beam</td><td>Far sidelobes</td><td>4.4 × 10<sup>-5</sup></td><td>B → B, E → B</td><td>R</td></tr><tr><td>Near sidelobes</td><td>5.7 × 10<sup>-6</sup></td><td>B → B, E → B</td><td>R</td></tr><tr><td>Main lobe</td><td>&lt;10<sup>-6</sup></td><td>E → B</td><td>E</td></tr><tr><td>Ghost</td><td>5.7 × 10<sup>-6</sup></td><td>E → B</td><td>R</td></tr><tr><td>Polarization and shape in band</td><td>&lt;10<sup>-6</sup></td><td>E → B</td><td>R</td></tr><tr><td>Cosmic ray</td><td>Cosmic-ray glitches</td><td>Noise</td><td>Power to B, E</td><td>E</td></tr><tr><td rowspan="4">HWP</td><td>Instrumental polarization</td><td>&lt;10<sup>-6</sup></td><td>T → B</td><td>E</td></tr><tr><td>Transparency in band</td><td>5.7 × 10<sup>-6</sup></td><td>E → B</td><td>R</td></tr><tr><td>Polarization efficiency in band</td><td>5.7 × 10<sup>-6</sup></td><td>B → B</td><td>R</td></tr><tr><td>Polarization angle in band</td><td>5.7 × 10<sup>-6</sup></td><td>E → B</td><td>R</td></tr><tr><td rowspan="3">Gain</td><td>Relative gain in time</td><td>5.7 × 10<sup>-6</sup></td><td>E → B</td><td>R</td></tr><tr><td>Relative gain in detectors</td><td>5.7 × 10<sup>-6</sup></td><td>E → B</td><td>R</td></tr><tr><td>Absolute gain</td><td>1.9 × 10<sup>-6</sup></td><td>B → B</td><td>E</td></tr><tr><td rowspan="4">Polarization angle</td><td>Absolute angle</td><td>9.1 × 10<sup>-6</sup></td><td>E → B</td><td>E</td></tr><tr><td>Relative angle</td><td>5.7 × 10<sup>-6</sup></td><td>E → B</td><td>E</td></tr><tr><td>HWP position</td><td>1.0 × 10<sup>-6</sup></td><td>E → B</td><td>E</td></tr><tr><td>Time variation</td><td>&lt;10<sup>-7</sup></td><td>E → B</td><td>E</td></tr><tr><td>Pol. efficiency</td><td>Efficiency</td><td>5.6 × 10<sup>-6</sup></td><td>B → B</td><td>E</td></tr><tr><td rowspan="3">Pointing</td><td>Offset</td><td>5.7 × 10<sup>-6</sup></td><td>E → B</td><td>R</td></tr><tr><td>Time variation</td><td>&lt;10<sup>-6</sup></td><td>E → B</td><td>E</td></tr><tr><td>HWP wedge</td><td>5.7 × 10<sup>-6</sup></td><td>E → B</td><td>R</td></tr><tr><td>Bandpass</td><td>Bandpass efficiency</td><td>5.3 × 10<sup>-6</sup></td><td>E → B</td><td>R</td></tr><tr><td rowspan="2">Transfer function</td><td>Crosstalk</td><td>5.7 × 10<sup>-6</sup></td><td>B → B</td><td>R</td></tr><tr><td>Detector time constant knowledge</td><td>5.7 × 10<sup>-6</sup></td><td>E → B</td><td>R</td></tr></tbody></table>

knowledge of the beam. We further divide the beam into three: the main lobe; the near sidelobe; and the far sidelobe. The main lobe is the center of the beam with width shown in Table 13. We define the near sidelobe as the region out to 3° with respect to the main-lobe peak direction. The remaining part of the beam is defined as the far sidelobes. The optical system, consisting of the HWP and lenses, is supposed to produce multiple reflections and hence ghosting images for bright sources. The requirement for the suppression of these images is obtained given the required systematic bias in $Δr$. The systematic effects of the beam variation and polarization (co-polar and cross-polar) in the observation band are also studied. For far-sidelobe effects, which give the most sizeable systematic bias, we study the effects for all frequency channels of the three telescopes. We study the HWP systematic effects in terms of the Mueller matrix, consisting of: the instrumental polarization that leaks from temperature into the polarization
signal in the measurement band; the transparency; the polarization efficiency due to the re-
tardance (phase difference) uncertainties of the HWP; and the polarization angle in multiple
bands due to mixing of the two orthogonal polarized signals. We study cosmic-ray effects based
on a newly developed simulation tool. Gain-variation systematic effects are also studied. Given
the systematic error budget, the requirements are obtained on the focal-plane temperature sta-
bility and the gain variation in detectors. Polarization angle systematic effects are sourced by
the uncertainties of the absolute angle, the relative angles among all bands, the HWP spinning
position, and the time variation due to the attitude and HWP position measurement errors.
Polarization efficiency systematic effects for the detector plane antenna and the optical system
are additionally studied. The pointing uncertainty is modeled as a constant offset, as well as a
random variation. The requirements on the HWP wedge angle and the bandpass calibration
uncertainties are also given. The transfer-function category includes systematic effects due to
the uncertainties of the detector time constant and the crosstalk among detectors. Details of
each of these items are given in the following sections. Some of the systematic error budget of
individual sources can be further divided into the sub-level requirements. One example is the
pointing uncertainty, which will be detailed in Sect. 5.3.7.<fnref n="20" />

We note that the estimated systematic bias is approximately a linear sum of individual sys-
tematic effect biases, when those power spectra are much smaller than that of the lensing effect.
This is because Δr is approximately proportional to the square of the error of the systematic
effects to first order; the details are given in Sect. 5.4. The total error budget assigned for the
systematic bias Δr is 0.57 × 10⁻³ [290]. We assign an error budget of 1% of the total systematic
bias 5.7 × 10⁻⁶ to individual systematic effects as a nominal error budget in order to account
for dozens of potential systematic sources. We note that there are some exceptions including
the systematic effects caused by the far-sidelobe uncertainty, however.

To give Δr, we employ one of two approaches. One way is to estimate Δr values based on the
expected uncertainties in the calibration in flight or at the expected precision of measurement
devices. The other approach is to impose requirements on the calibration accuracy or of the
measurement devices. In Table 15, we divide the systematic sources into two types by assigning
“E” (Expectation) and “R” (Requirement), respectively. We show the obtained angular power
spectra for systematic effects in Sect. 5.4.

Our analysis has shown that the employment of the PMU is crucially important to mitigate
systematic errors. Without the PMU, data from a pair of detectors mutually orthogonal in
the polarization orientation are usually combined, causing leakage from temperature to po-
larization when there are differences in any features between the two detectors. The features
giving outstanding systematic effects include the gain, beam, and bandpass. The systematic
errors for these are evaluated assuming that the effects could be mitigated by increasing the
number of detectors *N* (if the systematic effects on individual detectors are uncorrelated), with
the scaling Δ*r* ∝ 1/*N*. Existing studies include: bandpass effects [276], yielding Δ*r* on the order
of 10<sup>-4</sup> with an assumption that the bandpass boundaries have an uncertainty of 1% without
any corrections;<fnref n="21" /> the instrumental polarization caused by reflections on the two mirrors in the
LFT, giving at most 3 × 10<sup>-4</sup> for *M*<sub>QI</sub> components in the Mueller matrix [180]; and the carbon
monoxide (CO) line emission contamination, with Δ*r* ≃ 10<sup>-4</sup> for a bandpass response mismatch
of 10% without applying notch filters [292]. We conduct a simulation study of the systematic

<footnote n="20"><sup>20</sup>We identify about 70 systematic sources in total including the ones in the sub-levels.</footnote>

<footnote n="21"><sup>21</sup>Reference [291] demonstrates the correction of bandpass systematic effects.</footnote>
errors caused by differences for a pair of detectors in gain, pointing, beam ellipticity, beam width, cross-polarization, and sidelobes, and find the systematic bias $\Delta r$ to be on the order of $10^{-6}$. However, these effects would increase significantly when there is some cross-correlation among the detectors with a correlation coefficient larger than $10^{-3}$. A time-correlated noise contribution is expected to contaminate detector timestream data at long timescales. If the noise contamination is projected onto the sky map without adequate treatment, it can cause significant degradation of the sensitivity for the lower-multipole range. We generate 1/f-noise simulated data, project it onto the sky map, and estimate the degradation to be $\delta r \simeq 3 \times 10^{-4}$ for a 20-mHz knee frequency. All the effects above are found to be significantly reduced when we use the PMU. For example, Ref. [276] demonstrates this for the bandpass effects. In summary, with the PMU we can suppress almost all outstanding systematic sources that give rise to the leakage from temperature to polarization, as well as the 1/f noise that causes significant deterioration of the sensitivity in the lower-multipole region.

In the rest of this section, we focus on systematic effects with the PMU. Among those we will give details of systematic studies of the beam, cosmic rays, HWP, and gain in the following sections. Other systematic effects are summarized in Sect. 5.3.7. We estimate the systematic effects based on simulated sky maps of the CMB and foregrounds described in Sect. 5.2.2. The basic strategy is to obtain a residual sky map that is a difference between the reconstructed sky maps with and without a particular systematic effect, which mitigates the model dependence of the sky maps employed.

### 5.3.2. Formalism of the errors.
This section gives the definition of the bias $\Delta r$ and the total error $\delta r$. We define a likelihood function $L(r)$ expressed in the multipole domain due to the large sky coverage of the satellite mission. Assuming Gaussian stationary and isotropic fields and hence no coupling between multipoles (which is an approximation in the presence of a mask of the Galactic plane), the likelihood function is given by

$$ \log L(r) = \sum_{\ell=\ell_{\min}}^{\ell_{\max}} \log P_{\ell}(r), \qquad (36) $$

where $\ell_{\min} = 2$ and $\ell_{\max} = 191$.<fnref n="22" /> We define

$$ \log P_{\ell}(r) = -f_{\text{sky}} \frac{2\ell + 1}{2} \left[ \frac{\hat{C}_{\ell}}{C_{\ell}} + \log C_{\ell} - \frac{2\ell - 1}{2\ell + 1} \log \hat{C}_{\ell} \right] \qquad (37) $$

where $\hat{C}_{\ell}$ ($C_{\ell}$) is the measured (modeled) B-mode power spectrum [293]. We use a sky mask having $f_{\text{sky}} = 0.495$. In order to estimate the potential bias of individual systematic effects on $r$, which we call $\Delta r$, we represent the measured B-mode spectrum as a sum of the following contributions (assuming no primordial B-modes, i.e., $r = 0$):

$$ \hat{C}_{\ell} = C_{\ell}^{\text{sys}} + C_{\ell}^{\text{lens}} + N_{\ell}. \qquad (38) $$

Here $C_{\ell}^{\text{sys}}$ is the estimated systematic effects power spectrum, $C_{\ell}^{\text{lens}}$ is the lensing B-mode power spectrum, and $N_{\ell}$ is the expected noise power spectrum, including the residual noise after component separation (shown in Fig. 39) as well as the cosmic-ray contribution. The modeled power spectrum is given by

$$ C_{\ell} = rC_{\ell}^{\text{tens}} + C_{\ell}^{\text{lens}} + N_{\ell}, \qquad (39) $$

<footnote n="22">²²We use $N_{\text{side}} = 64$, giving $\ell_{\max} = 3 \times N_{\text{side}} - 1$.</footnote>
where $C_{\ell}^{\text{tens}}$ is the tensor mode with $r=1$. We estimate the systematic bias under the condition of the existence of the expected noise $N_{\ell}$ coming from instrumental and foreground residual uncertainties.

The potential systematic bias $\Delta r$ is defined as the value giving the maximum of the likelihood function:

$$ \left. \frac{dL(r)}{dr} \right|_{r=\Delta r} = 0. \qquad (40) $$

We also define the total error on $r$, $\delta r$, as the value covering 68% of the area of the total likelihood function:

$$ \frac{\int_{0}^{\delta r} L(r) dr}{\int_{0}^{\infty} L(r) dr} = 0.68. \qquad (41) $$

### 5.3.3. Beam systematic effects.
Imperfect knowledge of the beam, due to a combination of uncertainties in beam measurements and modeling, are important sources of systematic effects for the measurement of primordial B-modes. We model the sky signal $p(t, v, \hat{n}, \hat{s})$ observed at location $\hat{n}$ at time $t$ with frequency $v$, from a beam centered at location $\hat{s}$ in the presence of an ideal HWP, without coupling with other optical elements:

$$ p(t, v, \hat{n}, \hat{s}) = B_{\parallel}(t, v, \hat{n} - \hat{s})[I(v, \hat{n}) + Q(v, \hat{n}) \cos 2\varphi(t, v) + U(v, \hat{n}) \sin 2\varphi(t, v)] \\ + B_{\perp}(t, v, \hat{n} - \hat{s})[I(v, \hat{n}) - Q(v, \hat{n}) \cos 2\varphi(t, v) - U(v, \hat{n}) \sin 2\varphi(t, v)], \quad (42) $$

where $\varphi(t, v) = 4\rho(t, v) - 2\psi(t, v)$, $\rho(t, v)$ is the HWP rotation position, $\psi(t, v)$ is the polarizer measurement angle with respect to the axis fixed in the sky, and $B_{\parallel}(t, v, \hat{n} - \hat{s})$ and $B_{\perp}(t, v, \hat{n} - \hat{s})$ are the co-polar and cross-polar beams, respectively. The measured signal is the integral of $p(t, v, \hat{n}, \hat{s})$ over the direction $\hat{n}$ and over the frequency $v$ weighted by the bandpass function. Different regions of the beam kernel have different effects on the data and are also estimated differently during the calibration process (either in flight or from the ground). We divide the beam into three regions (as shown in Fig. 32): the main lobe; the near sidelobes; and the far sidelobes. Below we describe systematic effects for the far sidelobes, the near sidelobes, and the main lobe in the order that they are sourced by the uncertainties in our knowledge of $B_{\parallel}$. We also discuss the possible ghosting effect caused by reflections inside the telescope, and the polarization and beam shape in the observed frequency band.

#### 5.3.3.1 Far-sidelobe systematic effects

The main effect of uncertainties in the beam far-sidelobe (FSL) measurement or modeling is that it induces unexpected leakage of the Galactic signal coming from the Galactic plane to higher Galactic latitudes. The mismatch of this excess signal between different frequency bands leads to residuals in the recovered CMB Stokes parameter maps after component separation. With the help of the HWP, allowing quasi-instantaneous measurement of the polarization Stokes parameters ($Q$ and $U$) and of the symmetrization of the effective beam on the maps by the diversity of scanning angle in each pixel of the sky, the dominant systematic effect originates from residual B-mode Galactic signal leaking towards the CMB B-mode estimate. On top of this effect, uncorrected effective beam asymmetries might potentially lead to small E-to-B leakage effects.

Diffraction by optics elements located between the HWP and the detectors is not expected to induce a significant contribution to the cross-polarization beam $B_{\perp}(t, v, \hat{n} - \hat{s})$. This is because
the HWP does not modulate potential generated polarization signal and hence negligible levels of I-to-P leakage should be induced. However, potential contributions to $B_{\perp}(t, v, \hat{n} - \hat{s})$ are the HWP itself and external elements such as baffles. These might induce I-to-P leakage, which will be further studied in subsequent works.

In this subsection we evaluate the impact on the tensor-to-scalar ratio $r$ of the uncertainties on the beam knowledge by perturbing the shape of the beam response at different angles from the beam center. Before describing this procedure we will introduce the simulations that were carried out.

### 5.3.3.2 Description of simulations

For this study, we employ the beam-convolved sky maps generated by TOAST described in Sect. 5.2.2. In order to set the requirements on the far-sidelobe beam calibration, we have considered the beam convolution in the angle ranges $5^{\circ} < \theta < 10^{\circ}$, $10^{\circ} < \theta < 15^{\circ}$, and $15^{\circ} < \theta < 90^{\circ}$, where $\theta$ is the angle with respect to the main-lobe peak. This convolution has been performed for all sets of co-added foreground maps. We do not include noise in the input simulated maps because we are interested in the estimation of the requirements in terms of a bias on $r$ ($\Delta r$) for the beam systematic effects. The time-ordered data (TOD) are simulated for the 22 nominal frequency bands of LiteBIRD. A subset of detectors for each frequency band is considered, instead of the full focal plane. The selected detectors are well spread across the arrays so that we can assume that the results are representative of the full focal plane, and therefore a global rescaling of the results is applied at the end.

LiteBIRD's scanning strategy with nominal parameters of the IMo (see Sect. 5.1) is simulated and a full-sky convolution of the beams is performed. Input simulated reference beams are calculated using GRASP software with a simplified optical model; i.e., no multiple reflections have been included in the beam model. GRASP beams have been simulated at the central frequency for each band, neglecting the color effects and hence the differential effects on various sky components. We will describe the possible systematic effects that are not accounted for in the model above in the last paragraph of this subsection.

### 5.3.3.3 Impact of FSL systematic effects on B-modes

We now evaluate the impact of the lack of knowledge of the beam shape at large angular scales on the measurement of the tensor-to-scalar ratio $r$. With this goal in mind, we process the previously described simulated band-averaged maps produced with non-Gaussian beam convolution through all the steps leading to the measurement of $r$. The effect of uncertainties on the beam shape will then be translated into a bias $\Delta r$ on $r$.

For each perturbed case studied, we proceed with two sets of maps, with and without beam perturbations, which are both processed with component separation to calculate CMB $I$, $Q$, $U$ maps. The first set, including the unperturbed beam, provides reference CMB maps $m_0$. The second set of simulations includes perturbations of the FSL and are run through the same component-separation code with the same assumptions. The perturbations are performed in a specified angle range and for a given frequency band in the following manner:

$$B_{\text{pert}}(v, \theta, \phi) = \mu(\alpha)[1 + \alpha W(\theta)] B_0(v, \theta, \phi). \quad (43)$$

Here $B_0(v, \theta, \phi)$ is the reference GRASP beam for the frequency band centered at $v$, at radial and azimuthal angles $\theta$, $\phi$, $\alpha$ is an amplitude parameter of the perturbation that we are varying to provide the requirements, and $\mu(\alpha)$ is a normalization parameter. The window function $W(\theta)$
is chosen to include apodization in order to reduce ringing edge effects and is set to

$$W(\theta) = A_{\text{inf}} - A_{\text{sup}}, \quad \text{where} \quad A_i(\theta) = \begin{cases} 0, & \text{if } \theta < \frac{1}{2}\theta_i \\ \frac{1}{2}\left(1 - \cos\left(\frac{(2\theta - \theta_i)\pi}{\theta_i}\right)\right), & \text{if } \frac{1}{2}\theta_i < \theta < \theta_i \\ 1, & \text{if } \theta_i < \theta, \end{cases} \quad (44)$$

where $i = (\text{inf, sup})$, for three different angle ranges [$\theta_{\text{inf}}$, $\theta_{\text{sup}}$], namely [5°, 10°], [10°, 15°], and [15°, -] ($A_{\text{sup}} = 0$ for the third window). With the previous parameters, the window is actually effective, i.e., $W/W_{\text{max}} > 0.5$, in the angle ranges [3.75°, 7.5°], [7.21°, 11.69°], and [11.25°, 90°].

The normalization $\mu$ is calculated for each value of $\alpha$ such that the beam-convolved maps have the same dipole component ($\ell = 1$) amplitude as the unperturbed ones. This procedure allows us to isolate the effect of beam mismatch only from calibration uncertainties that are studied in Sect. 4.8.2.

In our approach, we find the limiting values $\alpha_{\text{lim}}$ independently for each frequency band; i.e., we perturb one beam at one frequency at a time, leaving the others constant, such that the induced bias on the tensor-to-scalar ratio is $\Delta r = 1.9 \times 10^{-5}$ divided by 66, the number of frequency channels times the number of windows. The systematic bias $\Delta r$ is obtained using Eq. (40) given $C_{\ell}^{\text{sys}}$, which is obtained from the residual sky maps $\delta \mathbf{m} = \mathbf{m} - \mathbf{m}_0$ after applying a Planck-HFI Galactic mask leaving 51% of sky area, where $\mathbf{m}$ ($\mathbf{m}_0$) is the sky map convolved by the far sidelobe with (without) the perturbation of individual $\alpha$ values.

The results for the three angle ranges are shown in Table 16. The $\alpha$ are just intermediate parameters and we indicate the values of the beam perturbation amplitudes in the three windows by calculating the following values:

$$\delta \bar{B}_{\text{lim}} = \frac{\int \delta B_{\nu}(\theta, \phi) W(\theta) d\Omega}{\int W(\theta) d\Omega} \quad \text{and} \quad \delta R_{\text{lim}} = \frac{\int \delta B_{\nu}(\theta, \phi) W(\theta) d\Omega}{\int B_{0}(\nu, \theta, \phi) d\Omega}. \quad (45)$$

Here the beam perturbation $\delta B_{\nu}(\theta, \phi) = \alpha_{\text{lim}}B_0(\nu, \theta, \phi)$, including the normalization of the input beam at the peak: $B_0(\nu, \theta = 0) = 1$. The quoted values $\delta \bar{B}_{\text{lim}}$ represent the precision that is required on the knowledge of the mean amplitude of the beam in the window $W(\theta)$, while $\delta R_{\text{lim}}$ represents the required precision on the beam's relative power in the window angle range specified by $W(\theta)$. Because the amplitude of the beam models used for this analysis drops drastically at angles larger than around 50°, we have cut the last window for the perturbation $W(\theta)$ above 70° (to be conservative) for the calculation of $\delta \bar{B}_{\text{lim}}$ and $\delta R_{\text{lim}}$, as well as the estimated calibration precision that we will discuss later.

The derived constraints can be translated into the required accuracy in beam-calibration measurements during ground testing. The quantity effectively measured during the beam calibration can be modeled as

$$P_{\text{cal}}(\vec{r}) = \int B_{\nu} \omega(\vec{r}' - \vec{r}) d\Omega' \frac{1}{\int B_{\nu} \omega(\vec{r}') d\Omega'} + n_{\text{cal}}, \quad (46)$$

with $\omega(\vec{r})$ a small integration window of the beam and $n_{\text{cal}}$ the noise in the beam calibration. We have estimated the precision required on the measured quantity $P_{\text{cal}}(\vec{r})$ assuming random uncorrelated errors in each measurement (and hence no systematic effects in the calibration), and a grid of measurements at many angles. In this case the calibration measurement uncertainty $\sigma_{\Omega_{\text{pix}}^{\text{cal}}} = \sqrt{\langle n_{\text{cal}}^{\text{T}} n_{\text{cal}} \rangle}$ is related to the uncertainty in the beam amplitude averaged over the
<table><caption>**Table 16.** Beam perturbation requirements for each frequency channel and each of the three angle ranges of the beam perturbations giving $\Delta r = 1.9 \times 10^{-5}/66$, fixing the other perturbations for the other frequency channels and angle ranges to 0. Units are dB for $\delta \bar{B}_{\text{lim}}$ and $\sigma_{0.25\square}^{\text{cal}}$, and dBi for $\delta R_{\text{lim}}$. The indicated values of angles in the first row, providing the typical range in which the beam is perturbed, are such that 95% of the perturbed beam power for a typical CMB channel is within this range (the value varies from channel to channel as the beam shape varies). This explains why the upper bound for the last window is only 50° while the window extends to 90° as the beam drops above 50–60°, depending on the channel. The calculated values of the required precision $\sigma_{0.25\square}^{\text{cal}}$ are assuming 0.5° separation between noise-limited measurements. An identical calibration precision of $\sigma_{0.25\square}^{\text{cal}} = -56.90$ dB for all 66 entries leads to a similar bias.</caption><thead><tr><th rowspan="2"></th><th colspan="3">3° &lt; θ &lt; 8°</th><th colspan="3">6° &lt; θ &lt; 13°</th><th colspan="3">10° &lt; θ &lt; 50°</th></tr><tr><th>δB̄<sub>lim</sub></th><th>δR<sub>lim</sub></th><th>σ<sup>cal</sup><sub>0.25□</sub></th><th>δB̄<sub>lim</sub></th><th>δR<sub>lim</sub></th><th>σ<sup>cal</sup><sub>0.25□</sub></th><th>δB̄<sub>lim</sub></th><th>δR<sub>lim</sub></th><th>σ<sup>cal</sup><sub>0.25□</sub></th></tr></thead><tbody><tr><td>L1-040</td><td>-42.55</td><td>-23.54</td><td>-28.20</td><td>-46.62</td><td>-25.41</td><td>-30.68</td><td>-66.40</td><td>-27.49</td><td>-42.79</td></tr><tr><td>L2-050</td><td>-34.09</td><td>-13.45</td><td>-19.73</td><td>-38.35</td><td>-15.51</td><td>-22.41</td><td>-57.98</td><td>-17.44</td><td>-34.37</td></tr><tr><td>L1-060</td><td>-39.46</td><td>-17.68</td><td>-25.11</td><td>-43.67</td><td>-19.68</td><td>-27.73</td><td>-63.20</td><td>-21.53</td><td>-39.60</td></tr><tr><td>L3-068</td><td>-36.70</td><td>-13.27</td><td>-22.34</td><td>-41.38</td><td>-15.75</td><td>-25.44</td><td>-61.52</td><td>-18.20</td><td>-37.92</td></tr><tr><td>L2-068</td><td>-30.45</td><td>-8.02</td><td>-16.10</td><td>-35.62</td><td>-10.99</td><td>-19.68</td><td>-56.00</td><td>-13.67</td><td>-32.39</td></tr><tr><td>L4-078</td><td>-42.82</td><td>-18.36</td><td>-28.46</td><td>-45.87</td><td>-19.21</td><td>-29.93</td><td>-61.76</td><td>-17.41</td><td>-38.16</td></tr><tr><td>L1-078</td><td>-39.42</td><td>-16.43</td><td>-25.07</td><td>-42.48</td><td>-17.28</td><td>-26.53</td><td>-59.07</td><td>-16.18</td><td>-35.46</td></tr><tr><td>L3-089</td><td>-48.56</td><td>-23.14</td><td>-34.20</td><td>-51.75</td><td>-24.13</td><td>-35.81</td><td>-68.70</td><td>-23.40</td><td>-45.10</td></tr><tr><td>L2-089</td><td>-38.84</td><td>-15.50</td><td>-24.49</td><td>-42.04</td><td>-16.49</td><td>-26.10</td><td>-59.51</td><td>-16.27</td><td>-35.91</td></tr><tr><td>L4-100</td><td>-51.80</td><td>-25.61</td><td>-37.45</td><td>-54.98</td><td>-26.59</td><td>-39.04</td><td>-72.09</td><td>-26.01</td><td>-48.49</td></tr><tr><td>L3-119</td><td>-54.82</td><td>-27.57</td><td>-40.47</td><td>-57.91</td><td>-28.46</td><td>-41.97</td><td>-75.00</td><td>-27.86</td><td>-51.39</td></tr><tr><td>L4-140</td><td>-51.25</td><td>-23.23</td><td>-36.90</td><td>-54.11</td><td>-23.88</td><td>-38.17</td><td>-74.60</td><td>-26.69</td><td>-51.00</td></tr><tr><td>M1-100</td><td>-50.65</td><td>-26.84</td><td>-36.30</td><td>-53.91</td><td>-27.90</td><td>-37.97</td><td>-69.44</td><td>-25.74</td><td>-45.84</td></tr><tr><td>M2-119</td><td>-54.58</td><td>-29.71</td><td>-40.23</td><td>-57.78</td><td>-30.71</td><td>-41.84</td><td>-73.96</td><td>-29.20</td><td>-50.36</td></tr><tr><td>M1-140</td><td>-49.55</td><td>-23.92</td><td>-35.20</td><td>-52.45</td><td>-24.61</td><td>-36.51</td><td>-75.36</td><td>-29.83</td><td>-51.76</td></tr><tr><td>M2-166</td><td>-60.87</td><td>-34.64</td><td>-46.52</td><td>-64.07</td><td>-35.63</td><td>-48.13</td><td>-80.06</td><td>-33.91</td><td>-56.45</td></tr><tr><td>M1-195</td><td>-63.56</td><td>-37.07</td><td>-49.20</td><td>-66.58</td><td>-37.89</td><td>-50.64</td><td>-80.69</td><td>-34.30</td><td>-57.08</td></tr><tr><td>H1-195</td><td>-59.01</td><td>-33.00</td><td>-44.65</td><td>-62.38</td><td>-34.17</td><td>-46.44</td><td>-78.63</td><td>-32.73</td><td>-55.03</td></tr><tr><td>H2-235</td><td>-62.92</td><td>-35.65</td><td>-48.57</td><td>-66.07</td><td>-36.60</td><td>-50.13</td><td>-81.55</td><td>-34.38</td><td>-57.95</td></tr><tr><td>H1-280</td><td>-60.57</td><td>-32.52</td><td>-46.22</td><td>-63.60</td><td>-33.34</td><td>-47.66</td><td>-75.23</td><td>-27.26</td><td>-51.63</td></tr><tr><td>H2-337</td><td>-70.48</td><td>-41.78</td><td>-56.12</td><td>-73.57</td><td>-42.66</td><td>-57.63</td><td>-86.41</td><td>-37.79</td><td>-62.81</td></tr><tr><td>H3-402</td><td>-67.91</td><td>-38.03</td><td>-53.56</td><td>-70.77</td><td>-38.68</td><td>-54.83</td><td>-88.72</td><td>-38.93</td><td>-65.12</td></tr></tbody></table>

window area $\delta \bar{B}_{\text{lim}}$ by

$$ \sigma_{\Omega_{\text{pix}}}^{\text{cal}} = \frac{\int W(\theta) d\Omega}{\sqrt{\sum_{ij} W^2(\theta_{ij}) \Delta\Omega_{\text{pix}}}} \delta \bar{B}_{\text{lim}}, \quad (47) $$

with $\Delta\Omega_{\text{pix}}$ the solid angle covered by one calibration measurement and $i, j$ the pixel numbers for a pixelized beam-calibration map. We assume that the calibration measurements are normalized to unity at the peak (for the beam center). From the previous equation, we can see that the factor $(\int W d\Omega / \sqrt{\sum_{ij} W^2(\theta_{ij}) \Delta\Omega_{\text{pix}}})^2$ is the effective number of pixels in the window area of the beam perturbation. Assuming a pixel width $\sqrt{\Delta\Omega_{\text{pix}}}$ = 0.5°, which is the expected beam measurement step size in the ground-calibration process, we have calculated the required precision for the calibration measurements and present the values in Table 16. The effective numbers of measurements for the three windows are 742.0, 1541.9, and 52943.2.
We can see that the requirements on calibration precision ($\sigma_{0.25\square}^{cal}$) are more stringent at the highest frequencies for the last window, which covers a larger area of the sphere. This shows that it is necessary to accurately calibrate the FSLs of the channels, providing a proxy for the Galactic dust for component separation and for which the fraction of power going into the main lobe is lower, since those have the smallest beam FWHM values. Assuming that we calibrate the whole area of the beam for all the frequency channels with a constant precision $\sigma_{0.25\square}^{cal}$, we have evaluated the desired precision as $-56.90$ dB to have requirements equivalent to that presented in Table 16. The total contribution of $\Delta r$ has been estimated with 100 realizations by simultaneously varying individual $\alpha$ values with a uniform distribution in the range $[\max(-1, -\sqrt{3}\alpha_{lim}); \sqrt{3}\alpha_{lim}]$, with $\alpha_{lim}$ calculated from Table 16. The $\sqrt{3}$ factor is to keep the same variance $\alpha_{lim}^2$ as the Gaussian distribution used for Table 16 to derive the requirements. We find the systematic bias resulting from the above assumptions to be $\Delta r = 2.2 \times 10^{-5}$.

We want to stress several important limitations of this analysis. Firstly, one of the key assumptions is that the scanning strategy induces enough symmetrization of the beam so that many calibration measurements can be averaged for large $\theta$ and that the beam systematic effects are then dominated by uncertainties in the transfer function. This explains why the number of effective calibration points is high for each window. Secondly, in order to be representative of the final results, we have assumed that the bias on $r$ is not too sensitive to the structure of the beam inside each window $W(\theta)$. Thirdly, the correlation in the overlapping region in Table 16 is assumed to be negligible. Lastly, we quote the required precision on the calibration measurements assuming that they are limited by calibration noise. This is somewhat unrealistic, since measurements are expected to be limited by systematic effects in the ground measurements. Nevertheless, this analysis still provides a framework for quantifying uncertainties in further studies that might be needed.

We have checked that the results are not sensitive to the component-separation parameters and the model, in particular the effective beam function assumed as input in Sect. 5.2. The variations introduced by the change of the recovered component-separation parameters due to the slight difference of the beam transfer functions in the recovered CMB maps vanish at first order in the difference $\mathbf{m} - \mathbf{m}_0$.

The beam model used for the simulations contains several uncertainties including the effects of beam asymmetry caused by the multiple reflections and scattering, especially for the detectors placed at the edge of the focal plane. In order to account for the uncertainties, we conservatively multiply by a factor of two the requirements to obtain the FSL total systematic bias of $\Delta r = 4.4 \times 10^{-5}$.

#### 5.3.3.4 Near-sidelobe effect

We pursue two approaches to estimate the systematic effects due to inaccuracy of the knowledge of the near sidelobe, i.e., the beam structures in the angle range up to about $3^\circ$ with respect to the main-lobe center. In the first method, we model the beam shape with an axially symmetric Gaussian main beam plus a power-law tail. We convolve the sky (including foregrounds) with the beam and deconvolve the sky map with a beam with a different amplitude of the power-law tail. Then we take the difference between the sky maps before and after this procedure for the 15 frequency channels, and apply the component-separation procedure to obtain the residual sky map and power spectrum. We find that $\pm 10\%$ variations in the amplitude of the near sidelobe, whose nominal magnitude is less than 1% of the main lobe, give rise to a systematic bias of $\Delta r$
$= 5.7 \times 10^{-6}$. In the second method, we model the beam as the sum of the Gaussian main beam and a ring shape near the sidelobe, whose cross-sectional shape is a Gaussian with a peak height of 1% of the main-lobe peak. The position, width, and height of the modeled near sidelobes are set to values close to those in the GRASP simulation. The typical values are 1.7° from the main-lobe center for the position, 0.6° for FWHM, and -20 dB with respect to the main lobe for the height. We repeat a similar analysis by changing the ring position, width, and height, and find that a ±20% variation from the nominal values corresponds to $\Delta r = 5.7 \times 10^{-6}$. In summary, we require a calibration accuracy of better than 10% for the near sidelobe in order to have a systematic bias $\Delta r < 5.7 \times 10^{-6}$.

### 5.3.3.5 Main-lobe effect

At the simplest level, the beam's main lobe is characterized by two quantities, namely the height and the shape. The height corresponds to the gain, which we describe in Sect. 5.3.6. The main-lobe beam has a shape close to elliptical for the detectors located far from the center of the focal plane. We quantify this effect using two parameters, namely the width and the ellipticity or flattening. We conduct a simulation study similar to the one for the near-sidelobe effects. With an expected calibration accuracy of 1% for both the width and the flattening in flight (see Sect. 4.8.2), we find that the resulting systematic bias is smaller than $10^{-6}$.

### 5.3.3.6 Ghosting effect

It is known that the multiple reflections between the focal plane and the optics components, such as lenses and HWP, produce a ghosting image [e.g., Ref. 294]. To investigate this, we conduct a simulation study to estimate the systematic effects due to foreground contamination leaking through off-boresight small-scale structure into the beam. We find that a leakage of 0.05% amplitude through a 30 arcmin diameter spot separated from the beam boresight direction by 900 arcmin produces a systematic bias equivalent to $\Delta r = 5.7 \times 10^{-6}$. The lens reflection effects are included in the discussion regarding the FSL. We note that the ghosting effect can be partly mitigated by tilting the HWP by an angle larger than half of the FoV.

### 5.3.3.7 Polarization and beam shape in band

Sinuous antennas are sensitive to polarization at specific angles due to their design. However, the angle of the polarization sensitivity has a dependence on the frequency of the incoming light (commonly referred to as the “wobbling effect”). In our case, this effect is of the order of 3° in the observing band [272]. It is known that the wobbling effect can be canceled using four detectors, i.e., two pairs of detectors with mutually mirrored sinuous patterns [272]. We study wobble cancelation effects with the beam patterns obtained from a GRASP simulation for 100-GHz detectors with 88, 100, and 112-GHz frequencies, corresponding to the central and edge frequencies in the band. We use the QuickPol [295] algorithm, in which the LiteBIRD scanning strategy is taken into account, to estimate the residual systematic effect when the two mutually mirrored patterned pairs have different polarization efficiencies. We find that a difference of 20 dB produces a residual power spectrum lower by 1% with respect to the one due to the lensing effect, corresponding to $\Delta r < 10^{-6}$. Note that the beam calibration will be conducted with a precision of 53 dB, as discussed for the FSL systematic effects. Considering this assumption, this systematic effect can be considered to be negligible.
5.3.4. Cosmic-ray systematic effects. For Planck-HFI, cosmic rays produced glitch signals at a rate of 5 hits cm⁻² s⁻¹ [296], giving rise to systematic errors and requiring a deglitching process. For LiteBIRD, we will use a different detector technology and expect that the outcome of cosmic-ray effects may differ from that of Planck-HFI. In order to predict and characterize the potential effects of cosmic-ray glitches, two studies have been conducted: one is to use actual laboratory measurements to estimate the response of the detectors to energetic particle impacts by irradiating the detectors with alpha particles, γ-rays, or cosmic-ray muons; and the second is to realize an end-to-end simulation to evaluate the propagation of cosmic-ray effects. The former is described in Sect. 4.4.2 (see also Ref. [225]). Here we describe details for the latter approach (see also Refs. [297,298]).

For LiteBIRD, we will employ TES arrays on a 10-inch silicon wafer. We estimate the expected cosmic-ray hit-rate on a single silicon substrate using PAMELA data fitted with an Usoskin model at L2 for the cosmic-ray flux determination [299] and GEANT4 for simulating the interaction of particles (mainly protons) in the substrate and the surrounding materials (mainly aluminum), including secondaries and electromagnetic showers. The cosmic-ray rate is expected to be 400 Hz on a single silicon wafer on which the bolometers are placed. Since the time between hits is comparable to the time constant of the TES (a few ms), individual cosmic-ray hits may not be identified. Hence the net effect is an increase in noise, which may be non-Gaussian.

Cosmic rays impacting the silicon wafer deposit part of their energy of 1.8 MeV on average as heat. The heat propagation in the wafer is simulated using a finite-element thermal model in the commercial software COMSOL, assuming that the wafer heat capacitance and the heat conductance to the refrigerator's thermal bath is kept at a temperature of 100 mK [298]. The heat produced by the cosmic-ray impact propagates to the TES through the SiN support structure, producing a transient fluctuation in the TES current. The sensitivity of the detector is used to convert the TES current to an equivalent signal power. The bolometer analog signal is digitized with a rate of 20 MHz and down-sampled to about 19 Hz with the application of digital low-pass filters. In this process, the higher-frequency components are dropped without any aliasing effects. The bolometers' signals are read out in the frequency domain using cryogenic resonators. The resonant peak overlap in the frequency domain may cause crosstalk between detectors, since the cosmic-ray glitch signal is fast and has a wide bandwidth. In the simulation study, we implement both effects in the TOAST simulator that generates the time-ordered data and projects them onto sky maps [297].

In this study, we do not include the athermal (ballistic) phonon effect, since this effect is hard to model at this stage. The modeling would require detailed comparisons between the measurements of the athermal phonon signals and the underlying physics models in a simulation. We plan to develop a ray-tracing technique for producing Monte Carlo simulations of athermal phonon propagation effects in the LiteBIRD detector wafer, which are currently ongoing. It should be noted, however, that with the expected impact rate of 400 hits s⁻¹ in the wafer, the total amount of energy that propagates as athermal phonons will contribute to the same amount of average noise across the wafer, as indicated by modeling the thermal phonons. Furthermore, we note that athermal phonons would produce short pulses smoothed by the TES time constant of a few ms, whose higher-frequency components are removed by the digital filters.
In Sect. 5.4 we show the band-weighted averaged angular power spectrum obtained with the technique developed in Refs. [297] and [298]. We take into account this effect as an additional noise in $N_ℓ$ in Eqs. (38) and (39). The magnitude of the cosmic-ray noise is about 1% of the total noise estimated with the component separation in Fig. 39.

5.3.5. *HWP systematic effects.* The fast rotating HWPs are key elements of LiteBIRD, allowing quasi-instantaneous estimation of the Stokes parameters I, Q, and U with individual polarized detectors. An ideal HWP rotates the incident polarization by a known angle. However, due to the complexity of the structure of transmissive broadband HWPs, requiring multiple stacked sapphire plates with ARC, dielectrically embedded in a multilayer structured plate, etc., the transformation of Stokes parameters is more complex, leading to potential artifacts in polarization maps (specifically, mixing of Stokes parameters). These imperfections can be described by the Jones formalism connecting the input and output electric fields, before and after the HWP, respectively [300]:

$$ \begin{pmatrix} E_{x;0} \\ E_{y;0} \end{pmatrix} = \begin{pmatrix} J_{xx} & J_{xy} \\ J_{yx} & J_{yy} \end{pmatrix} \begin{pmatrix} E_{x;i} \\ E_{y;i} \end{pmatrix}, \qquad (48) $$

as well as the Mueller matrix formalism relating the input and output Stokes parameters:

$$ \begin{pmatrix} I_o \\ Q_o \\ U_o \\ V_o \end{pmatrix} = \begin{pmatrix} M_{II} & M_{IQ} & M_{IU} & M_{IV} \\ M_{QI} & M_{QQ} & M_{QU} & M_{QV} \\ M_{UI} & M_{UQ} & M_{UU} & M_{UV} \\ M_{VI} & M_{VQ} & M_{VU} & M_{VV} \end{pmatrix} \begin{pmatrix} I_i \\ Q_i \\ U_i \\ V_i \end{pmatrix}. \qquad (49) $$

The relation between the Mueller matrix M and the Jones matrix J can be found in, e.g., Refs. [301] and [302].

We now introduce different types of perturbations in the coefficients describing the HWP in order to model imperfections. We have carried out two independent studies of the effect of imperfections on the final B-mode signal: the first study, based on the Mueller matrix formalism, is focused on the effect of instrumental polarization (IP) for a transmissive HWP; and the second study, based on the Jones formalism, addresses the effects of polarization efficiency and Q/U mixing. Other potential systematic effects include the HWP rotation-synchronous signal observed in Refs. [206,209], which will be studied in future.

We first focus on the effect of instrumental polarization. For oblique incidence angles, a small fraction of the unpolarized incident intensity is transferred into polarization (due to anisotropic transmission for a transmissive HWP), inducing non-zero $M_{QI}$ and $M_{UI}$ elements with azimuthal angle dependence in particular at four times the HWP rotation frequency. A model of the observations can be written as

$$ s = (1 \cos 2\psi_0 \sin 2\psi_0 0) M(\Theta, \rho - \psi) R(2\psi) \begin{pmatrix} I_i \\ Q_i \\ U_i \\ V_i \end{pmatrix} + n. \qquad (50) $$

The various quantities here are: s is the observed signal on a polarized detector at angle $\psi_0$ with respect to a reference axis in the focal-plane (FP) coordinates; $\psi$ is the rotation angle of the reference axis of the FP with respect to the sky reference frame in which the Stokes parameters $I_i$, $Q_i$, $U_i$, and $V_i$ of the incoming wave are calculated; $\rho = \omega_{HWP} t + \rho_0$ is the HWP rotation
angle with respect to the sky reference frame; $\Theta$ is the incident angle with respect to the optical axis of the HWP; and $R(2\psi)$ is the rotation matrix for Stokes parameters. The HWP Mueller matrix $M$ is written in the frame of the HWP.

The Mueller matrix of the HWP has been calculated using the electromagnetic wave propagation simulation tool rigorous coupled-wave analysis (RCWA) assuming a nine-layer HWP with pyramidal ARC optimized for LFT.<fnref n="23" /> In an electromagnetic simulation, we find that the Mueller matrix can be expanded into three components consisting of a constant one (0f), one modulated with an angular frequency twice that of the HWP rotation (2f), and one with four times that frequency (4f). Each element of the matrix $M_{ij}$ for a given incident angle $\Theta$ is decomposed into sine waves with respect to the angle $\rho - \psi$, and only the first three terms are shown to be significant:

$$ M_{ij}(\Theta, \rho - \psi) = M_{ij}^{(0f)}(\Theta) + M_{ij}^{(2f)}(\Theta) \cos(2\rho - 2\psi + \phi_{ij}^{(2f)}) \\ + M_{ij}^{(4f)}(\Theta) \cos(4\rho - 4\psi + \phi_{ij}^{(4f)}). \qquad (51) $$

Coefficients of the three matrices $M^{(0f)}$, $M^{(2f)}$, and $M^{(4f)}$, as well as the phases $\phi_{ij}^{(2f)}$ and $\phi_{ij}^{(4f)}$<fnref n="24" /> have been estimated for the whole band at 140 GHz and for many angles $\Theta$. The IP coefficients of the $M^{(4f)}$ matrix (four times the HWP spinning frequency) are the ones inducing a significant effect on the final B-mode power spectrum. The values of those coefficients $\epsilon = M_{QI}^{(4f)} = M_{UI}^{(4f)}$ are found to be $\simeq 4 \times 10^{-5}$ for $\Theta_0 = 10^\circ$. The empirical relation $\epsilon(\Theta) = \epsilon(\Theta_0)(\sin\Theta/\sin\Theta_0)^2$ is used to extrapolate the coefficients for any $\Theta$ values corresponding to various locations of detectors in the focal plane.

We have generated input TOD simulations implementing Eq. (50), which include polarized CMB anisotropies (with $r=0$), CMB dipole, and polarized Galactic foregrounds from PySM described in Sect. 5.2.1, with and without noise.<fnref n="25" /> We assume the nominal scanning strategy described in Table 12. Simulated data are projected into I, Q, U maps using the optimal GLS map-making method SANEPIC [303] (see also Ref. [304] for a similar derivation and a description of the application for mapping on a sphere) applied to each detector. The residual maps are calculated by differencing the recovered and input maps used for the simulations.

Some cancellation of the leakage contributions happens when combining several detectors at different locations around the HWP rotation axis. This is because, for detectors with different azimuthal angles around the HWP rotation axis, the ordinary axis of the HWP is orthogonal to the different detector lines of sight at different times, leading to different phases of the IP contribution for different detectors at a given time. This phase shift is accounted for in our multi-detector simulations. The resulting contribution of the IP to the B-mode spectra after combining all Stokes parameter maps for all detectors at 140 GHz is shown in Fig. 40. The contributions of the dipole and the other fluctuations are shown separately. An apodized Galactic mask has been applied, considering 50% of the sky. We have assumed two cases: (i) the optical axis of the HWP is centered in the LFT focal plane; and (ii) the HWP rotation axis is shifted

<footnote n="23">²³For the HWPs of MHFT, we give a mathematical model using the Jones matrix, which will be described in the following section.</footnote>

<footnote n="24">²⁴Those phases are defined for one reference location in the focal plane and are modified for different azimuthal angles around the HWP rotation axis in a fixed frame in the focal plane.</footnote>

<footnote n="25">²⁵Since the model and the map-making relations are linear with respect to the input TOD, the effect of HWP imperfection on the recovered signal can be estimated without noise in the simulations.</footnote>
<figure><img src="image_42.png" /><figcaption>**Fig. 40.** Left: *B*-mode power spectra of the residual IP leakage. The curves are: the residual power spectrum for one detector at 6.24° from the HWP rotation axis (red); the combined effect (all detectors at 140 GHz LFT) for the dipole only and a centered HWP (cyan); the combined effect of the dipole and 5° HWP tilt (violet); and the same as the cyan and violet cases, but for anisotropies without the dipole (yellow and black). Right: *B*-mode power spectra of residual IP after correction with the fit of individual detector IP amplitude in cyan (model 1), and using the $\epsilon(\Theta)$ model in violet (model 2). The power spectrum without correction (effect of the dipole and the anisotropies combined) is shown in red. We used the less favorable case of a 5° tilt of the HWP rotation axis.</figcaption></figure>

<table><caption>**Table 17.** Contribution of the IP to the $\Delta r$ parameter. The most pessimistic case of a shifted HWP is used for evaluation of the effect after correction with two different schemes: for model 1, when the $\epsilon_i$ coefficients are estimated independently for each detector; and for model 2, when we use the scaling relation $\epsilon(\Theta, \phi)$ accounting for the exact relative phase dependence of this parameter with respect to the azimuthal angle of the detector locations $\phi$ around the HWP rotation axis (two detectors at different locations in the FP see the same 4f IP effect with a known delay, as described further in the text).</caption><thead><tr><th></th><th>HWP centered</th><th>HWP shifted by 5°</th><th>Mitigated, model 1</th><th>Mitigated, model 2</th></tr></thead><tbody><tr><th>Δr</th><td>1.36 × 10<sup>-4</sup></td><td>1.47 × 10<sup>-3</sup></td><td>&lt;1.20 × 10<sup>-6</sup></td><td>&lt;1.84 × 10<sup>-7</sup></td></tr></tbody></table>

by 5°. There is less cancellation in the second case. The induced bias $\Delta r$ on the tensor-to-scalar ratio $r$ using the likelihood described in Sect. 5.3.2 is indicated in Table 17.

A correction method has been implemented and applied to the simulations including noise. With this method, the Stokes parameter maps, as well as the IP coefficients $\epsilon$, are jointly estimated by maximizing a global likelihood. The recovered polarization maps are then corrected from the effect. Two models used are: (1) when the $\epsilon_i$ coefficients are estimated independently for each detector; and (2) when we use the scaling relation $\epsilon(\Theta, \phi)$ with the relative phase dependence of this parameter with respect to the azimuthal angle of the detector locations $\phi$ around the HWP rotation axis. The method does not use any external template and exploits the large-intensity signal to recover the $\epsilon$ parameters. Because the monopole is not recovered within the map-making process, it is jointly fitted with the $\epsilon$ parameters, and so is marginalized over, in order to remove its leakage in the polarization maps. The residual *B*-mode spectra of the IP after correction are shown in Fig. 40 and the values of $\Delta r$ are indicated in Table 17. A three-year period of observation is assumed here. We observe that the bias is greatly reduced after mitigation. We argue that the phase of the residual effect in the map is random after mitigation because the errors on $\epsilon$ are limited by noise, so if more frequency channels are included in the analysis the net effect after component separation should be reduced. The quoted value
in Table 17 calculated for the 140-GHz LFT channel only is then expected to be pessimistic. We do not consider in this study the coupling with other systematic effects, such as gain variation that might affect the estimation of the parameters and the efficiency of the subtraction. This is postponed to future studies. In our second study, we only consider orthogonal incidence on the HWP, i.e., $\Theta = 0$ and HWP systematic effects other than the IP leakage.

Even if orthogonal propagation through the HWP is assumed, manufacturing imperfections could still lead to systematic effects in the observed CMB signal. The goal of our second study is to set requirements on the accuracy needed to constrain departures of the HWP from the ideal setup. Requirements are set by imposing a threshold on the maximum bias $\Delta r$ that we could tolerate being due to a combination of HWP systematic effects, including the effect of frequency-dependent HWP parameters [305,306].

We can make explicit the dependence of the Jones matrix in Eq. (48) (in the HWP frame) on the HWP non-idealities:

$$J_{\text{HWP}} = \begin{pmatrix} 1 + h_1 & \zeta_1 e^{i\chi_1} \\ \zeta_2 e^{i\chi_2} & -(1 + h_2)e^{i\beta} \end{pmatrix}, \qquad (52)$$

where the frequency dependence of each parameter is understood. In Eq. (52), $h_{1,2} < 0$ are loss terms, $\beta = \psi - \pi$, where $\psi$ is the phase shift between the orthogonal modes, and $\zeta_{1,2}$ and $\chi_{1,2}$ are amplitudes and phases responsible for x-y polarization mixing. The phases $\chi_{1,2}$ are set to zero in this work, being degenerate at first order with $\zeta_{1,2}$. From the Jones matrix, the equivalent Mueller matrix elements $M^{IX} = M^{IX}(h_{1,2}, \zeta_{1,2}, \beta)$ with $X = I, Q, U$ can be derived.

In the simulation we consider a pair of orthogonal detectors on the boresight with a noise contribution. Since we want to propagate to $r$, we make use of CMB bands, specifically four MFT bands and one band for LFT and HFT, respectively (the closest ones to the CMB channels). We make use of simulated mesh-HWP profiles in frequency for the parameters ($h_1(\nu)$, $h_2(\nu)$, $\beta(\nu)$) in the four MFT bands centered at 100, 119, 140, and 166 GHz, respectively. We assume $\zeta_{1,2} = 10^{-2}$, constant in frequency [211]. For the LFT/HFT band centered at 100/195 GHz, we assume only parameters that are constant in frequency, since our study is not strongly affected by the frequency dependence of the model profiles. The resolution of those profiles is 1 GHz.

When building the TOD, we perturb the HWP frequency profiles so as to simulate a mismatch with the nominal HWP profile. In particular, we have $M_{\text{TOD}}^{IX} \equiv M^{IX}(h_{12} + \Delta h, \zeta_{1,2} + \Delta \zeta, \beta + \Delta \beta)$, where, again, the frequency dependence in each term is understood. In each frequency bin within a band, the perturbations $\Delta x$ are drawn from band- and frequency-independent Gaussian distribution with variance $\sigma_{\Delta x}^2$, where $x$ stands for $h, \beta$, or $\zeta$. Instead, in the map-making procedure, we make use of the unperturbed, nominal HWP profile, i.e., $M_s^{IX} \equiv M^{IX}(h_{1,2}, \zeta_{1,2}, \beta)$.

In order not to be dependent on a specific realization of the systematic perturbations, we simulate 10 realizations for each $\sigma_{\Delta x}$ by always keeping CMB and foreground fixed. For each band, we compute a template map $m_{\text{templ}}$ with $M_{\text{TOD}}^{IX}(\nu) = M_s^{IX}(\nu)$ (and thus no error on the HWP parameters) from the same foreground and CMB maps as above. These templates are used to obtain residual maps; this yields maps of residuals $m_{\text{res}}$ that are minimally affected by the foreground color effect and mostly due to the mismatch between $M_{\text{TOD}}^{IX}(\nu)$ and $M_s^{IX}(\nu)$.

For each map of residuals, masked with a 70% Galactic mask, the corresponding B-mode spectrum is computed and added to the fiducial CMB B-mode spectrum. We use the simple anafast pseudo-$C_\ell$ for the residuals, since those are not very different from the cases in which
<table><caption>**Table 18.** Accuracy level needed for measurements of HWP parameters *h*, *β*, and *ζ* in order to keep the bias on *r* below Δ*r* = 5.7 × 10⁻⁶. Threshold values are given for individual *LiteBIRD* MFT frequency bands and one band for each of LFT and HFT (quoted with their band centers). The total MFT threshold is set by the lowest threshold in the MFT bands.</caption><thead><tr><th>Band</th><th>σ<sub>Δh</sub>(Δr = 5.7 × 10<sup>-6</sup>)<br>[√GHz]</th><th>σ<sub>Δβ</sub>(Δr = 5.7 × 10<sup>-6</sup>)<br>[°√GHz]</th><th>σ<sub>Δζ</sub>(Δr = 5.7 × 10<sup>-6</sup>)<br>[√GHz]</th></tr></thead><tbody><tr><td>100 GHz (LFT)</td><td>≤0.0022</td><td>≤2.5</td><td>≤0.0012</td></tr><tr><td>100 GHz (MFT)</td><td>≤0.0023</td><td>≤2.0</td><td>≤0.0013</td></tr><tr><td>119 GHz (MFT)</td><td>≤0.0031</td><td>≤1.6</td><td>≤0.0011</td></tr><tr><td>140 GHz (MFT)</td><td>≤0.0021</td><td>≤0.8</td><td>≤0.0012</td></tr><tr><td>166 GHz (MFT)</td><td>≤0.0014</td><td>≤1.1</td><td>≤0.0010</td></tr><tr><td>Total (MFT)</td><td>≤0.0014</td><td>≤0.83</td><td>≤0.0010</td></tr><tr><td>195 GHz (HFT)</td><td>≤0.0013</td><td>≤0.83</td><td>≤0.0008</td></tr></tbody></table>

corrections for partial sky and E-B mixing are applied. We use the angular power spectrum obtained from the sum of the sky maps of the fiducial CMB and the residual as the first term in Eq. (38) to the likelihood function in Eq. (36). In the exact likelihood approach, we compute the posterior probability distribution of $r_{\Delta x}$ given our data. The likelihood analysis is performed in the range $2 \le l \le 200$, also considering the foreground residuals and noise from component separation, which allow us to properly weight each multipole. The bias $\Delta r$ due to $\sigma_{\Delta x}$ is finally quantified, as defined in Eq. (40).

In Table 18 we report the threshold value of $\Delta x$ for each systematic parameter to have a bias $\Delta r = 5.7 \times 10^{-6}$, when perturbing each systematic one at a time. We have checked that, for small enough values of $\sigma_{\Delta x}$, we have $\Delta r \propto \sigma_{\Delta x}^2$. In cases where two systematic effects are introduced simultaneously with uncorrelated errors, we expect the resulting systematic bias to be the sum of the biases resulting from individual systematic errors, i.e., $\Delta r(\sigma_{\Delta x}, \sigma_{\Delta y}) \approx \Delta r(\sigma_{\Delta x}) + \Delta r(\sigma_{\Delta y})$, where *x* and *y* are selected from *h*, *β*, and *ζ*. This has been checked by taking 200 realizations for each $\sigma_{\Delta x}$, causing approximately the same level of $\Delta r$, and comparing them with the cases of two systematic effects perturbed at the same time with those $\sigma_{\Delta x}$ values.

We checked that if we consider the residual maps for each band generated with errors just smaller than or equal to those shown in Table 18 and add them to the corresponding weights from component separation, the $\Delta r$ values associated with that come out less than half of the error budget of $5.7 \times 10^{-6}$. So, we can assume that properly performing component separation might lead to smaller residuals.

The net effect of $\Delta \zeta$ integrated over the bandwidth produces the same effect as a shift of the angle of the HWP $\Delta \rho$ and so contributes to the global uncertainty on this parameter. Its impact at first order can be reduced by minimizing the EB correlation of the CMB (as described in Ref. [80]). Because of the frequency dependence of $\Delta \zeta$, color effects between Galactic dust and the CMB might induce residual systematic contributions.

5.3.6. *Gain systematic effects.* The gain drift, i.e., the variation of the gain in time, is mainly caused by fluctuations in the focal-plane temperature, which is kept around 100 mK. Changes in this temperature cause variation of the operational point of the TES placed on the focal plane and its responsivity. We model the gain for each detector, *i*, as a time-varying gain function,
<figure><img src="image_43.png" /><figcaption>Fig. 41. Left: Systematic bias in the tensor-to-scalar ratio, r, equivalent to a bias in the B-mode power spectrum between $\ell = 2$ and 200 (the range of the recombination peak), coming from the gain drift $A_{\text{bath}}$ and the gain miscalibration $\sigma_g$ introduced in Sect. 5.3.6 for the LiteBIRD LFT 100-GHz channel. The quantity $A_{\text{bath}}$ is the amplitude of thermal bath fluctuations of a focal-plane wafer at 1 Hz, while $\sigma_g$ is the relative detector gain miscalibration uncertainty. Right: B-mode power spectrum showing the bias from unmitigated gain systematic effects in each frequency channel of LiteBIRD. The $\Lambda$CDM lensing signal is shown as the dot-dashed black line and the tensor signals for $r = 10^{-3}$ and $r = 10^{-5}$ as the dotted and dashed black lines, respectively. The blue band shows the range of power of the gain systematics biases for each individual frequency channel. The separate spectra for each channel are shown as faint blue lines. The thick red line is the systematic power spectrum from the noise-weighted average over each frequency channel.</figcaption></figure>

$G_i(t)$, which is injected in a signal timestream as

$$d_i(t) = G_i(t) [I(t) + Q(t) \cos(4\rho - 2\psi) + U(t) \sin(4\rho - 2\psi)], \quad (53)$$

where $\rho$ is the HWP rotation angle and $\psi$ is the detector polarization orientation. We model the gain function $G_i(t)$ as a product of a term $g_i^0$ that is constant in time and varying from detector to detector, and another term $g_i(T_{\text{bath}}(t))$ that depends on the thermal bath temperature $T_{\text{bath}}(t)$, varying in time, i.e., $G_i(t) = g_i^0 g_i(T_{\text{bath}}(t))$. In a simulation we generate some time variation of $T_{\text{bath}}(t)$ using a model of the power spectrum $P_{\text{bath}}$ in the frequency domain:

$$P_{\text{bath}} = A_{\text{bath}}^2 \frac{1 \text{ Hz}}{f}, \quad (54)$$

where $A_{\text{bath}}$ is the amplitude of $P_{\text{bath}}$ at a frequency $f$ of 1 Hz. The model is a good match to measurements of the temperature power spectrum using a dilution refrigerator. We assume the bath temperature $T_{\text{bath}}(t)$ to be common for all detectors on the same wafer. We assess the systematic effects by simulating a realistic LiteBIRD detector response in order to convert the focal-plane bath temperature to the responsivity of the detectors. We assume that the detector constant gain $g_i^0$ varies as a Gaussian distribution with a mean value of unity and a standard deviation of $\sigma_g$. The time-varying term $g_i(T_{\text{bath}}(t))$ starts at a value of unity at $t = 0$. Finally, we simulate using the TOAST framework the time samples encoding both astrophysical signals and thermal bath fluctuations for the full LFT focal-plane frequency channels observing for one year with LiteBIRD's scanning strategy.

In Fig. 41 (left) we show the bias introduced in the B-mode power spectrum between $\ell = 2$ and 200 in terms of the equivalent tensor mode signal and the respective value of $r$ for a simulation of the LiteBIRD 100-GHz channel with 144 detectors. We note that thermal bath fluctuations below $A_{\text{bath}} = 1 \mu\text{K}$ and relative gain miscalibration of $\sigma_g = 10\%$ (which is consistent with hardware and calibration requirements) produce a bias of $r < 10^{-5}$. Under the assumption that
the gain systematic effects are mostly uncorrelated between frequency channels of LiteBIRD, we can take this result for one of the relevant CMB frequency channels as a conservative upper limit on the systematic bias on *r* from gain systematic effects.

Figure 41 (right panel) shows the bias from unmitigated gain systematic effects in the B-mode power spectrum. The chosen input parameters for the simulations are $A_{\text{bath}} = 1~\mu\text{K}$ and $\sigma_g = 10\%$. The bias in B-mode power for each frequency channel mostly resembles a white spectrum and lies within the blue band, while the bias from the noise-weighted average over all frequency bands is shown in red. Over the relevant multipole range this bias is well below the requirement of $r < 10^{-5}$.

To further lower the residuals several detection and mitigation techniques are now being implemented and tested. They rely on identifying the leakage residuals, estimating these using approximations on the underlying signal, and then mitigating the leakage by subtraction from the data. These mitigation techniques have been shown to reduce the residuals by up to two orders of magnitude [307]. This is mostly due to the benefit resulting from high-sensitivity data and enough cross-linking redundancy provided by the LiteBIRD scanning strategy. The other sub-dominant potential sources of gain variation in time include fluctuations of the magnetic field, loading power, bias current, and gains of the cold and warm readout electronics. Given the power spectrum models for these fluctuations, we are in principle able to estimate the systematic effects in a similar manner, but this is left for future work.

The absolute gain is a conversion factor from the recorded values, usually in volts, to the physical value in units of thermodynamic temperature $K_{\text{CMB}}$. The absolute gain can be calibrated using the solar dipole signal during scans of the sky.<fnref n="26" /> We estimate the expected accuracy of the absolute gain measurement assuming the existence of 1/f noise with a knee frequency of 20 mHz for the channels with the lowest white noise for one year. We obtain a value of $5.8 \times 10^{-5}$. We employ the required value of the gain measurement from Ref. [177], on the order of $10^{-4}$, to derive a systematic bias of $5.7 \times 10^{-6}$ for the most stringent case in the highest-frequency channels. Therefore the expected bias due to the absolute gain systematic effects is $5.7 \times 10^{-6} \times (5.8 \times 10^{-5}/10^{-4})^2 = 1.9 \times 10^{-6}$.

### 5.3.7. Other systematic effects.

#### 5.3.7.1 Polarization angle

The error in the polarization angle causes mixing of the Stokes parameters between *Q* and *U*, resulting in leakage from *E*-modes to *B*-modes. Since the power of the *E*-modes is significantly larger than that of the *B*-modes, even a small amount of mixing could cause significant contamination to the *B*-modes. The leakage can be described as $C_{\ell}^{EE} \times K_{\text{pol}}$, where $C_{\ell}^{EE}$ is the *E*-mode power spectrum and $K_{\text{pol}}$ is a factor describing the polarization angle homogeneous offset in the entire sky region. The effect is sourced by the global offset of the absolute polarization angle determined by the roll angle of the spacecraft attitude or the HWP rotation position. This effect is modeled as $K_{\text{pol}} = \sin^2(2\theta_g)$, where $\theta_g$ is the global offset. We assume that the expected calibration accuracy of $\theta_g$ will be 2.7 arcmin, as described in Ref. [252] and in Sect. 4.8, yielding a $\Delta r$ value of $9.1 \times 10^{-6}$.

<footnote n="26">²⁶The orbital dipole can also be used for the absolute gain calibration. We have not yet conducted a study with the existence of possible 1/f noise, which remains for future work.</footnote>
The systematic effect of relative polarization angle uncertainties between frequency channels is studied in Ref. [250] and described in Sect. 4.8. The requirement to give the systematic bias $\Delta r = 5.7 \times 10^{-6}$ is shown in Fig. 35. The expected precision of measuring the relative angle with the Crab Nebula is shown in Table 10. Since the expectation is comparable to the requirement, we assign a $\Delta r$ value of $5.7 \times 10^{-6}$ for the systematic uncertainty of the relative angle. The frequency dependence of the polarization angle in the observation band includes the parameter $\zeta$ for the HWP (described in Sect. 5.3.5) and the sinuous-antenna wobbling effect. The net effect of those effects with an integration within the bandwidth is calibrated by the global and relative angles as described above.

The wobble effect in individual detectors is expected to be smaller than the systematic uncertainties of the HWP $\zeta$ if we assume that the uncertainty of the wobble effect is uncorrelated among the detectors, resulting in a reduction of the bias $\Delta r$ proportional to the inverse of the number of detectors in the band. We also note that the wobble effect could be canceled further if we use four sinuous-antenna patterns [272], mutually rotated by $45^\circ$ and inverted. The current focal-plane design employs this technique.

The HWP rotation position is determined by the encoder described in Sect. 4.2.4, with a demonstrated accuracy of less than 1 arcmin, yielding an expected systematic bias of $\Delta r = 1.0 \times 10^{-6}$. The time variation of the polarization angle determination accuracy of the star trackers gives a negligible systematic effect $\Delta r < 10^{-7}$ if we assume that the variation is Gaussian, since the effect cancels out for longer observation times. This may not be true when there is a long time correlation in the error; however, this systematic effect is still found to be much smaller than the constant offset effects shown above. The identification of possible error correlations needs further study.

### 5.3.7.2 Polarization efficiency

The uncertainty of the polarization efficiency is modeled as an uncertainty of the couplings of the Stokes parameters of $Q$ or $U$, described in the Mueller matrix as $M_{QQ}$ and $M_{UU}$. The polarization efficiency uncertainties are sourced by the HWP and detectors. When the polarization efficiency errors of the detectors are uncorrelated, the effect may be scaled proportional to the inverse of the square root of the number of detectors, yielding effects of negligible amplitude. On the other hand, an uncertainty in the HWP gives a significant impact, since it is common to all the focal-plane detectors. The polarization efficiency is related to the absolute normalization of the $E$- and $B$-mode power spectra. The net effect of including the efficiency from the HWP with the frequency band-average is calibrated using the $E$-mode power spectrum for higher multipoles, while its frequency dependence in the band will be calibrated using the ground facility. The expected calibration accuracy is found to be 0.2% using the $E$-mode power spectrum for a higher-multipole region, given the detector noise. If we assume that the $B$-mode lensing effect converts to the tensor signal by the amount of the uncertainty, then we expect a systematic bias of $\Delta r = 5.6 \times 10^{-6}$. We note that the uncertainties in the frequency dependence of the polarization efficiency in the observation band are taken into account through the uncertainties in the parameter $\beta$ in Sect. 5.3.5.

### 5.3.7.3 Pointing

Pointing errors can arise for several reasons: measurement uncertainties in the start tracker; misalignment between the star tracker and the boresight direction of the telescopes; vibration caused by the refrigerators; and deformation of the optical system, which includes the mirrors,
lenses, and the support structure. The frequency and direction of the pointing deviation depend on the source of the error. Models of the possible effects on the pointing will be explored in future. Here we consider two simple cases: a static bias in the pointing and Gaussian random perturbations at a sampling rate of 19 Hz. The net effect causes polarization leakage from E- to B-modes, resulting in a residual power spectrum similar to that of gravitational lensing. This effect is also independent of the existence of the PMU. We perform a simulation study for the two cases. For the offset case, we find that a pointing bias in the direction orthogonal to the scanning orientation yields the most stringent requirement of 4.6 arcmin to give $\Delta r = 5.7 \times 10^{-6}$. On the other hand, the random disturbances expected from the star-tracker pointing accuracy of 0.23 arcmin yield a negligible bias of $\Delta r < 10^{-6}$.

The wedge shape of the spinning HWP also causes pointing disturbances, where we define the wedge shape as a lack of parallelism between the plate surfaces. We therefore also conduct a simulation study assuming that the wedge produces a pointing disturbance rotating with the HWP spinning rate around the original detector pointing. This gives a requirement on the allowable maximal wedge angle in the HWP wafer fabrication with the relation $\phi = (n - 1)\psi$, where $\psi$ is the HWP wedge angle, $n$ is the refractive index, and $\phi$ is the pointing disturbance angle. For sapphire, having $n = 3.1$ as an example, we set the requirement of $\psi$ to be smaller than 4.0 arcmin to give $\Delta r = 5.7 \times 10^{-6}$.

#### 5.3.7.4 Bandpass

We study the impact of uncertainties on the frequency bandpass determination producing inter-frequency mismatch on the measurement of $r$. These uncertainties propagate into uncertainties in the amplitude of components, including leakage, after applying foreground separation. The amplitude offsets are supposed to be an average effect over the full arrays of detectors; therefore the requirement on individual detectors is relaxed by a factor of $\sqrt{N}$, where $N$ is the number of detectors in the array, with an assumption that the offset measurements are uncorrelated. This assumption is not appropriate when the correlation coefficient is larger than $1/N$, however. We describe the bandpass uncertainty using a single parameter $\gamma_c$ defined in Eq. (22), where subscript $c$ is either the dust or synchrotron component. We employ the discussions in Ref. [177] to give the systematic bias $\Delta r$ for the bandpass uncertainty with the usage of the PMU. The systematic effects are evaluated as a function of the uncertainties of $\gamma_c$ for Lite-BIRD's 15 frequency channels individually with the application of the foreground subtraction procedure [177]. With a requirement of $\Delta r$ to be less than $5.7 \times 10^{-6}$, the requirement on the measurement accuracy of $\gamma_c$ is in the range of $10^{-4}-2.5 \times 10^{-3}$, depending on the frequency channels for the detector arrays. The most stringent requirement is given for the two highest-frequency channels. Here we assume that we measure the bandpass with the most stringent requirement for all the frequency channels and calculate the systematic residual B-mode power spectrum to obtain a systematic bias of $\Delta r = 5.3 \times 10^{-6}$.

#### 5.3.7.5 Crosstalk

The readout system makes use of frequency-domain multiplexing, which is described in Sect. 4.5.2. This scheme can reduce the total number of readout wires and the heat load to the cryogenic detectors, but could also introduce crosstalk effects due to the interference of the frequency comb. Such effects are modeled as a single matrix $W_{ij}$ that describes the leakage of the measured power from detector $j$ to detector $i$. We study the crosstalk effects with a model similar to the SPT-3G design [308]. We explore 10 models of detector readout orders in the
frequency domain and detector physical positions in the focal plane to give $W_{ij}$, and find that the systematic effect is almost independent of the models, but does depend on the magnitude of the uncertainties in $W_{ij}$. For the case of 0.1% uncertainties, the systematic error on the B-mode power (lensing) is 0.07%, corresponding to $\Delta r$ of less than $10^{-7}$ in the case of the true $r$ value being 0 without foregrounds. We assign a systematic error budget of $\Delta r = 5.7 \times 10^{-6}$ to this, requiring the crosstalk knowledge uncertainty to be less than 1%.

### 5.3.7.6 Detector time constant

The TES used for LiteBIRD is known to have a time constant $\tau$ of about 3 ms, which can be well modeled by a single exponential function. The net effect of the convolution in time due to this function is a rotation angle shift $\phi$ of the HWP for the 4f-modulated signals: $\tan(4\Delta\phi) = 4\omega_{HWP}\tau$, where $\omega_{HWP}$ is the angular speed of the spinning HWP. This shift remains after applying polarization angle corrections using the $C_l^{EB}$ power spectrum and has to be further corrected. The uncertainty of $\tau$, $\delta\tau$, causes a systematic effect, which is given as $\delta\tau \simeq \delta\phi/\omega_{HWP}$. For an HWP spinning rate of 1 Hz with the expectation of $\delta\phi < 1$ arcmin (corresponding to $\Delta r = 1.0 \times 10^{-6}$) we obtain $\delta\tau < 47$ µs. A measurement of $\tau$ may be conducted in flight using the 2f-modulated signals caused by the CMB monopole $I \rightarrow P$ leakage through the HWP. The LFT HWP is expected to give a 2f-modulated signal with an amplitude of $p_0 = 9 \simeq 27$ mK. The time duration $T$ required to measure $\tau$ with a precision of $\delta\tau$ is given by

$$T = 5.7 \times 10^2 \text{s} \left( \frac{\text{NET}}{50 \mu\text{K}\sqrt{\text{s}}} \right)^2 \left( \frac{10 \text{ mK}}{p_0} \right)^2 \left( \frac{1 \text{ Hz}}{\nu_{\text{HWP}}} \right)^2 \left( \frac{47 \text{ µs}}{\delta\tau} \right)^2, \quad (55)$$

where $\nu_{HWP} = \omega_{HWP}/2\pi$. The actual time variation in the time constant is not known yet and requires study in the future. We assign $\Delta r$ a value of $5.7 \times 10^{-6}$ as the requirement.

## 5.4. Total uncertainties of the tensor-to-scalar ratio

This section provides an evaluation of the total uncertainty on the tensor-to-scalar ratio $r$. The power spectrum of the systematic effects in Eq. (38) is given as

$$C_{\ell}^{\text{sys}} = \sum_{i} C_{\ell}^{\text{sys}(i)} + \sum_{i \neq j} \langle B_{\ell m}^{\text{sys}(i)} B_{\ell m}^{\text{sys}(j)*} \rangle, \quad (56)$$

where the first term describes the sum of the $i$th systematic effect power spectrum and the second term shows the potential correlations between two systematic effects. The factors $B_{\ell m}^{\text{sys}(i)}$ and $B_{\ell m}^{\text{sys}(j)}$ are the B-mode coefficients in the spin-harmonic expansion of the $i$th and $j$th systematic residual sky maps, where the residual map is the difference between the sky maps with and without the systematic effects. Therefore, the first term is given as

$$C_{\ell}^{\text{sys}(i)} = \langle B_{\ell m}^{\text{sys}(i)} B_{\ell m}^{\text{sys}(j)*} \rangle. \quad (57)$$

With the condition of Eq. (40), we may obtain an approximation of small systematic biases:

$$\Delta r \simeq \frac{\sum_{\ell} (2\ell + 1) C_{\ell}^{\text{tens}} C_{\ell}^{\text{sys}} / (C_{\ell}^{\text{n}})^2}{\sum_{\ell} (2\ell + 1) (C_{\ell}^{\text{tens}} / C_{\ell}^{\text{n}})^2}, \quad (58)$$

where $C_{\ell}^{\text{tens}}$ is the tensor mode with $r=1$ and $C_{\ell}^{\text{n}} = C_{\ell}^{\text{lens}} + N_{\ell}$ is the sum of the lensing and the noise contributions. The first term in Eq. (56) gives the total bias of $r$ as the sum of the individual systematic biases, i.e., $\Delta r = \sum_{i} \Delta r^{(i)}$, where $\Delta r^{(i)}$ is obtained from Eq. (58) when $C_{\ell}^{\text{sys}}$ is replaced by $C_{\ell}^{\text{sys}(i)}$. The second term is estimated at the map level by summing up the residual sky maps for individual systematic sources. We evaluate the second term for the systematic effects including the beam far sidelobes, the HWP, and the gain. We compare two $\Delta r$ values: one is the
<figure><img src="image_44.png" /><figcaption>Fig. 42. Angular power spectra of the sources of systematic effects and cosmological predictions for the lensing and tensor modes, where $D_{\ell} = \ell(\ell + 1)C_{\ell}/2\pi$. From the top of the legend downwards, we plot: the lensing B-modes; the tensor B-modes with $r = 0.001$; the systematic effects of the gain; the HWP instrumental polarization (IP); the HWP parameters of $h$, $\beta$, and $\zeta$; the far sidelobes (FSL) in the angle ranges of $5^{\circ}-10^{\circ}$, $10^{\circ}-15^{\circ}$, and $15^{\circ}-180^{\circ}$; the polarization angle (Pol. Ang.); the bandpass effects; the pointing effects; and the sum of all the systematic effects. The final line is the cosmic-ray (CR) contribution as an additional noise source.</figcaption></figure>

sum of individual systematic biases $\Delta r^{(i)}$ and the other is obtained by a single power spectrum from the superimposed residual sky map of individual systematic effects. With 20 realizations, we find that the difference is $\pm 20\%$. We therefore assume that the second term in Eq. (56) cancels. An evaluation of the correlation between all systematic effects needs a single combined simulation tool, accounting for all the effects, which is beyond the scope of the current work and remains for future study. In this paper, we assume that the total power spectrum of all systematic effects is given by the first term of Eq. (56). Figure 42 shows the power spectra of individual systematic effects, as well as their sum.

Using Eq. (56) for all sources of systematic effects, including the gain, the HWP IP, the HWP parameters ($h$, $\beta$, and $\zeta$), the beam far sidelobes, the polarization angle, the bandpass, pointing, component separation, noise contributions $N_{\ell}$ from foreground subtraction in Fig. 39, and noise arising from cosmic-ray effects in Fig. 42, we estimate the total error on $r$ to be $\delta r = 1.2 \times 10^{-3}$ and the systematic bias to be $\Delta r = 0.5 \times 10^{-3}$. Figure 43 shows the likelihood function obtained as the light blue curve. Other systematic effects including the beam near sidelobe, main lobe, ghosting, beam polarization and shape in band, absolute gain, polarization efficiency, detector time constant, and crosstalk give a bias of $0.03 \times 10^{-3}$, which is obtained by summing the $\Delta r$ values in Table 15.

We consider a possible debiasing of the systematic effects by marginalization. We introduce an additional term $\alpha M_{\ell}$ in Eq. (39):

$$C_{\ell} = rC_{\ell}^{\text{tens}} + C_{\ell}^{\text{lens}} + N_{\ell} + \alpha M_{\ell}, \qquad (59)$$
<figure><img src="image_45.png" /><figcaption>Fig. 43. Likelihoods as a function of the tensor-to-scalar ratio $r$. The light blue line is that obtained from Eq. (36). The dark blue line is obtained after debiasing using marginalization.</figcaption></figure>

where $M_\ell$ is a template power spectrum accounting for the systematic bias, and $\alpha$ is a scaling factor to marginalize. We multiply the likelihood $L(r)$ in Eq. (36) by a function $\exp(-\alpha^2/2\sigma_\alpha^2)$ and integrate over $\alpha$ to define a marginalized likelihood function $\tilde{L}(r)$:

$$ \tilde{L}(r) = \int_{\alpha_1}^{\alpha_u} d\alpha L(r) e^{-\alpha^2/2\sigma_\alpha^2}. \qquad (60) $$

The net effect of this procedure is to subtract the systematic bias and to inflate the width of the likelihood function by an increase of the covariance of the power spectrum. The larger $\sigma_\alpha$ is, the smaller the systematic bias is, since the larger $\sigma_\alpha$ gives almost no constraint in the $\alpha$ value a priori, and is equivalent to fitting $\alpha$. We employ the template spectrum $M_\ell$ as the dust spectrum, which is evaluated by the component-separation process. We justify this by the fact that the main systematic effects are caused by polarization leakage from the higher-frequency channels where the dust foreground dominates, and the shape of the resultant systematic power spectra for those effects is similar to that of the dust component. For some systematic effects, including the pointing and the gain, the spectra are similar to that of lensing, since those effects are mainly caused by leakage from E-modes to B-modes, as shown in Fig. 42. Those effects, however, give rise to power spectra much lower than the others, having a spectral shape similar to that of the dust component for the lower-multipole region. We conduct the marginalization assuming $\sigma_\alpha = 10$, $\alpha_1 = 0$, and $\alpha_u = +\infty$ in Eq. (60) and obtain $\delta r = 1.0 \times 10^{-3}$. Figure 43 shows the likelihood function obtained, $\tilde{L}(r)$, as the dark blue line.

We note that the $\delta r$ value is unchanged when we increase $\sigma_\alpha$ to values larger than 10. This is due to the fact that the large value of $\sigma_\alpha$ allows $\alpha$ to move freely, which completely marginalizes over the systematic biases. A more sophisticated method is to marginalize individual systematic effects using estimated power spectrum shapes. However, this requires more precise models of the instruments and the calibration as well as an estimation of $\sigma_\alpha$ that corresponds to uncertainties of the systematic biases. These issues remain for future work.
<figure><img src="image_46.png" /><figcaption>Fig. 44. $\chi^2$ distribution when $r_{\text{true}} = 0.01$. The dark blue dashed line is obtained using the multipole range of $11 \le \ell \le 191$ containing the recombination bump. The dark blue dot-dashed line is obtained using the multipole range of $2 \le \ell \le 10$, containing the reionization bump. The light blue curve is the combined $\chi^2$. The horizontal line at $\chi^2 = 25$ indicates $5\sigma$ significance.</figcaption></figure>

We estimate the expected significance of the tensor-to-scalar ratio assuming $r_{\text{true}} = 0.01$ as the true value. We modify the measured B-mode power spectrum in Eq. (38):

$$ \hat{C}_{\ell} = r_{\text{true}} C_{\ell}^{\text{tens}} + C_{\ell}^{\text{sys}} + C_{\ell}^{\text{lens}} + N_{\ell}, \qquad (61) $$

and calculate the likelihood function in Eq. (36). We do not apply the marginalization because the total bias $\Delta r$ is much smaller than $r_{\text{true}}$. We define

$$ \chi^2(r) = -2 \log \left( \frac{L(r)}{L(\tilde{r})} \right), \qquad (62) $$

where $\tilde{r}$ is the $r$ value giving the maximum of $L(r)$, and we check the requirement shown in Table 1, i.e., we calculate the $\chi^2(r)$ value for two multipole ranges of $\ell_{\min} = 2$ and $\ell_{\max} = 10$, and $\ell_{\min} = 11$ and $\ell_{\max} = 191$; the former range includes the reionization bump and the latter includes the recombination bump. Figure 44 shows the $\chi^2(r)$ distributions for the two cases, as well as the combined one. We find that the hypothesis of $r = 0$ is rejected at a significance more than $5\sigma$ for both cases, implying that the requirement $L_{\text{1.02}}$ in Table 1 is fulfilled. We note that the asymmetric shape of the $\chi^2(r)$ distribution for the reionization bump measurement is due to the cosmic-variance dependence on $r$. The constraint on $r$ of the recombination bump measurement is stronger than that of the reionization bump for $r_{\text{true}} = 0.01$. However, the reionization bump becomes significant for lower $r_{\text{true}}$ values due to the contribution of the lensing effect.

## 5.5. Enhanced science case

The criterion for full success of the mission has been defined conservatively and does not rely on any new data sets external to LiteBIRD. However, we can reduce $\delta r$ further using external data sets, which contribute to the “enhanced science case”.27

<footnote n="27">27 In JAXA, this is also referred to as “extra success”.</footnote>
<figure><img src="image_47.png" /><figcaption>**Fig. 45.** Marginalized posterior distributions of the tensor-to-scalar ratio parameter *r* with (blue) and without (black) a delensing factor, which could be obtained by a multi-tracer template obtained by *Planck* CIB and lensing maps reconstructed with the help of *WISE* [313]. The instrumental systematic uncertainty budget is not included here.</figcaption></figure>

In this section, we focus on reducing the statistical part of the total uncertainty, $\sigma(r = 0)$, which includes the cosmic variance of the gravitational lensing signal and the noise after the component separation. One way to reduce $\sigma(r = 0)$ is to “delens” using external data sets [85–89]. Delensing removes the lensing *B*-modes by subtraction at the map level, thus reducing the lensing *B*-mode cosmic-variance contribution described above, rather than simply characterizing its power spectrum [309,310]. Successful delensing using the internal CMB data requires a higher angular resolution than that of *LiteBIRD* [311,312], because a low-noise lensing reconstruction requires the imaging of a large number of small-scale modes.

There are several promising ways to delens *LiteBIRD* with external data sets and thus contribute to the enhanced science case. The most conservative option would be to delens using currently available data sets, e.g., using the multi-tracer template obtained by the *Planck* cosmic infrared background (CIB) map and the *Planck* lensing potential map reconstructed with the help of *WISE* [313]. See Fig. 45 for the expected improvements on the constraint in *r* obtained by applying a 43% reduction of the lensing *B*-mode power spectrum as obtained for a large fraction of the sky [313]. The resulting total uncertainty, including systematic effects, is $\delta r \simeq 0.9 \times 10^{-3}$, which enables us to distinguish no primordial gravitational waves from the Starobinsky model [14] with a significance greater than $5\sigma$. *LiteBIRD* can also help reduce the Galactic dust contamination of the CIB on large angular scales, which may further improve the constraint on *r*. Further reduction of the lensing *B*-mode power spectrum can be achieved with future lensing measurements from ground-based CMB surveys, such as the CMB-S4 experiment [59] and its precursors, but only in the sky region that overlaps with *LiteBIRD*.

The other way to reduce $\sigma(r = 0)$ is to use external data sets in frequency bands below the lower edge of *LiteBIRD*'s lowest-frequency band (34 GHz). Better modeling of the polarized synchrotron emission helps to reduce the foreground residuals. External low-frequency ground-based data sets such as QUIJOTE [314,315], C-BASS [316], and S-PASS [282] at frequency
bands outside those of LiteBIRD ($\nu < 34$ GHz) would be useful for potentially improving foreground cleaning and hence reducing $\sigma(r = 0)$.

## 5.6. Summary of the statistical and systematic uncertainties in the tensor-to-scalar ratio measurement

Let us now summarize the results obtained in Sect. 5. Based on the sky observation strategy with the detector sensitivity shown in Sect. 5.1, we conducted studies on the foreground component separation in Sect. 5.2 and the systematic uncertainties in Sect. 5.3.

As explained in Sect. 5.2.3, we estimate the statistical foreground residual bias $r_{FG}$ and its error $\sigma(r = 0)$ is obtained as

$$r_{FG} = (3.3 \pm 6.2) \times 10^{-4}, \qquad (63)$$

where the definition of $\sigma(r = 0)$ is given in Eq. (35), which is equivalent to the total statistical error $\sigma_{stat}$ in the science requirement Lv2.01 in Sect. 3.2. The obtained error satisfies the requirement.

In Sect. 5.3, we give details of the studies on the systematic effects which are summarized in Table 15. The power spectra of the systematic effects are presented in Fig. 42. The bias on $r$ is defined in Eq. (40). The total systematic bias is estimated to be

$$\Delta r_{\text{syst}} = 1.7 \times 10^{-4}, \qquad (64)$$

which is equivalent to the $\sigma_{\text{syst}}$ in the science requirement Lv2.02 in Sect. 3.2.

In Sect. 5.4, we estimate the total error $\delta r$ defined in Eq. (41) by accounting for the statistical uncertainties with the component separation and the biases produced by the uncertainties of the foregrounds and instrumental systematic effects with the sky observation fraction of $f_{\text{sky}} = 49.5\%$. We apply the marginalization defined in Eq. (60) to conduct a debiasing procedure and obtain the total error

$$\delta r = 1.0 \times 10^{-3}, \qquad (65)$$

which satisfies the requirement of Lv1.01.

We also estimate the expected significance of $r$ with the assumption that the true value of $r$ is 0.01. We define $\chi^2(r)$ in Eq. (62) and compute it for two cases: the multipole range of $\ell_{\min} = 2$ and $\ell_{\max} = 10$ covering the reionization bump, and the range of $\ell_{\min} = 11$ and $\ell_{\max} = 191$ covering the recombination bump. We obtain Fig. 44 and find that both cases reject the hypothesis of $r = 0$ with a significance more than $5\sigma$, fulfilling the requirement of Lv1.02.

We point out that we could further reduce $\delta r$ using external data sets as described in the enhanced science case shown in Sect. 5.5. We examine one example by applying delensing and obtain

$$\delta r = 0.9 \times 10^{-3}, \qquad (66)$$

with an assumption of 43% reduction of the lensing effect.

# 6. Scientific outcomes of LiteBIRD beyond primordial gravitational waves

The primary science goal of the LiteBIRD mission is to discover and characterize the signature of the primordial gravitational waves from cosmic inflation in the B-mode polarization of the CMB. We have described this goal in detail in Sect. 2 and presented our forecast in Sect. 5.4.

In addition to the primary goal, the full-sky maps in 15 microwave bands will offer rich new data sets, which will enable exciting breakthroughs in a variety of science areas. We now de-
scribe the expected outcomes of the LiteBIRD mission for some representative science topics, including the reionization of the Universe (Sect. 6.1); cosmic birefringence (Sect. 6.2); the hot gas in the Universe probed using the Sunyaev–Zeldovich (SZ) effect (Sect. 6.3); spatially varying deviations from a perfect Planckian blackbody CMB spectrum (Sect. 6.4); primordial magnetic fields (Sect. 6.5); tests using polarization of the so-called “anomalies” in the temperature data (Sect. 6.6); and Galactic astrophysics (Sect. 6.7).

## 6.1. Optical depth, reionization of the Universe, and neutrino masses

The hydrogen atoms in the intergalactic medium are fully ionized in the recent Universe ($z \lesssim 6$). We have multiple lines of evidence for this from the lack of saturated hydrogen absorption lines in the spectra of quasars and afterglows of γ-ray bursts [317–321]. Given that the Universe became almost completely neutral after the epoch of hydrogen recombination (at $z \simeq 1100$), the Universe must have “reionized” during some intermediate epoch. However, we still have no clear description of the history of this epoch of reionization (EoR). In the current picture, early galaxies reionized hydrogen atoms progressively throughout the entire history between $z \simeq 12$ and $z \simeq 6$, while quasars took over reionizing helium atoms from $z \simeq 6$ to $2$ [322].

There are several ways to observe the EoR: the number counts of star-forming galaxies and quasars at $z > 6$, which were presumably producing ionizing photons [323]; the shape of the luminosity functions of high-redshift galaxies in the ultraviolet (UV) bands [324]; redshifted 21-cm lines from hydrogen atoms before the completion of reionization [325]; Doppler shifts of CMB photons by the bulk motion of ionized gas [326] (called the kinetic Sunyaev–Zeldovich effect [327]); and finally, the polarization of the CMB produced by electrons scattering quadrupole temperature anisotropies in a reionized Universe [328], which will be probed by LiteBIRD.

Thomson scattering between the CMB photons and free electrons generates linear polarization from the quadrupole moment of the CMB radiation field at the scattering epoch. This occurs at recombination as well as during the EoR. Rescattering of the CMB photons at reionization generates an additional polarization anisotropy at large angular scales, because the horizon size at this epoch subtends a much larger angular size [36]. The CMB is affected by the total column density of free electrons along each line of sight, parametrized by its Thomson-scattering optical depth, $\tau$. The wavenumber of the fluctuations contributing to quadrupole temperature anisotropy, as seen by an electron at a redshift $z$, is given by $k \simeq 3l/r_L - r(z)$ where $r_L = 14$ Gpc is the comoving distance from Earth to our last-scattering surface and $r(z)$ is the comoving distance to redshift $z$. For example, a redshift of $z = 7.7$ gives $r(7.7) = 9.1$ Gpc. We observe this wavenumber at a multipole of $l \lesssim kr(7.7) \simeq 6$, which corresponds to the so-called “reionization bump” in the polarization power spectra. The effect on the E-modes is shown in Fig. 46.

The amplitude of the reionization bump is proportional to $\tau^2 A_s$ where $A_s$ is the amplitude of the scalar curvature power spectrum. Because scattering washes out small-scale power, which is proportional to $\exp(-2\tau)A_s$, increasing $\tau$ (while fixing $A_s$) enhances the reionization bump and suppresses the small-scale power. However, since the small-scale power has been measured precisely, the value of $\exp(-2\tau)A_s$ is fixed; thus, for a given measured value of the high-$l$ power spectrum, the amplitude of the reionization bump scales as $\tau^2 A_s \propto \tau^2 \exp(2\tau)$. We can use this to determine the value of $\tau$, which in turn provides an integrated constraint on the reionization history of the Universe because $\tau = \sigma_T N_e$, where the column density of electrons is
<figure><img src="image_48.png" /><figcaption>**Fig. 46.** Planck E-mode power spectrum [329] with optical depths of $\tau = 0.055$ and $\tau = 0.000$ for constant $\exp(-2\tau)A_s$. The cosmological parameters are taken from the best-fitting ΛCDM Planck 2018 model. The green band shows the expected LiteBIRD uncertainties for each multipole in the range $\ell = 2–200$ (the cosmic-variance-limited uncertainties with $f_{\text{sky}} = 0.7$), demonstrating significant improvements over the Planck data.</figcaption></figure>

given by $N_e = c \int dt n_e$ integrated from today to the beginning of reionization. This number can be compared with the expected number of electrons from ionization by star-forming galaxies and quasars (including X-ray emission from accretion disks around black holes); see Fig. 47.

Accurate measurements of $\tau$ through the CMB are challenging because of the foreground contamination and instrumental systematic uncertainties [330]. These are most problematic on large angular scales, where the bulk of the information on $\tau$ is constrained. This difficulty is illustrated by the evolution of the constraint over time from the first WMAP release in 2003 ($\tau = 0.17 \pm 0.06$ [331]) using the TE cross-correlation, up to the latest Planck collaboration results ($\tau = 0.051 \pm 0.006$ [329], obtained from the “NPIPE” reprocessing of the Planck legacy data) using polarized EE measurements (see Fig. 48 for a complication of $\tau$ estimates, where time increases from top to bottom).

LiteBIRD will provide a cosmic-variance-limited determination of $\tau$, i.e., the sampling variance in the limit of zero instrumental noise. This is the smallest possible error bar, limited only by the fraction of sky ($f_{\text{sky}}$) available for the cosmological analysis. While ground-based experiments will provide an independent measurement of $\tau$ from gravitational lensing, with an uncertainty similar to the current Planck estimates over the coming years, a new, improved large-scale measurement will come from LiteBIRD. Assuming a fiducial value of $\tau = 0.054$ and an available sky fraction of $f_{\text{sky}} = 0.7$, the expected 68% CL uncertainty on $\tau$ is 0.002 [336], which is three times better than today's tightest bounds from Planck [329,334,337]. Not only is this a significant improvement over current measurements, but it will also be the definitive and most accurate measurement of $\tau$ from the CMB.

Beyond a cosmic-variance-limited measurement of $\tau$, the E-mode measurements carried out by LiteBIRD will constrain the precise reionization history [338]. In particular, the “dip” in the E-mode power spectrum at $\ell \simeq 10–20$ in Fig. 46 can distinguish between instantaneous reionization at a redshift of $z_{\text{reion}}$ and more physical models of reionization [339,340], includ-
<figure><img src="image_49.png" /><figcaption>**Fig. 47.** Optical depths predicted from various models of the number counts of star-forming galaxies as a function of the maximum redshift $z$. This figure is adapted from Ref. [323], showing $\tau$ inferences from different star-formation rates: the best estimates from 2013 (blue); the updated 2015 maximum-likelihood (ML) constraints (red); a model forced to reproduce the 9-yr WMAP $\tau$ constraints (orange); and a model with a truncated star-formation rate (light blue). The green band shows the expected LiteBIRD 68% and 95% confidence-level (CL) constraints on $\tau = 0.054$. The other bands show the WMAP and Planck constraints.</figcaption></figure>

ing those with a reionization history extending to longer durations and earlier onsets at $z > z_{\text{reion}}$. Although the *Planck* measurements have reconciled the value of $\tau$ from CMB polar-ization on large angular scales with the reionization process fueled by star-forming galaxies [289,323,328,341], there is still considerable uncertainty on its history [342,343] and degenera-cies with other cosmological parameters [344]. The *LiteBIRD* data will provide a significant improvement in constraints on the reionization history models, breaking the remaining de-generacies [345,346] and taking the joint CMB constraints with complementary astrophysical probes [341] to a higher level of precision. For example, in the case of homogeneous reioniza-tion models, *LiteBIRD* will reduce the uncertainty in the duration of reionization, $\Delta z_{\text{reion}}$, for an asymmetric reionization history by approximately 35% [345] with respect to the most recent *Planck* measurements [344]. The *LiteBIRD* data will also allow for testing the inhomogeneity of the reionization process and will contribute to characterizing the patchiness of reioniza-tion [347].

Beyond reionization, $\tau$ impacts the important cosmological and particle physics science topic of determining the sum of the neutrino masses, $\sum m_\nu$ [48–52,59]. Such a measurement would allow us to establish the absolute scale of the neutrino masses, and also possibly to distinguish the inverted neutrino mass ordering (i.e., two heavy, one light) from the normal ordering (one heavy, two light) [348–352]. Massive neutrinos slow down structure formation [353] and con-
<figure><img src="image_50.png" /><figcaption>**Fig. 48.** Measurements of $\tau$ based on the CMB power spectra: WMAP (2012) [42]; Planck (2013) [332]; Planck (2015) [289]; Planck (2016) [328]; WMAP+LFI (2018) [333]; Planck (2018) [279]; Planck-HFI (2019) [334]; WMAP+LFI (2020) [335]; and Planck (2020) [329]. The gray band shows the LiteBIRD uncertainty forecast for $\tau$.</figcaption></figure>

sequently we can measure the neutrino mass by comparing the amplitude of fluctuations in the low-redshift Universe with that at the last-scattering surface (i.e., $A_s$). However, we cannot determine $A_s$ unless we know $\tau$. This means that cosmological searches for the neutrino masses can be limited by the inability to measure $\tau$ with sufficient precision. For example, exact knowledge of $\tau$ would reduce the uncertainty in the neutrino mass estimates from galaxy redshift surveys by more than a factor of 2 relative to that with Planck's $\tau$ measurement (see the appendix of Ref. [52]). Thus, an improved $\tau$ determination from LiteBIRD will play a major role in measuring the neutrino mass.

When combined with measurements of the amplitude of density fluctuations at low redshifts, such as those coming from the CMB lensing data and galaxy survey observations of large-scale structure (see below), and possibly with constraints on the expansion history, a cosmic-variance-limited measurement of $\tau$ from LiteBIRD will enable a statistically significant detection of the neutrino mass, even for the minimum value $\sum m_\nu \simeq 60$ meV allowed by flavor-oscillation experiments.

Information on the amplitude of fluctuations in the low-redshift Universe might be provided by the CMB lensing data from observations of the small-scale anisotropies (e.g., those of the future Simons Array [354], SO [58,355], and CMB-S4 [356,357] experiments), and/or by data from large-scale structure surveys tracing the matter distribution (e.g., galaxy surveys with the Dark Energy Spectroscopic Instrument (DESI) [358] and the Legacy Survey of Space and Time (LSST) [359]). Distance measurements, such as those coming from BAO data, will further im-
<figure><img src="image_51.png" /><figcaption>**Fig. 49.** 2D marginalized contour levels at 68% CL for $\tau$ and the sum of the neutrino masses as measured by future combinations of CMB and large-scale structure data, including, e.g. baryonic acoustic oscillation (BAO) from DESI or galaxy lensing and clustering from LSST, adapted from Ref. [360]. The contours are centered on the fiducial values $\tau = 0.054$ and $\sum m_\nu = 60$ meV, as indicated by the cross. A cosmic-variance-limited measurement of $\tau$ is reached with LiteBIRD (i.e., $\sigma(\tau) = 0.002$). This $\tau$ limit will enable a better neutrino mass measurement, reaching a $5\sigma$ detection when combined with DESI or LSST.</figcaption></figure>

prove the constraints by adding information about the cosmic expansion history [352]. To give a specific example, combining a cosmic-variance-limited measurement of $\tau$ ($\sigma(\tau) = 0.002$) with observations of the small-scale CMB anisotropies from CMB-S4 and either BAO data from the DESI galaxy survey [358] or galaxy lensing/clustering data from the LSST survey of the Vera Rubin Observatory [359] will in both cases yield $\sigma(\sum m_\nu) = 12$ meV [360]. This will result in a detection of the neutrino mass at the $5\sigma$ level for a minimum value of 60 meV or larger. Figure 49 shows the constraining power of these data combinations in the ($\tau, \sum m_\nu$) plane, highlighting that a cosmic-variance-limited measurement of $\tau$ from LiteBIRD will be necessary for reaching a statistically significant detection of the neutrino mass from cosmological data.

To complete the picture on the neutrino sector, the expected uncertainty on the effective number of relativistic species, $N_{\text{eff}}$, from LiteBIRD alone is $\sigma(N_{\text{eff}}) \simeq 0.15$ [361], which is of the same order of magnitude as that obtained by Planck [279]. On top of that, it would give an independent measurement and an important cross-check, since any systematic effects would be different from those relevant for the high-$\ell$ Planck measurements, e.g., the modeling of small-scale foregrounds (see, e.g., Ref. [362] for a discussion of how these might affect the estimate of $N_{\text{eff}}$ provided by Planck).

The relevance of measuring $N_{\text{eff}}$ goes beyond neutrino physics, since this parameter traces the presence of any radiation-like components at the time of CMB decoupling. Such candidates include thermal axions [363,364] and sterile neutrinos [365], as well as other light species that were in thermal equilibrium with the cosmological plasma at some early point in the history of
the Universe [366]. Measuring $N_{\text{eff}}$ also allows us to probe the thermal history of the Universe with non-standard particle contents, e.g., the decay of non-relativistic massive particles [367]. A more accurate value of $N_{\text{eff}}$ would also help constrain the energy density of the stochastic gravitational wave background [368–370] because gravitational waves behave as radiation.

## 6.2. Cosmic birefringence

There exists a distinct possibility that either dark matter or dark energy (or both) is a pseudoscalar field, $\varphi$, that changes sign under parity, i.e., inversion of spatial coordinates [371,372]. If this is the case then this field could leave a unique signature in the polarization of the CMB. More generally, if the new physics that generated the initial scalar and tensor fluctuations does not violate parity symmetry and the CMB photons do not experience any parity-violating processes as they propagate for 13.8 billion years from the surface of last scattering to us, then any parity-violating correlation functions such as the temperature-B-mode (TB) correlation and the EB correlation must vanish. This is because, under spatial inversion, the spherical harmonic coefficients transform as

$$a_{\ell m}^{T} \rightarrow +(-1)^{\ell} a_{\ell m}^{T}, \quad a_{\ell m}^{E} \rightarrow +(-1)^{\ell} a_{\ell m}^{E}, \quad a_{\ell m}^{B} \rightarrow -(-1)^{\ell} a_{\ell m}^{B}. \qquad (67)$$

Consequently, the expectation value of any parity-odd observable, such as the TB and EB correlations, must vanish if the underlying physics is parity conserving. On the other hand, if the underlying physics violates parity, these correlation functions can and generically do have non-vanishing expectation values [73–77].

In this section we describe a physical effect known as “cosmic birefringence” [373–376]; see Ref. [129] for a review. The basic idea is that a new parity-violating coupling of $\varphi$ to the electromagnetic tensor makes the phase velocities of right- and left-handed polarization modes of photons different. This results in rotation of the direction of linear polarization as the CMB photons propagate through space. In other words, space filled with this $\varphi$ field behaves as if it were a birefringent medium. A homogeneous $\varphi$ field coupled to the electromagnetic field via the Chern–Simons term (i.e., the Lagrangian contains a term $-\frac{1}{4}g_{\varphi\gamma}\varphi F_{\mu\nu}\tilde{F}^{\mu\nu}$) rotates the linear polarization direction uniformly over the sky by an angle $\beta = \frac{1}{2}g_{\varphi\gamma}\int_{t_L}^{t_0} dt \, \dot{\varphi}$ [373–376], converting E-mode into B-mode polarization. The rotation angle is defined such that $\beta > 0$ corresponds to clockwise rotation on the sky in right-handed coordinates with the z-axis taken in the direction of the observer’s lines of sight. Here, $g_{\varphi\gamma}$ is the coupling constant of the interaction of $\varphi$ and photons, $F_{\mu\nu}$ and $\tilde{F}_{\mu\nu}$ are the electromagnetic tensor and its dual, and $t_L$ and $t_0$ are the time of last scattering and the present time, respectively. We would therefore observe a B-mode polarization signal even if there were no B-mode polarization initially [73]. The $\varphi$ field can be dark matter, dark energy, or both [371,372,376–378]; thus, discovery of such a signal would be a major breakthrough in cosmology and fundamental physics.

When $\beta$ does not vanish, the observed E- and B-mode spherical harmonics coefficients of the CMB, $a_{\ell m}^{E, \text{CMB}, o}$ and $a_{\ell m}^{B, \text{CMB}, o}$, are related to those at the surface of last scattering as

$$a_{\ell m}^{E, \text{CMB}, o} = a_{\ell m}^{E, \text{CMB}} \cos(2\beta) - a_{\ell m}^{B, \text{CMB}} \sin(2\beta), \qquad (68)$$

$$a_{\ell m}^{B, \text{CMB}, o} = a_{\ell m}^{E, \text{CMB}} \sin(2\beta) + a_{\ell m}^{B, \text{CMB}} \cos(2\beta). \qquad (69)$$

However, this effect is degenerate with an instrumental miscalibration of polarization angles [91,379]. This means that, in the absence of any other information, we cannot tell whether the polarization angle of the CMB is rotated by the new physics (i.e., $\beta$) or the polarization
angle of detectors is rotated with respect to the sky coordinates by a miscalibration angle $\alpha$. As a result, in the absence of any other information, we can only determine the sum of the two angles, $\alpha + \beta$, which explains why the previous determinations of $\beta$ were spread over a wide range beyond the quoted statistical uncertainties (see Refs. [78–80] for summaries).

Since the miscalibration angle $\alpha$ generates a spurious B-mode power spectrum $\sin^2(2\alpha)C_{\ell\ell}^{EE}$, we must calibrate the angles with a precision sufficient to achieve the full science requirements of LiteBIRD. As described in Sect. 4.8.4, there are no strong polarized astrophysical sources on the sky with precisely known polarization angles; thus, the calibration must rely on measurements on the ground, which currently limit the accuracy of calibration to a bit better than 1°. This accuracy is much poorer than the requirement, which is of order 0.05°, demanding substantial improvements in the ground-calibration methodology. As a result, the option of using the TB and EB correlations (assumed to vanish) to calibrate the instrumental polarization angles [255] was considered in Sect. 4.8.4. This method is based on accurate knowledge of the cosmological TE and EE power spectra from the scalar mode, and it is straightforward to fit to the TB and EB power spectra to solve for $\alpha$ with arcminute precision. However, not only does this “self-calibration” procedure eliminate LiteBIRD’s sensitivity to the new physics of uniform rotation caused by $\varphi$, but it can also bias the angle calibration if non-zero $\beta$ is present in the data.

A new method [251–253] has been developed to mitigate this issue. The basic idea is to use the polarized Galactic foreground emission as the angle calibrator. Cosmic birefringence is a cumulative effect, whose magnitude is proportional to the path length of photons. Since the origin of the Galactic emission is much closer to Earth than the CMB, we can ignore $\beta$ in the Galactic emission (labeled ‘fg’ below). On the other hand, the miscalibration angle $\alpha$ affects both the CMB and the Galactic emission. We thus write

$$a_{\ell m}^{E,o} = a_{\ell m}^{E,fg} \cos(2\alpha) - a_{\ell m}^{B,fg} \sin(2\alpha) + a_{\ell m}^{E,CMB} \cos(2\alpha + 2\beta) - a_{\ell m}^{B,CMB} \sin(2\alpha + 2\beta), \quad (70)$$

$$a_{\ell m}^{B,o} = a_{\ell m}^{E,fg} \sin(2\alpha) + a_{\ell m}^{B,fg} \cos(2\alpha) + a_{\ell m}^{E,CMB} \sin(2\alpha + 2\beta) + a_{\ell m}^{B,CMB} \cos(2\alpha + 2\beta), \quad (71)$$

which gives the ensemble average of the observed EB power spectrum as [251]

$$\langle C_{\ell}^{EB,o} \rangle = \frac{\tan(4\alpha)}{2} \left( (C_{\ell}^{EE,o}) - (C_{\ell}^{BB,o}) \right) + \frac{\sin(4\beta)}{2\cos(4\alpha)} \left( \langle C_{\ell}^{EE,CMB} \rangle - \langle C_{\ell}^{BB,CMB} \rangle \right) \\ + \frac{1}{\cos(4\alpha)} \langle C_{\ell}^{EB,fg} \rangle + \frac{\cos(4\beta)}{\cos(4\alpha)} \langle C_{\ell}^{EB,CMB} \rangle. \quad (72)$$

This formula allows us to determine $\alpha$ and $\beta$ simultaneously. Note that we do not need to know the intrinsic foreground EE and BB power spectra, but only those of the CMB as well as the intrinsic EB spectra of the CMB and the foreground emission. The last term on the right-hand side is the intrinsic EB power spectrum of the CMB at the surface of last scattering. This can be produced by parity-violating gravitational waves from, e.g., the gravitational Chern–Simons term [73,76,77] and Abelian [136] or non-Abelian [137] gauge fields during inflation (see Sect. 2.5). Because these effects arise from gravitational waves, the shape of the EB power spectrum is completely different from that generated by cosmic birefringence; thus, we can simultaneously fit $\beta$ and $C_{\ell}^{EB,CMB}$ without losing sensitivity [137,380].

The intrinsic EB correlation of the foreground emission, $C_{\ell}^{EB,fg}$, requires careful treatment. The worst-case scenario is that $C_{\ell}^{EE,fg}$ is proportional to $C_{\ell}^{EE,fg} - C_{\ell}^{BB,fg}$. In this case, fitting the observed EB power spectrum gives $\beta - \gamma$ rather than $\beta$, where $\gamma$ is an additional effective angle
<figure><img src="image_52.png" /><figcaption>**Fig. 50.** The *B*-mode power spectrum from cosmic birefringence with $\beta = 0.35^{\circ} \pm 0.14^{\circ}$ (the blue shaded area). The dotted and dot-dashed lines show the *B*-mode power spectra from primordial gravitational waves with $r = 0.00461$ (Starobinsky's model [14]) and 0.001, respectively, both of which are detectable by *LiteBIRD* by design. The red line shows the *B*-mode power spectrum from the gravitational lensing effect of the CMB. The purple shaded regions show the expected constraints from *LiteBIRD*, derived in Sect. 5.2.</figcaption></figure>

from the intrinsic foreground *EB* correlation [251]. Since the *Planck* data show a positive *TE* and *TB* signal for thermal dust emission [60,78], we expect a positive *EB* as well [256], although this is yet to be found and its strength may depend on the Galactic mask [257,381].

Recently, the value $\beta - \gamma = 0.35^{\circ} \pm 0.14^{\circ}$ (68% CL) has been obtained from the *Planck* 2018 data [80]. If $\gamma > 0$, as implied by the positive *TE* and *TB* measurements of the polarized dust emission, the significance of $\beta$ would increase further. This has been confirmed by the latest analysis using the *Planck* Public Data Release 4 (PR4) “NPIPE” reprocessing [329]. Thanks to the lower noise and better-characterized systematics, an improved measurement, $\beta - \gamma = 0.30^{\circ} \pm 0.11^{\circ}$ (68% CL), is obtained for nearly full-sky data [381]. Accounting for $C_{\ell}^{EB,fg}$ using the *TE* and *TB* measurements of the polarized dust emission and a physical model of Ref. [257], $\beta = 0.36^{\circ}$ is found with the same statistical uncertainty. The impact of the known systematics of the *Planck* data on $\beta$ is found to be negligible compared to the statistical uncertainty. It is too early to tell if this tentative hint for $\beta$ is due to cosmic birefringence, some unknown systematics in the *Planck* data, some unexpected property of the intrinsic foreground *EB* correlation, or a combination of them. In any case, *LiteBIRD* can play a decisive role in understanding the origin of any such signal. With *LiteBIRD*, we can reduce the statistical uncertainty to below $0.1^{\circ}$ [253], potentially increasing the statistical significance from $3\sigma$ to the level of discovery.

If cosmic birefringence exists at the level of $\beta = 0.35^{\circ}$, it will produce *B*-mode polarization of $\sin^2(2\beta)C_{\ell}^{EE}$, which is within reach of *LiteBIRD*'s sensitivity (Fig. 50). Therefore, *LiteBIRD* can test for the presence of the cosmic birefringence signal, not only via the *EB* correlation, but also the *B*-mode power spectrum.

So far we have only discussed birefringence that is the same in every direction, but in principle we can also look for an *anisotropy* in $\beta$ [382,383]. This can arise from inhomogeneity of
<figure><img src="image_53.png" /><figcaption>**Fig. 51.** Spectrum of the thermal SZ effect. The solid line shows the example spectrum in units of kJy sr⁻¹ in the non-relativistic limit, $k_B T_e / m_e c^2 \ll 1$, with $y = \tau_h k_B T_e / m_e c^2 = 5 \times 10^{-6}$. The dashed line shows the spectrum of the relativistic correction for $\tau_h = 0.01$ and $k_B T_e = 5$ keV (which gives $y = 9.8 \times 10^{-5}$) relative to the non-relativistic SZ effect with $y = 9.5 \times 10^{-5}$, calculated by the `SZpack` package [398,399], i.e., the difference between the full and non-relativistic spectra in units of kJy sr⁻¹. We show the difference between the relativistic and non-relativistic SZ spectra with two different values of y to highlight the genuine effect of the relativistic correction, which cannot be absorbed by changing the value of y. The colored bars show the sensitivity of the 15 partially overlapping bands of the *LiteBIRD* detectors. For clarity we show half of the bands as positive and the other half as negative quantities, but only their absolute values are meaningful.</figcaption></figure>

dark matter (and possibly also dark energy) made of a pseudoscalar field coupled to the electromagnetic tensor via the Chern–Simons term. This introduces correlations between *T* and *B* and between *E* and *B* at different multipoles (i.e., $\langle a_{\ell m}^T a_{\ell' m'}^E \rangle$ and $\langle a_{\ell m}^E a_{\ell' m'}^B \rangle$) in a manner similar to the gravitational lensing effect of the CMB [384–386]. Such a property makes it possible to create a map of $\beta$ in each *LiteBIRD* pixel. This map will be useful not only for probing new parity-violating physics, but also for characterizing instrumental systematics (i.e., anisotropy in $\alpha$). So far there is no evidence for the anisotropic signal of $\beta$ [387–390], but *LiteBIRD* can certainly improve upon the constraints, particularly on the largest angular scales.

### 6.3. *Mapping the hot gas in the Universe*

Electrons in the hot ionized gas transfer their energy to CMB photons by inverse Compton scattering, leading to a characteristic distortion of the blackbody spectrum of the CMB (see Fig. 51). This phenomenon is known as the thermal SZ effect [391,392] and has been routinely detected in the directions of galaxy clusters [393–397]. The amplitude of the thermal SZ effect is characterized by the so-called “Compton y parameter”, which is given by $y = \tau_h k_B T_e / m_e c^2$ where $\tau_h$ is the optical depth of hot gas (which should be distinguished from the optical depth of reionization, $\tau$, discussed in Sect. 6.1) and $T_e$ and $m_e$ are the electron temperature and mass, respectively.
<figure><img src="image_54.png" /><figcaption>**Fig. 52.** Reconstructed power spectrum of the thermal SZ effect from a simulation of *LiteBIRD* (red line), compared with the input one (black line). Both agree well except at $\ell < 20$, which still shows residuals from the Galactic emission; however, such low multipoles suffer from large non-Gaussian cosmic-variance error bars in any case. The noise power spectrum of *LiteBIRD* (green dashed line) is much lower than that of *Planck* (gray dotted line), showing the substantially improved sensitivity and fidelity of the thermal SZ map from *LiteBIRD*.</figcaption></figure>

Using the so-called constrained internal linear combination (cILC) method [400–403], we can reconstruct an all-sky map of the thermal SZ signal and its angular power spectrum with minimal residual foreground contamination [404–406]. Applying the same component-separation algorithm that was used on the *Planck* data to the *LiteBIRD* simulations, we find that, while the *Planck* SZ map still contains contamination of various foreground sources due to the limited number of frequency bands, *LiteBIRD* can faithfully reconstruct the SZ map at $\ell > 20$. Figure 52 shows the power spectrum of the reconstructed SZ map from this simulation. We also find that, while the reconstruction noise in the *Planck* SZ map is comparable to the signal itself, *LiteBIRD* will reduce the noise by an order of magnitude relative to the *Planck* SZ map.

Exploiting the 15 *LiteBIRD* frequency bands will yield a much improved, high-fidelity SZ map over the full sky at $\ell \le 200$, essentially free of contamination. This full-sky map will show in projection all hot gas in the Universe and will have a lasting impact on astrophysics as legacy data from *LiteBIRD*. An important application of this full-sky thermal SZ map will be to cross-correlate it with 3D catalogs of galaxies with known redshifts, as demonstrated in Refs. [407–412]. This cross-correlation allows us to perform *tomography* of the hot gas in the Universe as a function of cosmic time, which can test theories of structure formation [413]. This full-sky thermal SZ map can also be used to search for the warm–hot intergalactic medium (WHIM) by stacking at the positions of known galaxy pairs [414].
The high-fidelity SZ map is also useful for studying an inhomogeneous reionization process via the cross-correlation of the SZ map and fluctuations of the CMB optical depth, $\delta\tau_h$ [415]. The SZ effect can be generated by hot electrons during the reionization epoch and the optical depth has fluctuations due to the spatial variations of the electron density during the reionization epoch. Thus, an inhomogeneous reionization process leads to a correlation between the SZ map and $\delta\tau_h$. This cross-correlation is much less affected by late-time galaxy cluster contributions than the SZ auto-power spectrum. In addition, this cross-correlation can constrain the temperature of ionized bubbles, while the kinetic SZ and $\delta\tau_h$ cannot. The cross-correlation signal would be detectable by comparing the *LiteBIRD* SZ map with a $\delta\tau_h$ map reconstructed from CMB-S4 for an interesting parameter space for the ionized bubbles, e.g., a characteristic size of 5 Mpc and a temperature of $5 \times 10^4$ K.

The solid line in Fig. 51 shows the SZ spectrum in the non-relativistic limit (where $k_B T_e / m_e c^2 \ll 1$) for an example value of $y = 5 \times 10^{-6}$. The shape is universal and depends only on the mean CMB temperature; however, small relativistic corrections to this shape exist [416–418] and are proportional to $k_B T_e / m_e c^2$ at leading order [419–421]. Detecting this relativistic correction averaged over a full-sky SZ map [422,423] can yield the mean gas temperature of the Universe, providing an “integral constraint” on physics of the intergalactic medium [424] and stringent and robust constraints on the energy feedback from supernovae and active galactic nuclei (AGN). This complements information that can be obtained about the relationship between the WHIM and halos from correlating the SZ map with gravitational lensing and other tracers of large-scale structure [425,426].

The dashed line in Fig. 51 shows the relativistic correction for $\tau_h = 0.01$ and $k_B T_e = 5$ keV, relative to the non-relativistic SZ effect with $y = 9.5 \times 10^{-5}$. Note that $\tau_h = 0.01$ and $k_B T_e = 5$ keV correspond to $y = 9.8 \times 10^{-5}$. One may wonder why we do not take the difference between the relativistic and non-relativistic SZ spectra with the same $y = 9.8 \times 10^{-5}$. This would give a larger difference at low frequencies, giving the impression that we can detect the relativistic correction without the high-frequency data; however, this is false, since the difference at low frequencies can be compensated by slightly changing $y$, as we have done here [427], whereas the distortion at high frequencies is genuine and cannot be compensated by changing $y$. Observing above 300 GHz with *LiteBIRD* will thus give a great advantage compared to lower-frequency ground-based surveys for detecting the relativistic correction, since the electron-temperature dependence of the relativistic correction manifests itself mostly at high frequencies, as shown in Fig. 51. Finally, the *LiteBIRD* sensitivity and high-frequency coverage could be used to apply the method proposed in Ref. [428] to improve constraints on the monopole of the $y$-type distortion of the CMB spectrum via the spectrum of the SZ effect.

## 6.4. Anisotropic CMB spectral distortions

Although *LiteBIRD* is not sensitive to the spatially uniform (i.e., monopole) component to the distortion of the Planckian spectrum of the CMB, it is very sensitive to any *spatially varying* component of the spectral distortion. The thermal SZ effect, described in Sect. 6.3, is one example of such a spectral distortion and can be used to map the distribution of hot gas in the Universe.

Another example of anisotropies with spectral dependence different from that of the CMB is Rayleigh scattering of the CMB photons [429–431]. The CMB decouples from electrons at a redshift of $z \simeq 1090$, leading to the usual “surface of last scattering”. Around this epoch the
electrons combined with protons to form neutral hydrogen atoms, which also scatter photons via Rayleigh scattering. The cross-section for Rayleigh scattering has a characteristic frequency dependence of $\sigma_R(\nu)\propto\nu^4$ [432], which leaves a frequency-dependent imprint in temperature and polarization anisotropies of the CMB. No monopole spectral distortion is produced by Rayleigh scattering.

The cross-section of Rayleigh scattering is given more explicitly by

$$ \sigma_R(\nu) = \sigma_T \left( \tilde{\nu}^4 + \frac{638}{243} \tilde{\nu}^6 + \dots \right), \qquad (73) $$

where $\sigma_T$ is the Thomson cross-section, $\tilde{\nu} \equiv \nu/(\sqrt{8}/9cR_\infty) \simeq \nu/(3.1 \times 10^6 \text{ GHz})$, and $R_\infty$ is the Rydberg constant. This additional scattering produces a frequency-dependent shift of the peak of the visibility function to $z < 1090$, which modifies the Silk damping process and shifts the locations of the acoustic peaks in the temperature and E-mode polarization power spectra [430,431].

The Rayleigh-scattering signal can be extracted from the data by cross-correlating the frequency-independent Thomson-scattering component and the frequency-dependent ($\nu^4$) Rayleigh-scattering component [431,433]. Since the frequency dependence of these two components is known precisely, we can separate them by explicitly rejecting the other component using the CILC method [400–403]. The remaining largest contamination is then the CIB. While the expected Rayleigh-scattering signal is small in the frequencies observed by LiteBIRD, the superb sensitivity of LiteBIRD allows for detection of the signal at a statistical significance of $25\sigma$ [433], which is due primarily to the temperature data. The first detection of this signal would not only be a significant achievement in cosmology (because this is a firm prediction of the standard model of cosmology that is yet to be confirmed), but could also improve determination of the cosmological parameters such as $N_{\text{eff}}$ and $\sum m_\nu$ [433], as well as primordial non-Gaussianity [434].

It may also be possible to probe epochs earlier than CMB last scattering using different spectral information. Let us review the well-known physics of the early Universe (see Refs. [435,436] for reviews). When the temperature of the Universe exceeded $5 \times 10^6 \text{ K}$ (or $z > 2 \times 10^6$ in terms of redshift), double-Compton scattering (which changes the total photon number) was efficient in relaxing the photon spectrum to a Planck spectrum with a vanishing chemical potential, even if some extra energy was injected into the plasma [437,438]. If the energy was injected in the range $5 \times 10^4 < z < 2 \times 10^6$, double-Compton scattering would no longer be fast enough to relax the photon spectrum to a blackbody with no chemical potential. Compton scattering is still efficient for redistributing the photon energies so as to maintain an equilibrium distribution (i.e., a Bose–Einstein distribution with non-zero chemical potential, also known as the “$\mu$ distortion”). If the spectral distortion was caused by an energy injection after $z = 5 \times 10^4$, it would not relax to an equilibrium distribution because the energy exchange due to Compton scattering would be inefficient, resulting in, e.g., a permanent spectral distortion such as the SZ effect described in Sect. 6.3.

While there exist many theoretical possibilities for energy injection in the early Universe before $z = 5 \times 10^4$ [435,436], one mechanism present in the standard model of cosmology is energy injection from the dissipation of sound waves [439,440]. Because this spectral distortion occurs at second order in the perturbation, the energy injection rate due to the dissipation of sound
waves is proportional to the sound wave amplitude squared. This property makes it possible to constrain the small-scale power of fluctuations from the chemical potential [441–443].

This phenomenon offers the third example of anisotropic spectral distortions of the CMB. While this signal is isotropic in the sky if the fluctuations obey Gaussian statistics, a specific type of non-Gaussian fluctuations, called “squeezed non-Gaussianity”, can be produced by certain physical mechanisms during inflation, such as multi-field effects and non-Bunch-Davies vacuum initial conditions, and would generate spectral distortions characterized by a *spatially varying* chemical potential [444,445]. *LiteBIRD* can look for this signal by cross-correlating the measured temperature anisotropies with a map of the chemical potential reconstructed from *LiteBIRD*'s multi-frequency data, since this cross-correlation measures a three-point function (temperature fluctuation on large scales correlated with the squared amplitude of sound waves on small scales). Although the multi-field effect of inflation yields only a small signal-to-noise ratio for *LiteBIRD*, given the constraints on this type of non-Gaussianity from the *Planck* data [64,65], non-vacuum effects can yield a large signal-to-noise ratio, offering a powerful test of the physics of inflation at its onset [445]. Additional sources of anisotropic spectral distortion arise when the background spacetime is itself anisotropic [446–448].

There may be additional spectral-spatial variations in the CMB that can be probed by *LiteBIRD*, including those that affect polarization. As an example, light axion-like particles are converted into photons in the presence of magnetic fields, generating anisotropic distortions in the CMB spectrum [449]. In particular, “resonant conversion” of axions into photons by the Galactic magnetic field yields *polarized* spectral distortions of the CMB with the spatial distribution of the signal tracking the Galactic magnetic field. *LiteBIRD* can search for signals of this type.

### 6.5. *Primordial magnetic fields*

Magnetic fields are ubiquitous in the Universe at all scales from planetary systems to clusters of galaxies, with hints of their presence also in filaments of the large-scale structure [450,451] as well as in intergalactic voids [452–454]. The dynamo is a popular explanation for the origin of the cosmic magnetic fields, but it requires an initial field. In other words, a dynamo does not generate a field from nothing, but amplifies the existing “seed” field (see Ref. [455] for a review). Even the simplest adiabatic compression of magnetic fields frozen to plasma could explain the observed strength of $\mu$G in galaxy clusters if the pre-compression strength were 0.1 nG ($= 10^{-10}$ G) in intergalactic space [456]. Lower bounds on the intergalactic magnetic fields of $\gtrsim 10^{-16}-10^{-18}$ G have been inferred from the lack of extended $\gamma$-ray halos towards blazars in the GeV energy bands (see Ref. [457] and references therein). If they exist in intergalactic space and within filaments of the large-scale structure, where do they come from (see Refs. [458–462] for reviews)?

Cosmic inflation may provide this origin, i.e., the primordial magnetic fields (PMFs). In the standard model of elementary particles and fields, massless gauge fields (such as electromagnetism) are conformally coupled to gravity; thus, no metric excitation of spin-1 fields is possible in an expanding Universe, unlike for scalar [20] and tensor [25,26] perturbations. In other words, we must break the conformal invariance of the massless gauge fields to generate the PMFs during inflation [463,464]. Although theoretical challenges remain [465,466], searching for the PMFs has profound implications for our understanding of physics of inflation beyond
the origin of the scalar and tensor perturbations that we have already discussed extensively in this paper.

The power spectrum of the stochastic PMFs, $P_B(k)$, is defined by

$$ \langle B_i(\mathbf{k}) B_j^*(\mathbf{q}) \rangle = (2\pi)^3 \delta_D^{(3)}(\mathbf{k} - \mathbf{q}) P_B(k) \left( \delta_{ij} - \frac{k_i k_j}{k^2} \right), \quad (74) $$

where $\delta_D^{(3)}(\mathbf{k})$ is the Dirac delta function and $\delta_{ij}$ is the Kronecker delta. A power-law power spectrum, $P_B(k) = A_B k^{n_B}$, is characterized by its spectral index, $n_B$, and amplitude, $A_B$. Spectral indices generated during inflation are usually negative [463], e.g., a nearly scale-invariant spectrum of $n_B = -2.9$ (the exact scale invariance $n_B = -3$ gives a diverging energy density of the magnetic field), whereas they are positive if generated by causal processes such as the electroweak phase transition [467–469], i.e., an even integer of $n_B \ge 2$ [470].

In the literature, the magnetic field strength is usually quoted as the value extrapolated to the present epoch assuming that the fields are frozen to the plasma. Moreover, we often quote the field value smoothed over a length scale $\lambda$, calculated from the smoothed energy density of the magnetic field:

$$ \rho_{B,\lambda} = \int_0^\infty \frac{d^3 k}{(2\pi)^3} P_B(k) e^{-k^2 \lambda^2} = A_B \frac{\Gamma[(n_B + 3)/2]}{4\pi^2 \lambda^{n_B + 3}}. \quad (75) $$

Here we choose $\lambda = 1$ Mpc and compute the field strength, $B_\lambda$, from the energy density $\rho_{B,\lambda} = B_\lambda^2/2$, in natural units. We then convert the units and quote the value of $B_\lambda$ in units of nG. The convenient quantity is the ratio of $\rho_{B,\lambda}$ to the CMB photon energy density, $\rho_\gamma$:

$$ \frac{\rho_{B,\lambda}}{\rho_\gamma} = 9.53 \times 10^{-8} \left( \frac{B_\lambda}{1 \text{ nG}} \right)^2, \quad (76) $$

for the present-day CMB temperature of 2.7255 K.

The sensitive CMB polarization measurements provided by LiteBIRD can constrain the PMFs with unprecedented precision, using three of the distinct ways that the PMFs affect the CMB:

(1) A stochastic background of PMFs provides the stress–energy source, which gravitationally induces scalar, vector, and tensor perturbations in the CMB and hence the TT, EE, and BB power spectra (see Refs. [456,471] for reviews). In addition, helical magnetic fields can produce parity-violating correlations such as the TB and EB power spectra [472–474]; Ref. [475] describes the constraints on helical fields from the Planck data. Because the temperature and polarization anisotropies from the magnetic-field stress–energy tensor are highly non-Gaussian, higher-order correlations are a powerful probe of the PMF [476–478].

(2) Dissipation of the fields heats gas in intergalactic space, altering the thermal history of the Universe [479]. This affects the CMB primarily via the optical depth of electron scattering (Sect. 6.1) [480] and the isotropic (monopole) distortion of the blackbody spectrum of the CMB [481,482].

(3) Faraday rotation induces frequency-dependent $(\alpha\nu^{-2})$ inhomogeneous rotation of polarization angles of the CMB over the sky [483–485], which has been constrained by Planck [390,475] and also by ground-based experiments [486,487].

LiteBIRD is sensitive to all of the above effects, except for the monopole spectral distortion of the CMB.
<figure><img src="image_55.png" /><figcaption>**Fig. 53.** B-mode power spectra from the gravitational effects of an inflationary PMF with $n_B = -2.9$ (red) and a causally generated PMF with $n_B = 2$ (blue), in comparison to those of the primordial tensor perturbation of the Starobinsky model of inflation with $r = 0.0461$ (solid black) [14] and gravitational lensing (dotted). The amplitudes of the PMF curves are chosen to match the current limits, $B(1 \text{ Mpc}) = 2 \text{ nG}$ and $0.003 \text{ nG}$ for $n_B = -2.9$ and 2, respectively. The blue shaded regions show the expected constraints from LiteBIRD, derived in Sect. 5.2. The inset shows the power spectra up to higher multipoles.</figcaption></figure>

For the first effect, the study in Ref. [488] shows how LiteBIRD's B-mode power spectrum can improve upon the constraints on the amplitude of a stochastic background of PMFs, and how in combination with future ground-based experiments it can reach the nG-level field strength for a nearly scale-invariant spectrum ($n_B = -2.9$) and the sub-pG level for causally generated fields ($n_B = 2$). However, since the gravitational contribution of PMFs to the B-mode power spectrum can be degenerate with those of primordial gravitational waves and lensing in the multipole range probed by LiteBIRD [488,489; see Fig. 53], we can use higher-order correlations (i.e., non-Gaussianity) to further constrain PMFs. Using LiteBIRD's B-mode information, the sensitivity to the magnetized tensor bispectrum is improved by more than two orders of magnitude compared to the Planck results, yielding a 1-nG-level constraint [138]. For the second effect, LiteBIRD's cosmic-variance-limited measurement of E-modes at low and intermediate multipoles (Sect. 6.1) will significantly tighten the constraints on the post-recombination heating due to PMFs [480,490,491], going beyond the nG-level constraints. For the third effect, LiteBIRD can target Faraday rotation with its unique frequency dependence, $\nu^{-2}$, which can be used to cross-check the PMF constraints from the B- and E-mode spectra. In particular, using the Faraday rotation angle power spectrum, LiteBIRD can improve upon the limit on the amplitude of a nearly scale-invariant PMF by more than an order of magnitude, reaching the nG threshold, provided that the contamination of Faraday rotation from our Galaxy is properly modeled [492].

While simultaneous analyses of all of these effects have not been performed in the literature, it is possible that a combination of all the information could push the limit down to the 0.1-nG level robustly for a power spectrum with power-law index $n_B = -2.9$. A comprehensive analysis of the forecasts for LiteBIRD including all these effects is left for future work.
6.6. *Elucidating spatial anomalies with polarization*

A six-parameter cosmological model seems to explain the observed structure of the CMB in remarkable detail [279,289,493]. Nevertheless, there exist some features in the data that exert a mild tension against such a model, potentially providing hints of new physics to be explored. Below, we summarize the most important of these spatial “anomalies”, found at modest levels of statistical significance (i.e., 2–3σ) in the *WMAP* and *Planck* temperature data [66,92–94].

(1) *Low-ℓ power deficit (low variance).* A simple statistical measure of the data is afforded by the variance of the signal as a function of angular scale. Analysis of low-resolution maps reveals a lack of variance when compared to simulations based on the best-fitting cosmological model. The map-based variance is dominated by contributions from large angular scales on the sky, whereas the cosmological parameter fits are insensitive to these modes and mostly determined by scales corresponding to $ℓ > 50$. The dearth of large-angular-scale power in the *Planck* power spectrum for $ℓ < 30$ results in an apparently anomalous variance, with a *p*-value (the probability that simulations yield a lower value than the data) of order 1%.

(2) *Lack of correlation on large angular scales.* A lack of structure in the angular two-point correlation function is observed for angular separations larger than 60°, with the observed values lying close to zero between ≃ 60° and ≃ 170°. This is captured by the use of a posteriori statistic, $S_{1/2}$, as proposed by the *WMAP* team and given by the integral of the squared correlation function for angular separations larger than 60°. A corresponding *p*-value of less than 1% was determined from the 2018 *Planck* data.

(3) *Alignments of low-multipole moments.* Detections of the alignment of the quadrupole and octupole moments of the CMB temperature distribution have been found since *WMAP*'s first release. This is somewhat surprising, given that the temperature anisotropies are expected to have random phases in the standard cosmological model, implying that the multipole moments should be uncorrelated. The analysis in Ref. [92] finds that the quadrupole and octupole orientations are aligned to within about 10°, with a *p*-value lower than 1%.

(4) *Hemispherical asymmetry.* The standard cosmological model predicts that the same power spectrum should be measured in different patches of the sky, except for variations connected with sample variance. However, *WMAP* and *Planck* data show evidence for a hemispherical asymmetry (or dipolar modulation) of power in a particular direction. For the *Planck* data in particular, several methods to test for such asymmetry have been applied and compared [93,94]. These are sensitive to either amplitude, directional-ity, or both, although they do differ in terms of their weighting of power on different scales. Nevertheless, the results are all consistent with a modulation of power of around 7% between two hemispheres defined by the preferred direction (l, b) = (209°, -15°), extending over scales to $ℓ_{max}$ ≃ 60 with a significance approaching 3σ. Interestingly, one such test, based on the anomalous clustering of directions within bands of multipoles, suggests that this asymmetry holds even on relatively small scales.

(5) *Parity asymmetry.* To test whether the CMB is symmetric with respect to reflections about the origin, $\hat{n} \rightarrow -\hat{n}$, the CMB anisotropy field can be divided into symmetric and anti-symmetric functions with even and odd parity, corresponding to spherical harmonic modes with even and odd $ℓ$ values, respectively. On the largest angular scales, the Uni-
verse should be parity neutral, yet an odd-parity preference has been established from analysis of the *WMAP* and *Planck* data sets. By computing the ratio between the sum of even and sum of odd modes up to a given $\ell_{max}$, it was found that the significance varies with the maximum multipole chosen, with *p*-values of about 1% for $\ell_{max}$ = 20–30 [94].

(6) *Cold spot*. A particularly large cold region was originally discovered in the *WMAP* first-year data from the study of the kurtosis of spherical Mexican-hat wavelet (SMHW) coefficients over a range of angular scales. It corresponds to an anomalous temperature feature in these coefficients on angular scales of ≃ 10°, with the structure centered at Galactic coordinates (l, b) = (210°, -57°). Less than 1% of simulations based on the standard ΛCDM cosmological model yield a kurtosis at least as large as that seen in the data [93], with some dependence on a posteriori choices. Novel theories have been invoked to explain the cold spot, including the gravitational effect produced by a collapsing cosmic texture.

A conservative explanation for these anomalies, given their claimed levels of significance, is that they are mild statistical excursions whose significance is overestimated due to the application of a posteriori statistics. However, whether this is the case, or, alternatively, that they reflect real physical properties of the Universe cannot be elucidated further using the temperature anisotropies, which are already cosmic-variance limited. Instead, new observations are needed that independently probe the fluctuations that source the temperature field. Maps of the CMB polarization provide exactly such information.

The obvious premise is that much of the progress expected from *LiteBIRD* in this area will stem from working with *E*-mode data (though statistical tests of the *B*-mode maps should also prove very valuable near the cosmic-variance limit). Indeed, *LiteBIRD* should almost double the statistical information concerning these anomalies. Of course, given the modest significance of the temperature anomalies, high significance detections in polarization will still prove challenging; however, only an experiment like *LiteBIRD* has the possibility of providing statistically independent information on the largest scales.

The first comprehensive search for anomalies in polarization on large angular scales was presented in Ref. [94]. No definitive evidence was found in the polarization data for anomalous features corresponding to those observed in the temperature data. Nevertheless, several tests related to dipolar modulation showed hints of asymmetry on scales up to $\ell_{max} \simeq 250$. More specifically, a variance asymmetry estimator found an alignment between the preferred directions of the temperature and *E*-mode dipolar modulation at a modest significance. However, although residual systematics did not dominate the signal, as was the case for the 2015 polarization data set [494], it was apparent that the various tests of isotropy continued to be limited by their presence, as well as the dominance of noise over signal in the large-scale *Planck* polarization data. Indeed, a notable feature of the analyses was the variation in results with the four component-separated maps studied, presumably related to their different responses to noise and systematic residuals, and an incomplete understanding of the noise properties of the data. Such effects should not be significant in the case of *LiteBIRD* since the *E*-mode signal in particular is expected to be measured at close to the cosmic-variance level, with the actual sensitivity defined by the sky fraction available for analysis.

A simple measure of whether *LiteBIRD* *E*-mode data are approaching this level of sensitivity can be provided by the inferred error on $\tau$ since, in the sample-variance limit, this should be
of order 0.002. Estimates from parameter fits to simulated *LiteBIRD* data for 70% sky coverage and including foreground residuals related to component separation indicate an error consistent with this expectation (see Sect. 6.1). Nevertheless, making inferences about the amplitude of anomalous features that might be observed in the *LiteBIRD* polarization data based on what is seen in temperature is non-trivial. Specifically, predictions must be based on models constructed in 3D position space then propagated to spherical harmonic space. This mapping is different for temperature and polarization [495].

Despite extensive work, it remains the case that no theoretical model of primordial perturbations has been constructed that can explain all of the temperature anomalies. Hence we need to consider general approaches for testing the hints of anomalies in the temperature data. This can be achieved by comparing the distribution of a specific statistic in polarization built from constrained simulations (where the part of the *E*-mode anisotropy correlated with temperature is fixed by observations of the latter) with that constructed from unconstrained realizations. For example, Ref. [496] used constraints from the *WMAP* 7-yr temperature power spectrum to compute such distributions for the *S*-statistic due to the cross-correlation between temperature and polarization. They determined that a value of the measured $S^{TQ}$ statistic over the angular separation range [48°, 120°] exceeding the value 1.403 $\mu$K$^4$ would allow the hypothesis to be ruled out at the 99% confidence level. Similarly, Ref. [497] considered the variance of the polarization amplitude when the *E*-mode signal is constrained by the *Planck* 2015 SMICA reconstruction of the CMB temperature anisotropy. In this case, the temperature data reveal an anomalously low variance in the northern ecliptic hemisphere, but the constrained realizations show no evidence of a low amplitude variance in polarization. The measurement of such a signal would argue against the temperature anomaly being merely a statistical excursion.

Irrespective of our expectations, it remains important to search for characteristic signatures of spatial anomalies in the *LiteBIRD* data, whether or not they are related to interesting features of the temperature field. Any detection of anomalies in the polarized sky signal will inevitably hint at physics beyond that captured by the standard model of cosmology.

## 6.7. Galactic astrophysics

Observations of Galactic polarization were among the main outcomes of the *Planck* space mission. Spectacular images combining the intensity of dust emission with the magnetic field orientation derived from polarization data (so-called “drapery” patterns, also called “line integral convolution” [498]) have received worldwide attention and have become part of the general scientific culture [499]. Beyond this popular impact, the *Planck* polarization maps have represented a big step forward for Galactic astrophysics [61]. We anticipate a comparable breakthrough with *LiteBIRD*, which will provide full-sky maps of Galactic polarized emission with a sensitivity many times better than that of *Planck*, both for dust and synchrotron polarization. The data will complement the rich array of other polarization observations including: (1) stellar polarization surveys to be combined with *Gaia* astrometry [500] to map Galactic dust polarization in 3D; (2) synchrotron observations, together with Faraday-rotation measurements at radio wavelengths with the Square Kilometer Array (SKA) and its precursors [501]; and (3) ground-based CMB experiments that probe smaller angular scales over restricted parts of the sky with more modest frequency coverage [502,503]. We now sketch out the expected contributions from *LiteBIRD* to two main directions of Galactic research.
<figure><img src="image_56.png" /><figcaption>**Fig. 54.** LiteBIRD's spectral coverage and sensitivity to polarized dust emission at different angular scales. This plot shows how many LiteBIRD bands will yield S/N > 1 as a function of gas column density and angular scale. For comparison, the same information is provided for Planck, with dots, stars, and hatches delineating regions where dust polarization is detected only at 353 GHz, at both 217 and 353 GHz, and at the three frequencies 143, 217 and 353 GHz, respectively.</figcaption></figure>

6.7.1. Magnetic fields. Dust polarization probes the magnetic field orientation in dusty regions, mostly in the cold and warm neutral phases of the interstellar medium (ISM) that account for the bulk of the gas mass and turbulent energy [504]. Among the various means available to map the structure of interstellar magnetic fields, dust polarization is particularly good for tracing the dynamical interplay between magnetic fields, turbulence, and gravity in the ISM. This interplay helps define the structure of interstellar matter and star formation. The multiphase magnetized ISM is too complex to be described by an analytic theory. Our understanding in this research field progresses through observations, magnetohydrodynamic (MHD) simulations, and phenomenological models. Spectroscopic observations obtained from ground-based telescopes give access to the gas column density and its kinematics. LiteBIRD will provide complementary data about magnetic fields.

LiteBIRD observations at its five highest frequencies, from 166 to 448 GHz, will improve on the Planck 353-GHz sensitivity to dust polarization by an order of magnitude (Fig. 54), increasing the dynamic range of observations by a comparable factor. While the analysis of dust polarization from the diffuse ISM at high Galactic latitudes with Planck has been limited to an effective angular resolution of 80 arcmin, LiteBIRD has the required sensitivity to map almost the whole sky down to the 17.9 arcmin beam size at 402 GHz (Fig. 55). The interplay between magnetized turbulence and gas-phase transitions builds the structure of the ISM and seeds the formation of molecular gas; however, these processes are not well understood. LiteBIRD will contribute unique polarization maps that will become the new gold standard for analyses of magnetic fields in the diffuse ISM and in the outskirts of nearby molecular clouds. In particular, we expect LiteBIRD to reveal coherent magnetic structures that result from the non-linear interplay between turbulent gas motions and magnetic fields, and the dissipation of turbulence [505,506]. Much of our current understanding of the turbulent energy cascade in
<figure><img src="image_57.png" /><figcaption>**Fig. 55.** Dust polarization images. Left: Simulated *Planck* 353-GHz polarization intensity map of the diffuse ISM at high Galactic latitude computed with the model of Ref. [508] and centered close to the Polaris Flare [509]. The local resolution is adjusted to guarantee S/N > 3 in polarization from pixel to pixel, spanning 5 arcmin (in the brightest parts) to 140 arcmin (in the more diffuse regions). The lines represent the magnetic field orientation derived from the dust polarization angle and the colors show the dust polarized intensity. Right: Simulation of the 337-GHz *LiteBIRD* polarization intensity map computed with the same model, at the native resolution of 20.9 arcmin, which guarantees S/N > 3 in the whole region.</figcaption></figure>

the ISM derives from MHD simulations that are very far from reproducing its high magnetic Reynolds number [507]. The *LiteBIRD* data will be crucial to test how well these simulations match the observations.

Ongoing developments [510–512] promise to yield powerful statistical tools for characterizing the *LiteBIRD* dust data and its comparison with simulations, ultimately leading to improved simulations that fully describe the physical processes in our Galaxy. A specific scientific objective will be to elucidate the origin of parity violation (the TB correlation [60]) of polarized dust emission. Today with the *Planck* data, we cannot decide whether this is a generic feature of interstellar turbulence rather than a random statistical fluctuation [257].

*LiteBIRD* will improve on *Planck*'s sensitivity to polarized synchrotron emission by a factor of five at 40 GHz. This gain will extend the range of scales over which the correlation between dust and synchrotron polarization is characterized [513,514]. More generally, the analysis of *LiteBIRD* dust and synchrotron data will contribute to a community effort directed towards modeling of the 3D structure of the Galactic magnetic field [515], in particular within the Solar neighborhood.

6.7.2. *Interstellar dust.* The analysis of the *LiteBIRD* data will also yield the spectral characterization of Galactic polarization. The current description of the Galactic contribution to the *Planck* and *WMAP* polarization data as arising from two components, namely thermal dust emission and synchrotron [40,60], is likely to prove inadequate for *LiteBIRD*. By providing data at 15 frequency bands between 34 and 448 GHz *LiteBIRD* will challenge our current understanding of Galactic emission. The gain is most significant for studying the nature of dust grains and the origin of the so-called anomalous microwave emission (AME) [516,517].