Tensor Non-Gaussianity Search: Current Status and Future Prospects

TL;DR

This paper reviews the current status and future prospects of detecting tensor non-Gaussianities in the CMB, emphasizing the importance of B-mode polarization data for improving detection sensitivity and guiding future surveys.

Contribution

It summarizes existing methods for tensor non-Gaussianity detection and forecasts the potential improvements with upcoming CMB B-mode polarization measurements.

Findings

No significant tensor non-Gaussianity detected so far.

Adding B-mode polarization data greatly enhances detection limits.

Established bispectrum estimation methods support future CMB B-mode surveys.

Abstract

Primordial gravitational waves (GWs) are said to be a smoking gun in cosmic inflation, while, even if they are detected, the specification of their origins are still required for establishing a true inflationary model. Testing non-Gaussianity in the tensor-mode anisotropies of the cosmic microwave background (CMB) is one of the most powerful ways to identify sources of GW signals. In this paper, we review studies searching for tensor non-Gaussianities employing the CMB bispectrum and forecast future developments. No significant signal has so far been found from temperature and E-mode polarization data, while orders-of-magnitude improvements in detection limits can be achieved by adding the information of B-mode polarization. There is already an established methodology for bispectrum estimation, which encourages a follow-up investigation with next-decadal CMB B-mode surveys.

Tables

Table 1. Table 1 : CMB constraints on three tensor non-linearity parameters: f NL t ​ t ​ t , eq superscript subscript 𝑓 NL 𝑡 𝑡 𝑡 eq f_{\rm NL}^{ttt,\rm eq} ( 6 ), f NL t ​ t ​ t , sq superscript subscript 𝑓 NL 𝑡 𝑡 𝑡 sq f_{\rm NL}^{ttt,\rm sq} ( 10 ), and f NL t ​ s ​ s , sq superscript subscript 𝑓 NL 𝑡 𝑠 𝑠 sq f_{\rm NL}^{tss,\rm sq} ( 18 ). The central values and 1 ​ σ 1 𝜎 1\sigma errors estimated from the WMAP temperature data only (first row), the Planck temperature data only (second one), the Planck E 𝐸 E -mode polarization data only (third one), and the Planck temperature and E 𝐸 E -mode polarization data jointly (fourth one) are shown. The results of f NL t ​ t ​ t , eq superscript subscript 𝑓 NL 𝑡 𝑡 𝑡 eq f_{\rm NL}^{ttt,\rm eq} are obtained by analyzing both even and odd ℓ 1 + ℓ 2 + ℓ 3 subscript ℓ 1 subscript ℓ 2 subscript ℓ 3 \ell_{1}+\ell_{2}+\ell_{3} multipoles, while those of f NL t ​ t ​ t , sq superscript subscript 𝑓 NL 𝑡 𝑡 𝑡 sq f_{\rm NL}^{ttt,\rm sq} and f NL t ​ s ​ s , sq superscript subscript 𝑓 NL 𝑡 𝑠 𝑠 sq f_{\rm NL}^{tss,\rm sq} come from only even ℓ 1 + ℓ 2 + ℓ 3 subscript ℓ 1 subscript ℓ 2 subscript ℓ 3 \ell_{1}+\ell_{2}+\ell_{3} modes.

	$f_{\mathrm{NL}}^{ttt,\mathrm{eq}}$	$f_{\mathrm{NL}}^{ttt,\mathrm{sq}}$	$f_{\mathrm{NL}}^{tss,\mathrm{sq}}$
WMAP $T$ only	$600\pm 1500$ [40]	$220\pm 170$ [39]	$84\pm 49$ [41]
Planck $T$ only	$600\pm 1600$ [44]	$290\pm 180$ [42]	–
Planck $E$ only	$2900\pm 6700$ [44]	–	–
Planck $T+E$	$800\pm 1100$ [44]	–	–

Equations

$$⟨\prod _{n=1}^3{\xi }_{{\mathbf{k}}_n}^{({\lambda }_n)}⟩=(2\pi {)}^3{\delta }^{(3)}\left(\sum _{n=1}^3{\mathbf{k}}_n\right){⟨\prod _{n=1}^3{\xi }_{{\mathbf{k}}_n}^{({\lambda }_n)}⟩}^{\prime }.$$

$${⟨\prod _{n=1}^3{\xi }_{{\mathbf{k}}_n}^{({\lambda }_n)}⟩}^{\prime }={\mathcal{B}}_{k_1{k}_2{k}_3}^{{\lambda }_1{\lambda }_2{\lambda }_3}{A}_{{\hat{k}}_1{\hat{k}}_2{\hat{k}}_3}^{{\lambda }_1{\lambda }_2{\lambda }_3}.$$

$${\mathcal{B}}_{k_1{k}_2{k}_3}^{{\lambda }_1{\lambda }_2{\lambda }_3}$$

$$A_{{\hat{k}}_1{\hat{k}}_2{\hat{k}}_3}^{{\lambda }_1{\lambda }_2{\lambda }_3}$$

$$S_{k_1{k}_2{k}_3}^{\mathrm{eq}}\equiv \frac{18}{5}{\left(2{\pi }^2{\mathcal{P}}_{\zeta }\right)}^2\left[-\left(\frac{1}{k_1^3{k}_2^3}+2\mathrm{perm}\right)-\frac{2}{k_1^2{k}_2^2{k}_3^2}+\left(\frac{1}{k_1{k}_2^2{k}_3^3}+5\mathrm{perm}\right)\right].$$

$$f_{\mathrm{NL}}^{ttt,\mathrm{eq}}\equiv \underset{k_i\to k}{\lim }\frac{{\mathcal{B}}_{k_1{k}_2{k}_3}^{+2+2+2}{A}_{{\hat{k}}_1{\hat{k}}_2{\hat{k}}_3}^{+2+2+2}}{S_{k_1{k}_2{k}_3}^{\mathrm{eq}}}.$$

$$f_{\mathrm{NL}}^{ttt,\mathrm{eq}}\sim 2.5\frac{r^2}{{\Omega }_A}.$$

$${\mathcal{B}}_{k_1{k}_2{k}_3}^{{\lambda }_1{\lambda }_2{\lambda }_3}=\sqrt{2}{f}_{\mathrm{NL}}^{ttt,\mathrm{sq}}{S}_{k_1{k}_2{k}_3}^{\mathrm{loc}},$$

$$S_{k_1{k}_2{k}_3}^{\mathrm{loc}}\equiv \frac{6}{5}{\left(2{\pi }^2{\mathcal{P}}_{\zeta }\right)}^2\left(\frac{1}{k_1^3{k}_2^3}+2\mathrm{perm}\right),$$

$$f_{\mathrm{NL}}^{ttt,\mathrm{sq}}\equiv \underset{\begin{array}{c}{k}_1\to 0\\ k_2\to k_3\end{array}}{\lim }\frac{{\mathcal{B}}_{k_1{k}_2{k}_3}^{+2+2+2}{A}_{{\hat{k}}_1{\hat{k}}_2{\hat{k}}_3}^{+2+2+2}}{S_{k_1{k}_2{k}_3}^{\mathrm{loc}}}.$$

$$f_{\mathrm{NL}}^{ttt,\mathrm{sq}}\sim {\left(\frac{B_{1\mathrm{Mpc}}}{1\mathrm{nG}}\right)}^6,$$

$$A_{{\hat{k}}_1{\hat{k}}_2{\hat{k}}_3}^{{\lambda }_1{\lambda }_{2}0}=({\hat{k}}_1\cdot {\hat{k}}_2)e_{ij}^{(-{\lambda }_1)}({\hat{k}}_1)e_{ij}^{(-{\lambda }_2)}({\hat{k}}_2),$$

$${\mathcal{B}}_{k_1{k}_2{k}_3}^{{\lambda }_1{\lambda }_{2}0}\propto \frac{{\lambda }_1}{2}{\delta }_{{\lambda }_1,{\lambda }_2}.$$

$${\mathcal{B}}_{k_1{k}_2{k}_3}^{{\lambda }_1{\lambda }_2{\lambda }_3}\propto {\delta }_{{\lambda }_1,{\lambda }_2}{\delta }_{{\lambda }_2,{\lambda }_3},$$

$${\mathcal{B}}_{k_1{k}_2{k}_3}^{{\lambda }_1{\lambda }_2{\lambda }_3}\propto \frac{{\lambda }_1}{2}{\delta }_{{\lambda }_1,{\lambda }_2}{\delta }_{{\lambda }_2,{\lambda }_3}.$$

$${\mathcal{B}}_{k_1{k}_2{k}_3}^{{\lambda }_{1}00}$$

$$A_{{\hat{k}}_1{\hat{k}}_2{\hat{k}}_3}^{{\lambda }_{1}00}$$

$$f_{\mathrm{NL}}^{tss,\mathrm{sq}}\equiv \underset{\begin{array}{c}{k}_1\to 0\\ k_2\to k_3\end{array}}{\lim }\frac{{\mathcal{B}}_{k_1{k}_2{k}_3}^{+200}{A}_{{\hat{k}}_1{\hat{k}}_2{\hat{k}}_3}^{+200}}{S_{k_1{k}_2{k}_3}^{\mathrm{loc}}}.$$

$$f_{\mathrm{NL}}^{tss,\mathrm{sq}}\sim 0.1r{\lambda }_{sst},$$

$$X(\hat{n})=\sum _{\ell m}{a}_{\ell m}^X{Y}_{\ell m}(\hat{n}),$$

$$a_{\ell m}^X=a_{\ell m}^{X(s)}+a_{\ell m}^{X(v)}+a_{\ell m}^{X(t)}.$$

$$a_{\ell m}^{X(z)}=4\pi (-i{)}^{\ell }\int \frac{d^3\mathbf{k}}{(2\pi {)}^3}{\mathcal{T}}_{\ell ,k}^{X(z)}\sum _{\lambda }[\mathrm{sgn}(\lambda ){]}^{\lambda +x}{\xi }_{\mathbf{k}}^{(\lambda )}{}_{-\lambda }{Y}_{\ell m}^{\ast }(\hat{k}),$$

$$⟨\prod _{n=1}^3{a}_{{\ell }_n{m}_n}^{X_n(z_n)}⟩=\left[\prod _{n=1}^3(-i{)}^{{\ell }_n}\int \frac{d^3{\mathbf{k}}_n}{2{\pi }^2}{\mathcal{T}}_{{\ell }_n,k_n}^{X_n(z_n)}\sum _{{\lambda }_n}[\mathrm{sgn}({\lambda }_n){]}^{{\lambda }_n+x_n}{}_{-{\lambda }_n}{Y}_{{\ell }_n{m}_n}^{\ast }({\hat{k}}_n)\right]⟨\prod _{n=1}^3{\xi }_{{\mathbf{k}}_n}^{({\lambda }_n)}⟩.$$

$${\xi }_{\mathbf{k}}^{(\lambda )}=\sum _{\ell m}{\xi }_{\ell m}^{(\lambda )}(k){}_{-\lambda }{Y}_{\ell m}(\hat{k}),$$

$$⟨\prod _{n=1}^3{a}_{{\ell }_n{m}_n}^{X_n(z_n)}⟩=\left[\prod _{n=1}^3(-i{)}^{{\ell }_n}{\int }_0^{\infty }\frac{k_n^{2}d{k}_n}{2{\pi }^2}{\mathcal{T}}_{{\ell }_n,k_n}^{X_n(z_n)}\sum _{{\lambda }_n}[\mathrm{sgn}({\lambda }_n){]}^{{\lambda }_n+x_n}\right]⟨\prod _{n=1}^3{\xi }_{{\ell }_n{m}_n}^{({\lambda }_n)}(k_n)⟩,$$

$$⟨\prod _{n=1}^3{\xi }_{{\ell }_n{m}_n}^{({\lambda }_n)}(k_n)⟩=\left[\prod _{n=1}^3\int d^2{\hat{k}}_n{}_{-{\lambda }_n}{Y}_{{\ell }_n{m}_n}^{\ast }({\hat{k}}_n)\right]⟨\prod _{n=1}^3{\xi }_{{\mathbf{k}}_n}^{({\lambda }_n)}⟩.$$

$$\Braket{\prod_{n=1}^{3}a^{X_{n}(z_{n})}_{\ell_{n}m_{n}}}=\left(\begin{array}[]{ccc}\ell_{1}&\ell_{2}&\ell_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right)B_{(z_{1}z_{2}z_{3})\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}.$$

$$|{\ell }_1-{\ell }_2|\le {\ell }_3\le {\ell }_1+{\ell }_2.$$

$$e_{ij}^{(-{\lambda }_1)}({\hat{k}}_1)e_{jk}^{(-{\lambda }_2)}({\hat{k}}_2)e_{ki}^{(-{\lambda }_3)}({\hat{k}}_3)$$

$${\delta }^{(3)}\left(\sum _{n=1}^3{\mathbf{k}}_n\right)$$

$$h^{s_{1}s_{2}s_{3}}_{l_{1}l_{2}l_{3}}\equiv\sqrt{\frac{(2l_{1}+1)(2l_{2}+1)(2l_{3}+1)}{4\pi}}\left(\begin{array}[]{ccc}l_{1}&l_{2}&l_{3}\\ s_{1}&s_{2}&s_{3}\end{array}\right).$$

$$\int d^{2}\hat{k}_{n}\,{}_{-\lambda_{n}}Y_{\ell_{n}m_{n}}^{*}(\hat{k}_{n})\,Y_{L_{n}M_{n}}^{*}(\hat{k}_{n})\,{}_{\lambda_{n}}Y_{2\mu_{n}}^{*}(\hat{k}_{n})=h_{\ell_{n}L_{n}2}^{\lambda_{n}0-\lambda_{n}}\left(\begin{array}[]{ccc}\ell_{n}&L_{n}&2\\ m_{n}&M_{n}&\mu_{n}\end{array}\right).$$

Full text

\helveticabold

1 Keywords:

CMB, non-Gaussianity, gravitational wave, inflation, observational constraints, B-mode

2 Introduction

Recent cosmic-variance-limited-level measurements of the scalar sector in primordial fluctuations using the temperature and $E$-mode polarization fields of the cosmic microwave background (CMB), indicate that the Universe experienced an inflationary expansion at very early stages [1, 2]. Now, compelling evidence is expected to lie in the tensor sector; namely, the primordial gravitational wave (GW). Even without any particular source, primordial GWs arise naturally from inflationary vacuum fluctuations. Since their amplitude directly reflects the inflationary energy scale, various observational projects aim at hunting them. There, the measurements of the large-scale CMB $B$-mode polarization have attracted everyone’s attention as it is a distinctive observable of the tensor mode (see e.g., Refs. [3, 4] and references therein).

In order to investigate the primordial tensor sector using the CMB anisotropies, the power spectrum (2-point correlation) is primarily employed. In single-field inflation with Einstein gravity, GWs are nearly Gaussian and hence the higher-order correlations vanish [5, 6]. In this case, only the power spectrum becomes an informative observable. This paper, however, will address the potential of the bispectrum (3-point correlation) as a probe of inflationary physics.

One possible way, making it more informative, is the addition of some extra source fields. In the presence of the gauge field, the anisotropic stress fluctuations are formed, generating GWs. Because this is a non-linear process, the resulting GWs are non-Gaussian (NG) and therefore a large tensor bispectrum can be realized. In particular, the production of NG GWs becomes more efficient when the gauge field is coupled to the axion [7, 8]. Moreover, production due to higher-spinning particles has been studied intensively in recent years [9, 10]. There is also the possibility of a post-inflationary generation of NG GWs due to magnetic fields [11]. Another possible way to produce large tensor NGs comes from non-trivial non-linear gravitational interactions predicted in some modified theories of gravity (e.g. [12, 13, 14, 15, 16, 17, 18, 19, 20]). Then, the induced bispectra can have distinctive shapes compared with that from the usual Einstein term.

In this sense, the CMB tensor-mode bispectrum is a key indicator of the inflationary particle content and/or high energy gravity. Since a general formalism was established in Refs. [21, 22], various signatures of various theoretical models in the CMB bispectra have been studied (see Refs. [11, 23, 24, 25] for primordial magnetic field (PMF) models, Refs. [26, 7, 27] for the axion inflation ones, Refs. [22, 28, 16, 29, 30] for modified gravity ones and Refs. [31, 32] for other motivations). One of the most interesting phenomena discovered there (unseen in the scalar bispectrum analysis) is that the tensor bispectrum can yield the non-vanishing signal in not only even but also odd $\ell_{1}+\ell_{2}+\ell_{3}$ multipoles. The GW bispectra generate the CMB auto and cross bispectra including the $B$-mode polarization. In the usual theories where the GW bispectra are parity-invariant, the odd $\ell_{1}+\ell_{2}+\ell_{3}$ signal can arise in $TTB$, $TEB$, $EEB$, and $BBB$ [25, 31, 16, 29]. Moreover, parity-odd GW bispectra motivated by parity-breaking theories can also induce the non-zero odd $\ell_{1}+\ell_{2}+\ell_{3}$ signal in $TTT$, $TTE$, $TEE$, $TBB$, $EEE$, and $EBB$ [28, 24, 26, 33, 27, 30].

Now, any type of theoretical bispectrum template is testable with the temperature and E/B-mode polarization data by use of a general bispectrum estimator for the even [34, 35, 36, 37] and odd $\ell_{1}+\ell_{2}+\ell_{3}$ domain [38, 37]. The magnitudes of some representative tensor NG templates have already been estimated from the temperature map of WMAP [39, 40, 41] and the temperature and $E$-mode polarization from Planck [42, 43, 44]. Absence of NG reported there constrains the various inflationary models mentioned above.

In this paper, besides reviewing previous studies, we discuss future prospects of the tensor NG search assuming the detection of $B$-mode polarization in next-generation CMB experiments. We found that, in terms of observational limits on the size of some scale-invariant templates, an improvement up to three orders of magnitude is expected assuming LiteBIRD-level sensitivities [45, 46, 47].

This paper is organized as follows. In the next section, we summarize theoretical scenarios predicting characteristic primordial tensor bispectra. In Sec. 4, we compute the CMB bispectra and show their behavior. In Sec. 5, we present the current observational constraints on tensor NGs and discuss future prospects assuming $B$-mode polarization surveys. The final section is devoted to the conclusions of this paper.

3 Theoretical motivations

In this section, we briefly review the primordial tensor bispectra predicted in some theoretical scenarios. For convenience, let us work with the helicity basis; thus, the primordial perturbation is represented as $\xi_{\bf k}^{(\lambda)}$ [22]. Here, $\lambda$ denotes the helicity and takes [math], $\pm 1$, and $\pm 2$ for the scalar, vector, and tensor modes, respectively. In practice, we identify $\xi_{\bf k}^{(0)}$ and $\xi_{\bf k}^{(\pm 2)}$ with the curvature perturbation $\zeta_{{\bf k}}$ and the GW $h_{{\bf k}}^{(\pm 2)}$, respectively.

Assuming statistical homogeneity, the primordial bispectrum takes the form

[TABLE]

Then, let us split this in two

[TABLE]

In what follows, how the shape, represented by $\mathcal{B}_{k_{1}k_{2}k_{3}}^{\lambda_{1}\lambda_{2}\lambda_{3}}$, and the angular structure, represented by $A_{\hat{k}_{1}\hat{k}_{2}\hat{k}_{3}}^{\lambda_{1}\lambda_{2}\lambda_{3}}$, changes depending on the model is argued.

3.1 Extra sources

Efficient tensor NG production is realized by adding extra source fields to the theory. In what follows, we discuss the vector field as a source.

3.1.1 Inflationary axion-gauge coupling

A characteristic tensor NG is realized in inflationary models involving the coupling between the axion $\phi$ and a gauge field like $f(\phi)\widetilde{F}F$, with $F$ and $\widetilde{F}$ the gauge field strength tensor and its dual, respectively [48, 7, 49, 50, 8, 51]. In this case, the chirality of the gauge field is transferred into the GW sector and hence only the plus mode of GW, $h^{(+2)}$, survives. The resulting production is efficient at specific scales, and the tensor-tensor-tensor bispectrum peaks for equilateral configurations, $k_{1}\sim k_{2}\sim k_{3}$ [48].

The detailed spectral feature varies with the shape of the coupling $f(\phi)$ or the axion potential [7, 51]. For example, adopting $f(\phi)\propto\phi$ and a nearly flat potential, one can obtain a scale-invariant shape as [26]

[TABLE]

where $e_{ij}^{(\lambda)}(\hat{k})$ is the polarization tensor defined in the helicity basis obeying $e_{ii}^{(\lambda)}(\hat{k})=\hat{k}_{i}e_{ij}^{(\lambda)}(\hat{k})=0$, $e_{ij}^{(\lambda)*}(\hat{k})=e_{ij}^{(-\lambda)}(\hat{k})=e_{ij}^{(\lambda)}(-\hat{k})$, and $e_{ij}^{(\lambda)}(\hat{k})e_{ij}^{(\lambda^{\prime})}(\hat{k})=2\delta_{\lambda,-\lambda^{\prime}}$ [22], and $S_{k_{1}k_{2}k_{3}}^{\rm eq}$ is the usual scalar equilateral bispectrum template, given by:

[TABLE]

Here $\mathcal{P}_{\zeta}$ is the amplitude of the curvature power spectrum. A newly introduced parameter $f_{\rm NL}^{ttt,\rm eq}$ quantifies the relative size of the $\lambda=+2$ bispectrum to $S_{k_{1}k_{2}k_{3}}^{\rm eq}$ and satisfies111This $f_{\rm NL}^{ttt,\rm eq}$ is equivalent to $f_{\rm NL}^{P}$ in Ref. [40] and $f_{\rm NL}^{\rm tens}$ in Ref. [43, 44].

[TABLE]

The expected size of the GW bispectrum, of course, depends on the case. If the axion is identified with the inflaton field, the amount of scalar production exceeds the GW one [52]. Then, the scalar bispectrum measurement gives much tighter constraints on the model [43, 44]. In contrast, in multifield models where the inflaton field and the axion coexist, more effective GW production occurs. In the model where the axion is coupled to the U(1) gauge field, a characteristic bump appears in the GW bispectrum, and, depending on its location, it is detectable by the CMB $BBB$ bispectrum measurement [7, 27]. However, at the same time, a similar bump in the GW power spectrum is measured with a higher signal-to-noise ratio from the CMB $BB$ power spectrum. The strongest constraints are obtained through the GW power spectrum measurement, and the GW bispectrum provides complementary information. In the model including a SU(2)-gauge coupling, the GW bispectrum can dominate the scalar one, and more interestingly, in some regions of the parameter space, the GW bispectrum has high detectability compared to the GW power spectrum [8, 51]. The GW bispectrum is nearly scale-invariant and therefore parametrized by $f_{\rm NL}^{ttt,\rm eq}$, which is related to the tensor-to-scalar ratio $r$ and the energy density fraction of the SU(2) gauge field $\Omega_{A}$ according to [8]

[TABLE]

This indicates that $f_{\rm NL}^{ttt,\rm eq}$ can take detectably large values, i.e., $f_{\rm NL}^{ttt,\rm eq}>1$, when $\Omega_{A}$ is smaller than $r^{2}$.

3.1.2 Primordial magnetic fields

If PMFs, $B_{i}$, are created via some non-trivial mechanism before inflation, they are stretched beyond the horizon by the inflationary expansion and there remain relic fields at large scales. Such a scenario has been widely argued to explain galactic or extragalactic magnetic fields observed at present (for reviews see Refs. [53, 54]). In this case, PMFs form the anisotropic stress fluctuations, and they source the GW at superhorizon scales until neutrino decoupling. The induced GW is called a passive mode and behaves as the initial condition of the CMB fluctuations after reentering the horizon [55].

Assuming Gaussianity of PMFs, the induced GW, $h_{ij}\propto(B_{i}B_{j})^{TT}$, becomes a chi-square random field [56]. When the PMF power spectrum has a scale-invariant shape, in analogy with the usual local type NG, the tensor-tensor-tensor bispectrum has the squeezed shape [11, 23]:

[TABLE]

where the usual scalar local bispectrum template reads:

[TABLE]

and $A_{\hat{k}_{1}\hat{k}_{2}\hat{k}_{3}}^{\lambda_{1}\lambda_{2}\lambda_{3}}$ is given in the same form as Eq. (4). A newly introduced parameter $f_{\rm NL}^{ttt,\rm sq}$ represents the relative size of the $\lambda=+2$ bispectrum to $S_{k_{1}k_{2}k_{3}}^{\rm loc}$ and satisfies:

[TABLE]

As the GW bispectrum has a sextic dependence on PMFs, $f_{\rm NL}^{ttt,\rm sq}$ takes the form:

[TABLE]

where $B_{1\,\rm Mpc}$ is the PMF strength smoothed on $1\,{\rm Mpc}$. This means that $f_{\rm NL}^{ttt,\rm sq}>1$ is realized if nano-Gauss-level PMFs exist. Note that, in the same mechanism, the scalar local type bispectrum is generated, however, it is subdominant compared with the GW one [25].

This is the result under the absence of a helical (parity-odd) term in the PMF power spectrum, however, the bispectrum shape is modified if it exists [24].

3.2 Modified gravity

The above scenarios rely on some extra source fields. Let us discuss another NG GW production through the modification of Einstein gravity.

3.2.1 Weyl gravity

In the context of the Weyl gravity, there may exist the Chern-Simons term $f(\phi)\widetilde{W}W$, where $W$ is the Weyl tensor and $\widetilde{W}$ is its dual. Like the axion inflation models, this term also sources chiral GWs [57, 58]. The tensor-tensor-scalar bispectrum, whose angular structure reads:

[TABLE]

surpasses the tensor-tensor-tensor and tensor-scalar-scalar ones in amplitude [17]. Due to the parity-violating nature in the GW sector, the helical bispectrum obeys

[TABLE]

One can also construct the cubic actions $f(\phi)W^{3}$ and $f(\phi)\widetilde{W}W^{2}$ [12, 59, 28]. The tensor-tensor-tensor bispectra generated from both terms have identical angular dependence as Eq. (4). In contrast, the helicity dependence is completely different [28]: the former gives rise to the non-helical (parity-even) contribution:

[TABLE]

while the latter induces the helical (parity-odd) one

[TABLE]

3.2.2 Massive gravity

In single-field slow-roll inflation with Einstein gravity, the tensor-scalar-scalar bispectrum is affected by slow-roll suppression [6]. However, in a massive gravity model [16], the size of the bispectrum is controlled by the graviton mass. The bispectrum, taking the form:

[TABLE]

with $k_{t}\equiv k_{1}+k_{2}+k_{3}$, is amplified at long-wavelength tensor and short-wavelength scalar configurations ($k_{1}\ll k_{2}\sim k_{3}$). A newly introduced parameter $f_{\rm NL}^{tss,\rm sq}$ denotes the relative size of the $\lambda=+2$ bispectrum to $S_{k_{1}k_{2}k_{3}}^{\rm loc}$, obeying:

[TABLE]

With an upper bound on the tensor-to-scalar ratio $r\lesssim 0.1$, the prediction in Einstein gravity, $f_{\rm NL}^{tss,\rm sq}\sim 0.1\,r$, indicates $f_{\rm NL}^{tss,\rm sq}\lesssim 0.01$. In contrast, massive gravity model modifies this as:

[TABLE]

and $f_{\rm NL}^{tss,\rm sq}>1$ can then be realized depending on the strength parameter of a non-linear interaction $\lambda_{sst}$. Note that the spectral index of the GW power spectrum is also a possible observable of this model [16], but has been unconstrained so far.

4 CMB bispectra from tensor non-Gaussianities

Next, let us discuss the CMB bispectra generated from GW NGs.

4.1 General formalism for CMB angular bispectrum

We start by reviewing how to compute the CMB bispectra generated from the primordial scalar, vector, and tensor NGs based on the general formalism developed in Refs. [21, 22].

The CMB field is characterized by temperature and two linear polarizations called $E$ and $B$ modes. All of these are spin-0 fields, and the temperature and $E$-mode fields have even parity, while the $B$ mode is parity-odd. Since temperature ($X=T$), $E$-mode ($X=E$), and $B$-mode ($X=B$) fields are distributed on the 2D sphere, they can be expanded using the spherical harmonic basis according to

[TABLE]

where we have assumed the spatial flatness of the Universe. The spherical harmonic coefficients $a_{\ell m}^{X}$ are given by the sum of the scalar, vector, and tensor modes as:

[TABLE]

Note that the scalar mode cannot generate a $B$-mode signal due to zero helicity. The coefficients of the scalar ($z=s$), vector ($z=v$), and tensor ($z=t$) modes have the structure:

[TABLE]

where $x$ changes depending on the parity of the CMB field according to $x=0$ for $X=T,E$ and $x=1$ for $X=B$, and recall that $\lambda$ represents helicity as $\lambda=0$ for $z=s$, $\lambda=\pm 1$ for $z=v$, and $\lambda=\pm 2$ for $z=t$. The linear transfer function $\mathcal{T}_{\ell,k}^{X(z)}$ represents the time evolution of the CMB fluctuations originating from the primordial perturbation $\xi_{\bf k}^{(\lambda)}$. Using Eq. (22), one can form the CMB bispectrum as:

[TABLE]

The primordial curvature perturbation and the primordial GW act as initial conditions of the scalar and tensor CMB anisotropies, respectively; thus, $\xi_{{\bf k}}^{(0)}=\zeta_{{\bf k}}$ and $\xi_{{\bf k}}^{(\pm 2)}=h_{{\bf k}}^{(\pm 2)}$. Employing the harmonic expansion:

[TABLE]

the above formula is rewritten as:

[TABLE]

where,

[TABLE]

This is further simplified once the explicit formula of $\Braket{\prod_{n=1}^{3}\xi_{{\bf k}_{n}}^{(\lambda_{n})}}$ is given. The computation procedure is as follows:

Expand $\Braket{\prod_{n=1}^{3}\xi_{{\bf k}_{n}}^{(\lambda_{n})}}$ using the (spin-weighted) spherical harmonics in terms of $\hat{k}_{1}$, $\hat{k}_{2}$, and $\hat{k}_{3}$. 2. 2.

Perform the angular integrals, $\int d^{2}\hat{k}_{1}$, $\int d^{2}\hat{k}_{2}$, and $\int d^{2}\hat{k}_{3}$, of the resultant spherical harmonics and convert them into products of Wigner symbols. 3. 3.

Simplify the resultant products of the Wigner symbols by adding angular momenta.

If the primordial bispectrum $\Braket{\prod_{n=1}^{3}\xi_{{\bf k}_{n}}^{(\lambda_{n})}}$ respects rotational invariance, via the above computation, the Wigner $3j$ symbol \left(\begin{array}[]{ccc}\ell_{1}&\ell_{2}&\ell_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right) is singled out in $\Braket{\prod_{n=1}^{3}\xi_{\ell_{n}m_{n}}^{(\lambda_{n})}(k_{n})}$ and hence the resulting CMB bispectrum takes the form:

[TABLE]

In this case, the nonzero signal is confined to the tetrahedral domain:

[TABLE]

4.2 Practical examples

In the following, we demonstrate the CMB bispectrum computation in the context of the practical cases discussed in Sec. 3.

4.2.1 Three tensors

First, we discuss the CMB bispectra sourced by primordial tensors where we examine the axion model $f_{\rm NL}^{ttt,\rm eq}$ template (3) [26], the $W^{3}$ model (14) [28], the $\widetilde{W}W^{2}$ model (15) [28], and the PMF model $f_{\rm NL}^{ttt,\rm sq}$ (8) [11]. Since $A_{\hat{k}_{1}\hat{k}_{2}\hat{k}_{3}}^{\lambda_{1}\lambda_{2}\lambda_{3}}$ takes the identical form, $\Braket{\prod_{n=1}^{3}h_{\ell_{n}m_{n}}^{(\lambda_{n})}(k_{n})}$ can be computed in a similar way.

For the first step, in order to simplify the R.H.S. of Eq. (26), we expand $A_{\hat{k}_{1}\hat{k}_{2}\hat{k}_{3}}^{\lambda_{1}\lambda_{2}\lambda_{3}}$ and the Dirac delta function:

[TABLE]

where

[TABLE]

For the second step, the angular integrals of the resultant spherical harmonics are performed, such that:

[TABLE]

As a final step, the summation of the Wigner $3j$ symbols appearing above over angular momenta is computed according to:

[TABLE]

Consequently, we obtain:

[TABLE]

Here, the prefactor \left(\begin{array}[]{ccc}\ell_{1}&\ell_{2}&\ell_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right) is due to rotational invariance of the primordial bispectrum. Inserting this into Eq. (25) and simplifying the resultant equation, we find that the CMB bispectra take the rotational-invariant form (27) with:

[TABLE]

The $\ell$-space domain containing the non-zero signal is determined by the summation over the helicities $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$. In the $\widetilde{W}W^{2}$ model, performing the summation of Eq. (15) leads to the non-vanishing condition: [28]

[TABLE]

On the other hand, in the $W^{3}$ and PMF models, due to the helicity dependence in Eqs. (14) and (8), the non-vanishing signal is confined to [11, 28]:

[TABLE]

In the axion model, however, the helicity dependence in Eq. (3) does not yield any restriction, allowing a non-vanishing signal in the whole tetrahedral domain [26].

The left and center panels of Fig. 1 show the intensity distributions of the temperature bispectra from the axion model $f_{\rm NL}^{ttt,\rm eq}$ template (3) and the PMF model $f_{\rm NL}^{ttt,\rm sq}$ template (8), respectively. As expected, it is confirmed that the signal comes mostly from the equilateral ($\ell_{1}\sim\ell_{2}\sim\ell_{3}$) and squeezed ($\ell_{1}\ll\ell_{2}\sim\ell_{3}$, $\ell_{2}\ll\ell_{3}\sim\ell_{1}$ and $\ell_{3}\ll\ell_{1}\sim\ell_{2}$) configurations, respectively. The decaying nature for $\ell\gtrsim 100$ due to the lack of the tensor-mode integrated Sachs-Wolfe contribution [11, 26] is also visually apparent.

4.2.2 Two tensors and one scalar

Here, we compute the CMB tensor-tensor-scalar bispectrum from the $\widetilde{W}W$ model (13) [30]. The angular dependence in $A_{\hat{k}_{1}\hat{k}_{2}\hat{k}_{3}}^{\lambda_{1}\lambda_{2}0}$ (12) is decomposed according to:

[TABLE]

By means of the above methodology, we obtain [30]:

[TABLE]

Performing the summation over $\lambda_{1}$ and $\lambda_{2}$ with Eq. (13), we find that the non-vanishing signal obeys [30]:

[TABLE]

4.2.3 One tensor and two scalars

We here compute the CMB bispectrum from the massive gravity model $f_{\rm NL}^{tss,\rm sq}$ template (16) [16]. The spherical harmonic expansion of $A_{\hat{k}_{1}\hat{k}_{2}\hat{k}_{3}}^{\lambda_{1}~{}0~{}0}$ (17) then reads:

[TABLE]

In the same manner, we can derive the CMB bispectrum as [22, 16]:

[TABLE]

By summing over $\lambda_{1}$, we confirm that, because of the absence of the helicity dependence in Eq. (16), the non-vanishing signal is distributed on [22, 16]:

[TABLE]

The right panel of Fig. 1 describes the intensity distribution of the resultant temperature bispectrum showing explicitly that it is amplified in the squeezed limit. Differently from the $f_{\rm NL}^{ttt,\rm sq}$ case, the squeezed signal survives even for $\ell\gtrsim 100$ because of the presence of the scalar-mode acoustic oscillation.

5 Measurements of CMB tensor-mode bispectra

In this section, we examine the present and future constraints on GW NGs through the CMB bispectrum measurements.

5.1 Optimal tensor-mode bispectrum estimation

Let us begin by explaining how to estimate the overall magnitude of the primordial bispectrum, dubbed as $f_{\rm NL}$, from the CMB data.

An optimal $f_{\rm NL}$ estimator takes the form [60]:

[TABLE]

where $C_{\ell}$ is the CMB angular power spectrum, $B_{\ell_{1}\ell_{2}\ell_{3}}^{\rm th}$ is the theoretical template of the angle-averaged CMB bispectrum one wants to measure [corresponding to $B_{(z_{1}z_{2}z_{3})\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}$ in Eq. (27) for $f_{\rm NL}=1$],

[TABLE]

is the bispectrum reconstructed from the observed $a_{\ell m}$, and

[TABLE]

is the Fisher matrix. The bestfit value and the error bar on $f_{\rm NL}$ is obtained computing $\hat{f}_{\rm NL}$ from the observed data and the simulation maps, respectively. The straightforward computation of Eq. (81) is not a wise manner since the $\mathcal{O}(\ell_{\rm max}^{5})$ numerical operation required here is quite time-consuming and practically unfeasible for $\ell_{\rm max}\gtrsim 100$.

To reduce the computational cost, factorizing the summation as $\sum_{\ell_{1}\ell_{2}\ell_{3}}f(\ell_{1},\ell_{2},\ell_{3})=[\sum_{\ell_{1}}a(\ell_{1})][\sum_{\ell_{2}}b(\ell_{2})][\sum_{\ell_{3}}c(\ell_{3})]$ is a tested approach.222 See Ref. [61] for another approach by binning the $\ell$-space domain.

This idea was implemented in Ref. [60] for the first time, enabling the analysis for originally separable templates such as local, equilateral and orthogonal NGs. Later, it was generalized by developing the modal decomposition technique of any originally non-separable $B_{\ell_{1}\ell_{2}\ell_{3}}^{\rm th}$ [34, 35]. The original version of the modal methodology is for the analysis of even $\ell_{1}+\ell_{2}+\ell_{3}$ domain, while the spin-weighted version developed recently [38, 40] can also cover odd $\ell_{1}+\ell_{2}+\ell_{3}$ domain. As seen in Sec. 4, the tensor-mode CMB bispectrum is generally non-separable and does not always vanish in the odd $\ell_{1}+\ell_{2}+\ell_{3}$ domain. In this sense, the modal methodology is an indispensable tool for the tensor NG search. In fact, except for the results obtained bruteforcing the $\ell$-space summation with $\ell_{\rm max}=100$ [39], all observational constraints reported so far have been obtained by it [40, 43, 42, 41, 44]. In what follows, for simplicity, we explain the modal methodology for the auto bispectra ($TTT$, $EEE$, and $BBB$) based on Refs. [34, 35, 38, 40], while one can deal with the cross bispectra ($TTE$, $TTB$, $TEE$, $TEB$, $TBB$, $EEB$, and $EBB$) in the very similar way [36, 37].

Let us introduce the reduced bispectrum according to:

[TABLE]

where,

[TABLE]

Here and hereinafter, we follow the notation for the symmetric operation $A^{\{x}A^{y}A^{z\}}\equiv\frac{1}{6}A^{x}A^{y}A^{z}+{\rm 5~{}perm~{}in~{}}x,y,z$. The spinned weight $(x,y,z)$ is fixed to be, e.g., $(0,0,0)$ and $(1,1,-2)$ for even and odd $\ell_{1}+\ell_{2}+\ell_{3}$ analysis, respectively. The reduced bispectrum of the theoretical template is decomposed in the odd/even $\ell_{1}+\ell_{2}+\ell_{3}$ domain separately, according to:

[TABLE]

where the separable modal basis $Q_{ijk}$ is composed of the products of the eigenfunctions $q_{i}(\ell)\in\mathbb{R}$ as:

[TABLE]

The $\sigma$ factor, defined by:

[TABLE]

makes the L.H.S of Eq. (86) real, so the modal coefficients $\alpha_{ijk}^{(o/e)}$ are always real. The convergence speed of the modal decomposition relies on the choice of the $v_{\ell}$ weighting and $q_{i}(\ell)$. Employing the $Q_{ijk}$ templates, we define a matrix according to:

[TABLE]

where the triples $ijk$ in $Q_{ijk}$ are labeled by means of a single index $n$, and the bracket $\braket{\cdots}_{o/e}$ denotes the summation in terms of the odd/even $\ell_{1}+\ell_{2}+\ell_{3}$ domain as:

[TABLE]

Then, the modal coefficients are computed according to the inverse operation:

[TABLE]

Using Eq. (86) and the identity:

[TABLE]

the sum over $\ell_{1}$, $\ell_{2}$, and $\ell_{3}$ in the estimator (81) is rewritten into the sum over finite eigenmodes as:

[TABLE]

where the $\beta_{n}^{(o/e)}$ coefficients, reading:

[TABLE]

are computed from the filtered maps:

[TABLE]

The most time-consuming task is the computation of $\braket{\cdots}_{o/e}$ in $\alpha_{n}^{(o/e)}$, however, it requires, at most, $\mathcal{O}(\ell_{\rm max}^{3})$ operations (assuming that the modal decomposition converges within a reasonable time), and hence the separable estimator (93) makes the tensor bispectrum analysis feasible.

5.2 Current observational limits

Table 1 summarizes the current CMB limits on $f_{\rm NL}^{ttt,\rm eq}$, $f_{\rm NL}^{ttt,\rm sq}$, and $f_{\rm NL}^{tss,\rm sq}$ obtained from the WMAP and ${\it Planck}$ maps.333See Ref. [40] for the CMB limits on a few other shapes. Regarding $f_{\rm NL}^{ttt,\rm eq}$ and $f_{\rm NL}^{ttt,\rm sq}$, the consistency between the WMAP $T$ only results and Planck $T$ only ones is confirmed there. Note that the unimproved constraints are expected because, as inferred from Fig. 1, $B_{\ell_{1}\ell_{2}\ell_{3}}^{\rm th}$ decays rapidly for $\ell\gtrsim 100$ and hence the estimator sum (81) saturates well below the WMAP resolution ($\ell_{\rm max}=500$). In contrast, for $f_{\rm NL}^{tss,\rm sq}$ the current constraint from WMAP is expected to be much improved by Planck since the squeezed-limit signal survives even for high $\ell$ [22, 16]. We also confirm the improvement of the constraint on $f_{\rm NL}^{ttt,\rm eq}$ by analyzing the temperature and $E$-mode polarization data jointly.

As seen in Table 1, no significant deviation from Gaussianity has been found so far.

5.3 Future prospects

Finally, we discuss future prospects of detecting these tensor NGs by adding $B$-mode polarization to the data analysis. Here, we evaluate expected $1\sigma$ errors $\Delta f_{\rm NL}^{ttt,\rm eq}$ and $\Delta f_{\rm NL}^{ttt,\rm sq}$ from $BBB$, and $\Delta f_{\rm NL}^{tss,\rm sq}$ from $BTT$, through the computation of the Fisher matrix. The covariance matrices [corresponding to the denominator of Eq. (83)] are given by $C_{\ell_{1}}^{BB}C_{\ell_{2}}^{BB}C_{\ell_{3}}^{BB}$ and $C_{\ell_{1}}^{BB}C_{\ell_{2}}^{TT}C_{\ell_{3}}^{TT}$, respectively. Note that these are estimated under the diagonal covariance matrix approximation. The $B$-mode power spectrum is computed by summing up the primordial contribution parametrized by the tensor-to-scalar ratio $r$, the lensing $B$-mode one and the experimental noise spectrum as $C_{\ell}^{BB}=C_{\ell}^{\rm prim}(r)+C_{\ell}^{\rm lens}+N_{\ell}$. Then, let us examine three cleanliness levels of the $B$-mode data: a perfectly-delensed and noiseless full-sky case (i.e., $C_{\ell}^{\rm lens}=N_{\ell}=0$), a non-delensed and noiseless full-sky one (i.e., $N_{\ell}=0$), and a non-delensed realistic case where experimental uncertainties due to beam, noise, mask and residual foreground in LiteBIRD [45, 46, 47] are considered.444 In the third case, we compute $N_{\ell}$ by taking into account the contamination due to residual foregrounds in the same manner as Ref. [27]; namely, we assume that foregrounds due to galactic dust emission and synchrotron radiation are subtracted using 9 channels (corresponding to 40-89 GHz and 280-402 GHz) and therefore reduce to the 2% level in CMB maps. Any higher-order contribution is not considered for simplicity.

The first example provides the theoretical limits. In this case, the size of the covariance matrix is determined by $r$ alone and, accordingly, the errors simply scale like:

[TABLE]

Comparing this with Eq. (7) or Eq. (19), one can see that the error on a model parameter of the axion model, $\Omega_{A}^{-1}$, or that of the massive gravity model, $\lambda_{sst}$, is proportional to $r^{-1/2}$.

The left panel of Fig. 2 describes the $\ell_{\rm max}$ dependence of $\Delta f_{\rm NL}$ in the perfectly-delensed and non-delensed cases with $r=10^{-3}$ and $N_{\ell}=0$. From the former results (corresponding to the dot-dashed lines), the sensitivity improvement by increasing $\ell_{\rm max}$ is confirmed. From the latter results (corresponding to the solid lines), it is visually apparent that $\Delta f_{\rm NL}^{ttt,\rm eq}$ and $\Delta f_{\rm NL}^{ttt,\rm sq}$ immediately saturate because $C_{\ell}^{\rm lens}$ dominates over $C_{\ell}^{BB}$ for $\ell>\mathcal{O}(10)$. In contrast, $\Delta f_{\rm NL}^{ttt,\rm sq}$ is free from such a degradation since the dominant signal in the Fisher matrix comes from large-scale $B$ modes, more precisely, long-wavelength $B$ and short-wavelength $T$ squeezed configurations ($\ell_{1}^{B}\ll\ell_{2}^{T}\sim\ell_{3}^{T}$) [16, 31]. However, of course, $\Delta f_{\rm NL}^{ttt,\rm sq}$ would also saturate at very small scales, i.e., $\ell^{T}\gtrsim 3000$, where the scalar-mode lensing contamination dominates $C_{\ell}^{TT}$. Higher-order lensing contributions introduce non-vanishing off-diagonal components in the covariance matrix; thus, the above simple Fisher matrix computation would no longer be credible. In the left panel of Fig. 2, the results of the LiteBIRD-like experiment are not shown, however, they have similar $\ell_{\rm max}$ scalings and the slightly larger overall sizes in comparison with the no-delensed and noiseless full-sky results.

In the right panel of Fig. 2, the $r$ dependence of $\Delta f_{\rm NL}$ for the no-delensed and noiseless case (corresponding to the solid lines) and a LiteBIRD-like experiment (corresponding to the dashed lines) is presented. For large $r$, $C_{\ell}^{\rm prim}$ dominates over $C_{\ell}^{BB}$ within a wide range of $\ell$ and hence the $r$ dependence reaches the ideal case (96). On the other hand, for decreasing $r$, $C_{\ell}^{\rm lens}$ or $N_{\ell}$ starts to dominate and errors converge. We find that LiteBIRD could measure an $\mathcal{O}(1)$ signal of $f_{\rm NL}^{ttt,\rm eq}$, $f_{\rm NL}^{ttt,\rm sq}$, or $f_{\rm NL}^{tss,\rm sq}$.555 Comparable detectability is expected in the other next-generation CMB experiments such as CMB-S4 [62] and CORE [63].

If detected, as discussed in Sec. 3, we will have compelling evidence for a deviation from Einstein gravity or the existence of extra source fields, which will help establish an improved picture of the early Universe.

Finally, it is worth mentioning impacts of late-time secondary bispectra. Via gravitational lensing, primordial temperature and $E$-mode polarization fields induce secondary $BTT$ signals [64, 65]. Fortunately, this has a very weak correlation with the $f_{\rm NL}^{tss,\rm sq}$ bispectrum template [16]. Similarly, secondary $BBB$ is also induced, while this is subdominant compared to primordial $BBB$ at interesting $f_{\rm NL}^{ttt,\rm eq}$ or $f_{\rm NL}^{ttt,\rm sq}$ [27]. Besides of these leading-order contributions, higher-order ones exist due to post-Born lensing [66, 67]. Secondary polarized bispectra can also be generated via the Sunyaev-Zel’dovich effect [68] and extragalactic foreground [69, 70]. These impacts are still non-trivial and should therefore be solved for more precise discussions.

6 Conclusions

The search for tensor NG is theoretically well-motivated, and the methodology for testing with the CMB temperature and E/B-mode polarization data has already been established. Observational constraints on some templates from the WMAP temperature and the Planck temperature and $E$-mode polarization data do exist. Any significant signal has not yet been discovered, implying the smallness of tensor NGs. There is still a chance of approaching such a small signal by including $B$-mode polarization in the data analysis. As we have found, an $\mathcal{O}(1)$ signal of $f_{\rm NL}^{ttt,\rm eq}$, $f_{\rm NL}^{ttt,\rm sq}$, or $f_{\rm NL}^{tss,\rm sq}$ would be detectable in LiteBIRD; current constraints are expected to be updated in the next decade.

Throughout this paper, we have focused mainly on scale-invariant cases, while specific scale dependences can be generated for tensor NGs depending on the shapes of the inflationary potentials and non-linear interactions. The detectability is then enhanced at specific multipoles [27] or outside the CMB scales. In the latter case, the information from other probes such as an intereferometric GW survey (e.g., [71, 72]), the galaxy one (e.g., statistical anisotropy in the galaxy power spectrum [73]) and the 21-cm line one also becomes indispensable. A more comprehensive analysis based on multi-wavelength observations remains an interesting and important future issue.

Acknowledgments

I thank the editors of “Status and Prospects of Cosmic Microwave Background Analysis” and especially, Maurizio Tomasi and Michele Liguori for their invitation to contribute to this special issue and for all their support. I also appreciate the careful reviews of the manuscripts by the referees. This work was supported by the JSPS Grant-in-Aid for Research Activity Start-up Grant Number 17H07319 and the JSPS Grant-in-Aid for Early-Career Scientists Grant Number 19K14718. Numerical computations were in part carried out on Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.