of the source position along the polar direction is

$\displaystyle\bar{r}_{z}~{}\delta\theta_{v}$ $\displaystyle=$
$\displaystyle-\bar{r}_{z}\theta_{\alpha}V^{\alpha}_{o}-\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,(\bar{r}_{z}-\bar{r})\theta_{\alpha}(\alpha_{\chi}-\varphi_{\chi})^{,\alpha}+\theta_{\alpha}\delta
x^{\alpha}_{v}~{}$ $\displaystyle=$
$\displaystyle-{\bar{r}_{z}}\theta^{\alpha}(D_{V}{\cal
R}_{,\alpha})\big{|}_{o}-2\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,({\bar{r}_{z}}-\bar{r})D_{\Psi}(\bar{r})\theta^{\alpha}{\cal
R}_{,\alpha}+\int_{0}^{\bar{\eta}_{o}}\mathrm{d}\eta~{}D_{V}(\bar{\eta})\theta^{\alpha}{\cal
R}_{,\alpha}\bigg{|}_{\bm{x}=0}~{},$

where we ignored the orientation $\Omega^{i}$ of the observer frame in the
full expression of $\delta\theta$, as it is not correlated with ${\cal R}$.
Similarly, the angular distortion along the azimuthal direction is identical
with ${\bar{r}_{z}}\sin\theta\delta\phi_{v}$ in the left-hand side and with
$\theta_{\alpha}$ replaced by $\phi_{\alpha}$ in the right-hand side.

Now that we have expressed all the individual components of the matter density
fluctuation in Eq. (4.13) in terms of the initial condition ${\cal R}$, we put
them together here to show our main analytic equation for the matter density
fluctuation in the standard cosmology, and we will use this equation to
compute the bispectrum in the squeezed limit:

$\displaystyle\delta_{m}^{\rm obs}$ $\displaystyle=$ $\displaystyle
D_{1}\left(-\Delta{\cal R}+\frac{3}{2}{\cal R}^{,\alpha}{\cal
R}_{,\alpha}+4{\cal R}\Delta{\cal
R}\right)+\frac{5}{7}D_{A}\nabla_{\alpha}\left({\cal R}^{,\alpha}\Delta{\cal
R}\right)+\frac{1}{7}D_{B}\Delta\left({\cal R}^{,\alpha}{\cal
R}_{,\alpha}\right)$ $\displaystyle+\frac{1}{\mathcal{H}_{z}}D_{V}\Delta{\cal
R}\left[{\cal R}-{\cal R}_{o}-n^{\alpha}\left(D_{V}{\cal
R}_{,\alpha}-D_{V}{\cal
R}_{,\alpha}\big{|}_{o}\right)+2\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,D_{\Psi}(\bar{r})~{}\partial_{r}{\cal
R}(\bar{r}\hat{n})\right]$ $\displaystyle-D_{1}{\partial\over\partial
r}\Delta{\cal R}\bigg{\\{}D_{V}(x_{o}){\cal
R}_{o}-\frac{1}{\mathcal{H}_{z}}\bigg{[}D_{\Psi}{\cal R}+{\cal
R}_{o}+n^{\alpha}\left(D_{V}{\cal R}_{,\alpha}-D_{V}{\cal
R}_{,\alpha}\big{|}_{o}\right)$
$\displaystyle\qquad\qquad-2\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,D_{\Psi}(\bar{r})~{}\partial_{r}{\cal
R}(\bar{r}\hat{n})\bigg{]}+2\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,D_{\Psi}(\bar{r}){\cal
R}(\bar{r}\hat{n})\bigg{\\}}+D_{1}\left(D_{V}\nabla^{\alpha}{\cal
R}\right)_{o}\hat{\nabla}_{\alpha}\left(\Delta{\cal R}\right)$
$\displaystyle+2D_{1}\hat{\nabla}_{\alpha}\left(\Delta{\cal
R}\right)\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,\left(\frac{{\bar{r}_{z}}-\bar{r}}{{\bar{r}_{z}}\bar{r}}\right)D_{\Psi}(\bar{r})\hat{\nabla}^{\alpha}{\cal
R}-D_{1}\nabla^{\alpha}\left(\Delta{\cal
R}\right)\int_{0}^{\bar{\eta}_{o}}\mathrm{d}\eta\,D_{V}\nabla_{\alpha}{\cal
R}\bigg{|}_{\bm{x}=0}~{},$

where $\hat{\nabla}$ is the angular gradient, the first line is the second-
order matter density in Eq. (2.4), and the remaining terms are from $\Delta
x^{\mu}\partial_{\mu}\delta_{v}$ in sequential order. Three components from
the spatial shift of the observer position are combined into one in the last
line, and the product of two angular gradients in the last line represents

$2D_{1}\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,\left(\frac{{\bar{r}_{z}}-\bar{r}}{{\bar{r}_{z}}\bar{r}}\right)D_{\Psi}(\bar{r})\left[\left({\partial\over\partial\theta}\Delta{\cal
R}\right)_{{\bar{r}_{z}}}\left({\partial\over\partial\theta}{\cal
R}\right)_{\bar{r}}+{1\over\sin^{2}\theta}\left({\partial\over\partial\phi}\Delta{\cal
R}\right)_{{\bar{r}_{z}}}\left({\partial\over\partial\phi}{\cal
R}\right)_{\bar{r}}\right]~{}.$ (A.29)

Equation (A.1) is the main equation shown in Eq. (4.2) for our computation.

Further simplification can be made in the limit that the Universe is
approximated as the Einstein-de Sitter universe, or matter-dominated universe,
in which $\Lambda=0$ and the matter density parameter $\Omega_{m}=1$ ($f=1$).
Due to the simplicity, we can derive the analytic solutions, regarding the
Hubble parameter and the angular diameter distance

$\displaystyle a$ $\displaystyle=$ $\displaystyle\left({t\over
t_{o}}\right)^{2/3}=\left({\eta\over\eta_{o}}\right)^{2}~{},\qquad\qquad{t\over
t_{o}}=\left({\eta\over\eta_{o}}\right)^{3}~{},\qquad\qquad\eta_{o}=3t_{o}~{},$
(A.30) $\displaystyle H$ $\displaystyle=$ $\displaystyle{2\over
3t}~{},\qquad\qquad\mathcal{H}={2\over\eta}~{},\qquad\qquad\bar{r}_{z}=\bar{\eta}_{o}-\bar{\eta}={2\over
H_{0}}\left(1-{1\over\sqrt{1+z}}\right)~{},$ (A.31)

and the relation among the perturbation variables

$\Sigma=\frac{5}{2}\,,\qquad\qquad D_{\Psi}=-\frac{3}{5}\,,\qquad\qquad
D_{1}=\frac{\eta^{2}}{10}\,,\qquad\qquad D_{V}=\frac{\eta}{5}\,,\qquad\qquad
D_{A}=D_{B}=D_{1}^{2}={\eta^{4}\over 100}~{}.$ (A.32)

###  A.2 Fourier decomposition and one-point ensemble average

To facilitate the subsequent calculations, we introduce a set of second-order
Fourier kernels $F(\mathbf{k}_{1},\mathbf{k}_{2})$ that capture various
contributions to the matter density fluctuation in Eq. (A.1) at the observed
redshift and describe their time evolution from the initial condition ${\cal
R}$. Schematically, the matter density fluctuation at the second order will be
expressed as

$\delta(\mathbf{x})\propto\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}~{}e^{i\mathbf{k}_{1}\cdot\mathbf{x}}e^{i\mathbf{k}_{2}\cdot\mathbf{x}}F(\mathbf{k}_{1},\mathbf{k}_{2};{\bf\hat{n}},\bar{\eta}_{z}){\cal
R}(\mathbf{k}_{1}){\cal R}(\mathbf{k}_{2})\,,$ (A.33)

where the Fourier kernel is dimensionless. As shown in Eqs. (4.13) and (A.1),
the observed matter density fluctuation at the second order is coupled with
the contributions at the source position, at the observer position, or along
the line-of-sight direction. So we discuss three different types of Fourier
kernels in turn. The derivation of Fourier kernels for the individual
contributions is presented in Section B.

Furthermore, it proves useful for the computational convenience and also
conceptual clarity to define Fourier counter parts, given the individual
contributions to the matter density fluctuation. These Fourier counter parts
in configuration space and Fourier space are defined in a hypersurface set by
the observed redshift $z$ in the usual way:

$\delta(\mathbf{x})=\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\,e^{i\mathbf{k}\cdot\mathbf{x}}~{}\delta(\mathbf{k};{\bf\hat{n}},\bar{\eta}_{z})\,,\qquad\qquad\delta(\mathbf{k};{\bf\hat{n}},\bar{\eta}_{z})=\int
d^{3}x\,e^{-i\mathbf{k}\cdot\mathbf{x}}~{}\delta(\mathbf{x})~{},$ (A.34)

where the volume integral over the position $\mathbf{x}$ is all over the
infinite hypersurface, not over the light cone volume. Note that this
hypersurface encompasses a volume outside the observed light cone volume and
hence these Fourier counter parts are not directly observable, except at the
intersection with the light cone volume. However, we can imagine that
fictitious observers at different spatial position in the Universe perform the
same observations. In other words, for the same observed redshift and angle,
these fictitious observers can construct the fluctuation field
$\delta(\mathbf{x})$ with the observed position
$\mathbf{x}=\bar{r}_{z}{\bf\hat{n}}+\mathbf{x}_{o}$, where $\mathbf{x}_{o}$ is
the spatial position of the fictitious observer and can be set zero for the
real observer (us). Noting that our position in the Universe is not special,
it is conceptually useful to think of such fictitious observations and to
construct the observed matter density field and its Fourier counter part
outside our own light cone volume. However, it is noted that since the
observed matter density fluctuation depends not only on the redshift, but also
on the angle, the observed angle should be specified, when the fictitious
observations are considered, and hence the angular dependence in Eq. (A.34).
We refer the reader to the work [65] for more detailed discussion and
computation.

The Fourier kernels are also useful in computing the one-point ensemble
average. The ensemble average $\langle\delta(\mathbf{x})\rangle$ is an average
of the fluctuation $\delta(\mathbf{x})$ over many realizations of the
Universe, as discussed in Section 3.3. With the ergodic theorem, the ensemble
average is equivalent to the Euclidean average over the hypersurface, which
can be readily computed in our formalism with the observations by the
fictitious observers. Note that the cosmic variance in practice arises due to
our limitation to one observer position or the lack of average over
translation in the Euclidean average (see, e.g., [83]). The average over the
hypersurface yields

$\langle\delta(\mathbf{x})\rangle\underset{{\rm
Erg.}}{=}\lim_{V\rightarrow\infty}{1\over V}\int
d^{3}x\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\,e^{i\mathbf{k}\cdot\mathbf{x}}~{}\delta(\mathbf{k};{\bf\hat{n}},\bar{\eta}_{z})={1\over
V}\big{\langle}\delta(\mathbf{k}\equiv
0;{\bf\hat{n}},\bar{\eta}_{z})\big{\rangle}~{}.$ (A.35)

In general, the ensemble average of a fluctuation is zero for a Gaussian
distribution. However, at the second order, and in particular with the
relativistic contributions, the ensemble average of the matter density
fluctuation is non-vanishing, as we show in the following. More importantly,
it is apparent by the definition of the average over a hypersurface that the
ensemble average depends on a choice of hypersurface.

$\bullet$ Coupling terms with contributions at the source position.— First, we
consider the contributions $\delta_{s}({\bf x})$ at the source position only,
and these terms in Eq. (A.1) can be expressed as

$\delta_{s}(\mathbf{x})=:\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}~{}e^{i\mathbf{k}_{1}\cdot\mathbf{x}}e^{i\mathbf{k}_{2}\cdot\mathbf{x}}F_{s}(\mathbf{k}_{1},\mathbf{k}_{2};{\bf\hat{n}},\bar{\eta}_{z}){\cal
R}(\mathbf{k}_{1}){\cal R}(\mathbf{k}_{2})\,,$ (A.36)

and this yields the Fourier mode according to Eq. (A.34)

$\delta_{s}(\mathbf{k};{\bf\hat{n}},\bar{\eta}_{z})=\int\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}~{}F_{s}(\mathbf{q},\mathbf{k}-\mathbf{q};{\bf\hat{n}},\bar{\eta}_{z})~{}{\cal
R}(\mathbf{q}){\cal R}(\mathbf{k}-\mathbf{q})\,,$ (A.37)

where the Fourier kernel $F_{s}$ is subject to

$F^{\ast}_{s}(\mathbf{k}_{1},\mathbf{k}_{2};{\bf\hat{n}},\bar{\eta}_{z})=F^{\ast}_{s}(-\mathbf{k}_{1},-\mathbf{k}_{2};{\bf\hat{n}},\bar{\eta}_{z})\,.$
(A.38)

Note that the Fourier kernels include the time-dependence. For example, the
linear-order matter density contribution to the Fourier kernel is

$\delta^{(1)}_{m}(\mathbf{x},t)=-D_{1}(t)\Delta{\cal
R}(\mathbf{x})~{},\qquad\qquad F_{s}\ni~{}D_{1}(t)k^{2}~{}.$ (A.39)

For later convenience, we introduce another notation for Fourier kernels
${\cal F}$ that are directly related to the bispectrum, and for the
contributions $\delta_{s}({\bf x})$ at the source position we simply have

${\cal
F}_{s}(\mathbf{k}_{1},\mathbf{k}_{2};{\bf\hat{n}},\bar{\eta}_{z}):=F_{s}(\mathbf{k}_{1},\mathbf{k}_{2},{\bf\hat{n}},\bar{\eta}_{z})~{}.$
(A.40)

Finally, the ensemble average of the coupling terms with contributions at the
source position is obtained by using Eq. (A.35) as

$\Xi:=\langle\delta_{s}(\mathbf{x})\rangle=\int\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}F_{s}(\mathbf{q},-\mathbf{q};{\bf\hat{n}},\bar{\eta}_{z})P_{\cal
R}(q)=\int d\ln q~{}\Delta_{\cal R}^{2}(q)\int{d\mu_{q}\over
2}~{}F_{s}(\mathbf{q},-\mathbf{q};{\bf\hat{n}},\bar{\eta}_{z})~{},$ (A.41)

where $\mu_{q}=\mathbf{q}\cdot{\bf\hat{n}}/|\mathbf{q}|$ is the cosine angle
with respect to the observed direction.

$\bullet$ Coupling terms with contributions at the observer position.— Next we
consider the contributions $\delta_{o}({\bf x})$ to the matter density
fluctuation coupled with those at the observer position. To compute the field
in a hypersurface outside the light cone, we consider fictitious observers at
different spatial positions ($\mathbf{x}_{o}\neq 0$). Those coupling terms in
Eq. (A.1) can be expressed in terms of Fourier kernel $F_{o}$ as

$\delta_{o}(\mathbf{x})=:\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}~{}e^{i\mathbf{k}_{1}\cdot\mathbf{x}_{o}}e^{i\mathbf{k}_{2}\cdot\mathbf{x}}F_{o}(\mathbf{k}_{1},\mathbf{k}_{2};{\bf\hat{n}},\bar{\eta}_{o},\bar{\eta}_{z}){\cal
R}(\mathbf{k}_{1}){\cal R}(\mathbf{k}_{2})\,,$ (A.42)

where the Fourier kernel also has the time-dependence on $\bar{\eta}_{o}$ set
at the observer position ${\bf x}_{o}$. Note that since $\bar{r}_{z}$ and
${\bf\hat{n}}$ are fixed in a hypersurface, the position ${\bf x}_{o}$ for the
fictitious observer is a function of $\mathbf{x}$ in consideration. Also mind
that the Fourier kernel $F_{o}$ is not symmetric over the arguments
$\mathbf{k}_{1}$ and $\mathbf{k}_{2}$ in our convention, where the wave vector
$\mathbf{k}_{1}$ belongs to the Fourier transform of the contribution
evaluated at the observer position and $\mathbf{k}_{2}$ belongs to the
contribution at the source position. The Fourier counter part is then

$\delta_{o}(\mathbf{k};{\bf\hat{n}},\bar{\eta}_{o},\bar{\eta}_{z})=\int\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}~{}e^{-i\mathbf{q}\cdot\bar{r}_{z}{\bf\hat{n}}}~{}F_{o}(\mathbf{q},\mathbf{k}-\mathbf{q};{\bf\hat{n}},\bar{\eta}_{o},\bar{\eta}_{z})~{}{\cal
R}(\mathbf{q}){\cal R}(\mathbf{k}-\mathbf{q})\,,$ (A.43)

and its ensemble average is

$\displaystyle\Xi$ $\displaystyle=$
$\displaystyle\langle\delta_{o}(\mathbf{x})\rangle=\int\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}~{}e^{-i{\bar{r}_{z}}\mathbf{q}\cdot{\bf\hat{n}}}F_{o}(\mathbf{q},-\mathbf{q};{\bf\hat{n}},\bar{\eta}_{o},\bar{\eta}_{z})P_{\cal
R}(q)$ $\displaystyle=$ $\displaystyle\int d\ln q~{}\Delta_{\cal
R}^{2}(q)\int{d\mu_{q}\over
2}~{}e^{-iq\bar{r}_{z}\mu_{q}}F_{o}(\mathbf{q},-\mathbf{q};{\bf\hat{n}},\bar{\eta}_{o},\bar{\eta}_{z})~{}.$

We define the Fourier kernel

${\cal
F}_{o}(\mathbf{k}_{1},\mathbf{k}_{2};{\bf\hat{n}},\bar{\eta}_{o},\bar{\eta}_{z}):=e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}F_{o}(\mathbf{k}_{1},\mathbf{k}_{2};{\bf\hat{n}},\bar{\eta}_{o},\bar{\eta}_{z})~{},$
(A.45)

and note that it includes exponential factor.

$\bullet$ Coupling terms with contributions along the line-of-sight
direction.— Finally, we consider the contributions $\delta_{nl}({\bf x})$ to
the matter density fluctuation coupled with those along the line-of-sight
direction:

$\delta_{nl}(\mathbf{x})=:\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,~{}e^{i\bar{r}\mathbf{k}_{1}\cdot{\bf\hat{n}}}e^{i{\bar{r}_{z}}\mathbf{k}_{2}\cdot{\bf\hat{n}}}~{}F_{nl}(\mathbf{k}_{1},\mathbf{k}_{2};{\bf\hat{n}},\bar{r},\bar{\eta}_{z}){\cal
R}(\mathbf{k}_{1}){\cal R}(\mathbf{k}_{2})\,,$ (A.46)

where the time-dependence of the Fourier kernel is specified by $\bar{r}$ in
the line-of-sight integration. We use the subscript $nl$ to refer to this as
non-local Fourier kernel. The kernel is also not symmetric in the arguments,
and the wave vector $\mathbf{k}_{1}$ belongs to the contribution evaluated
along the line-of-sight, while $\mathbf{k}_{2}$ to the one evaluated at the
source position. The Fourier counter part is then

$\delta_{nl}(\mathbf{k};{\bf\hat{n}},\bar{\eta}_{z})=\int\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,e^{-i\Delta
r\mathbf{q}\cdot{\bf\hat{n}}}~{}F_{nl}(\mathbf{q},\mathbf{k}-\mathbf{q};\bar{r},\bar{\eta}_{z})~{}{\cal
R}(\mathbf{q}){\cal R}(\mathbf{k}-\mathbf{q})\,,$ (A.47)

and its ensemble average is

$\displaystyle\Xi$ $\displaystyle=$
$\displaystyle\langle\delta_{nl}(\mathbf{x})\rangle=\int\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,e^{-i\Delta
r\mathbf{q}\cdot{\bf\hat{n}}}~{}F_{nl}(\mathbf{q},-\mathbf{q};{\bf\hat{n}},\bar{r},\bar{\eta}_{z})P_{\cal
R}(q)$ $\displaystyle=$ $\displaystyle\int d\ln q~{}\Delta_{\cal
R}^{2}(q)\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,\int{d\mu_{q}\over
2}~{}e^{-iq\Delta
r\mu_{q}}F_{nl}(\mathbf{q},-\mathbf{q};{\bf\hat{n}},\bar{r},\bar{\eta}_{z})~{}.$

where $\Delta r:=\bar{r}_{z}-\bar{r}$. The Fourier kernel for the non-local
contributions is then defined as

${\cal
F}_{nl}(\mathbf{k}_{1},\mathbf{k}_{2};{\bf\hat{n}},\bar{\eta}_{z}):=\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,~{}e^{-i\mathbf{k}_{1}\cdot\Delta
r{\bf\hat{n}}}F_{nl}(\mathbf{k}_{1},\mathbf{k}_{2};{\bf\hat{n}},\bar{r},\bar{\eta}_{z})~{},$
(A.49)

and it includes the line-of-sight integration and the exponential factor in
the kernel. While $F_{nl}$ is dimensionful, the Fourier kernel ${\cal F}_{nl}$
is dimensionless due to the line-of-sight integral.

###  A.3 Three-point correlation and the bispectrum

Here we compute the three-point correlation function and its bispectrum. The
leading-order contribution to the three-point statistics arises from the
contraction of two linear-order contributions and a second-order contribution
of the matter density fluctuation. Since the Fourier mode of the observed
matter density fluctuation at the second order in perturbations is
collectively expressed in terms of Fourier kernels as

$\delta(\mathbf{k};{\bf\hat{n}},\bar{\eta}_{z})=\int\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}~{}{\cal
F}(\mathbf{q},\mathbf{k}-\mathbf{q};{\bf\hat{n}},\bar{\eta}_{z})~{}{\cal
R}(\mathbf{q}){\cal R}(\mathbf{k}-\mathbf{q})\,,$ (A.50)

and the linear-order matter density fluctuation is

$\delta^{(1)}(\mathbf{k})=D_{1}(z)k^{2}{\cal R}(k)~{},$ (A.51)

the leading-order contribution to the bispectrum is then

$\displaystyle\Big{\langle}\delta^{(1)}(\mathbf{k}_{1})\delta^{(1)}(\mathbf{k}_{2})\delta^{(2)}(\mathbf{k}_{3})\Big{\rangle}$
$\displaystyle=$ $\displaystyle D_{1}^{2}k_{1}^{2}k_{2}^{2}P_{\cal
R}(k_{1})P_{\cal
R}(k_{2})(2\pi)^{3}\delta^{D}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})$
$\displaystyle\times\left[{\cal
F}(-\mathbf{k}_{1},-\mathbf{k}_{2};{\bf\hat{n}}_{3})+{\cal
F}(-\mathbf{k}_{2},-\mathbf{k}_{1};{\bf\hat{n}}_{3})\right]+\Big{\langle}\delta^{(1)}(\mathbf{k}_{1})\delta^{(1)}(\mathbf{k}_{2})\Big{\rangle}\Big{\langle}\delta^{(2)}(\mathbf{k}_{3})\Big{\rangle}~{},$

where we omitted the time dependence in the growth function $D_{1}$ and the
Fourier kernels ${\cal F}$. Note that the tadpole term exists because the one-
point ensemble average $\left\langle\delta^{(2)}(\mathbf{k}_{3})\right\rangle$
is not vanishing. Accounting for the permutation of three leading-order
contributions and for three different observed positions at the same observed
redshift, we derive the observed three-point correlation function

$\displaystyle\Big{\langle}\delta(\mathbf{x}_{1})\delta(\mathbf{x}_{2})\delta(\mathbf{x}_{3})\Big{\rangle}=\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}k_{3}}{(2\pi)^{3}}e^{i\mathbf{k}_{1}\cdot\mathbf{x}_{1}}e^{i\mathbf{k}_{2}\cdot\mathbf{x}_{2}}e^{i\mathbf{k}_{3}\cdot\mathbf{x}_{3}}\Big{\langle}\delta(\mathbf{k}_{1})\delta(\mathbf{k}_{2})\delta(\mathbf{k}_{3})\Big{\rangle}$
(A.53)
$\displaystyle\qquad\qquad=\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}e^{i\mathbf{k}_{1}\cdot(\mathbf{x}_{1}-\mathbf{x}_{3})}e^{i\mathbf{k}_{2}\cdot(\mathbf{x}_{2}-\mathbf{x}_{3})}\bigg{\\{}D_{1}^{2}k_{1}^{2}k_{2}^{2}P_{\cal
R}(k_{1})P_{\cal R}(k_{2})\bigg{[}{\cal
F}(-\mathbf{k}_{1},-\mathbf{k}_{2};{\bf\hat{n}}_{3})~{}~{}~{}~{}$
$\displaystyle\qquad\qquad\qquad\qquad+{\cal
F}(-\mathbf{k}_{2},-\mathbf{k}_{1};{\bf\hat{n}}_{3})\bigg{]}+{\rm
perm.}\bigg{\\}}+\Big{[}\xi_{m}(\mathbf{x}_{1},\mathbf{x}_{2})+\xi_{m}(\mathbf{x}_{2},\mathbf{x}_{3})+\xi_{m}(\mathbf{x}_{3},\mathbf{x}_{1})\Big{]}\Xi~{},$

where the permutation terms contain the Fourier kernels ${\cal F}$ with two
other angular positions ${\bf\hat{n}}_{1}$ and ${\bf\hat{n}}_{2}$ and the two-
point correlation function is the just the linear-order matter correlation
function

$\xi_{m}(\mathbf{x}_{1},\mathbf{x}_{2})=\left\langle\delta^{(1)}(\mathbf{x}_{1})\delta^{(1)}(\mathbf{x}_{2})\right\rangle~{}.$
(A.54)

The three-point correlation is non-vanishing, even for a configuration
$\mathbf{x}_{1}=\mathbf{x}_{2}$ and $\mathbf{x}_{3}\rightarrow\infty$ due to
the non-vanishing constant contribution
$\Xi=\langle\delta(\mathbf{x})\rangle\neq 0$, or the difference between the
ensemble average and the background at a given time.

To address this subtlety, we first define a dimensionless fluctuation

${\mathfrak{D}}(\mathbf{x}):=\delta(\mathbf{x})-\Xi~{},\qquad\qquad\Big{\langle}{\mathfrak{D}}(\mathbf{x})\Big{\rangle}=0~{},$
(A.55)

though the correct fluctuation should be with extra constant $(1+\Xi)$. For
instance, the matter density can be split as

$\rho(\mathbf{x})=\bar{\rho}(t)\left(1+\Xi\right)\left(1+{{\mathfrak{D}}(\mathbf{x})\over
1+\Xi}\right)~{},\qquad\qquad\left\langle\rho({\bf
x})\right\rangle=\bar{\rho}(t)(1+\Xi)~{}.$ (A.56)

The three-point correlation function of ${\mathfrak{D}}(\mathbf{x})$ is then

$\Big{\langle}{\mathfrak{D}}(\mathbf{x}_{1}){\mathfrak{D}}(\mathbf{x}_{2}){\mathfrak{D}}(\mathbf{x}_{3})\Big{\rangle}=\Big{\langle}\delta(\mathbf{x}_{1})\delta(\mathbf{x}_{2})\delta(\mathbf{x}_{3})\Big{\rangle}-\bigg{[}\Big{\langle}\delta(\mathbf{x}_{1})\delta(\mathbf{x}_{2})\Big{\rangle}~{}\Big{\langle}\delta(\mathbf{x}_{3})\Big{\rangle}+{\rm
perm.}\bigg{]}+2\big{\langle}\delta\big{\rangle}^{3}~{},$ (A.57)

where the extra terms in the three-point correlation including the constant
contribution ensure that the tadpole terms do not contribute and only the
connected contribution to the three-point correlation function remains.
Computing to the second-order in the matter density fluctuation, we derive the
three-point correlation function

$\displaystyle\Big{\langle}{\mathfrak{D}}(\mathbf{x}_{1}){\mathfrak{D}}(\mathbf{x}_{2}){\mathfrak{D}}(\mathbf{x}_{3})\Big{\rangle}$
$\displaystyle=$ $\displaystyle
2~{}\Xi^{3}+\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}~{}e^{i\mathbf{k}_{1}\cdot(\mathbf{x}_{1}-\mathbf{x}_{3})}e^{i\mathbf{k}_{2}\cdot(\mathbf{x}_{2}-\mathbf{x}_{3})}$
$\displaystyle\times\left\\{D_{1}^{2}k_{1}^{2}k_{2}^{2}P_{\cal
R}(k_{1})P_{\cal R}(k_{2})\bigg{[}{\cal
F}(-\mathbf{k}_{1},-\mathbf{k}_{2};{\bf\hat{n}}_{3})+{\cal
F}(-\mathbf{k}_{2},-\mathbf{k}_{1};{\bf\hat{n}}_{3})\bigg{]}+{\rm
perm.}\right\\}~{},$

where the tadpoles in Eq. (A.53) are all removed. Furthermore, since we only
considered the field up to the second order in perturbations or keep terms up
to $P_{\cal R}^{2}$, the constant term $\Xi^{3}\sim P_{\cal R}^{3}$ will be
ignored for consistency, as a proper treatment of such term would need one-
loop contributions $P_{\cal R}^{3}$ in the bispectrum.

Having derived the general expression for the three-point statistics, we now
consider a special triangular configuration, in which two observed positions
are identical and the third position points in the opposite side of the
sky555While these two identical points are then subject to large non-linear
corrections, these two points for the squeezed triangular configuration only
need to be close to each other, compared to the third position, avoiding any
extra complication due to nonlinearity on small scales.

$z:=z_{1}=z_{2}=z_{3}~{},\qquad\qquad{\bf\hat{n}}:={\bf\hat{n}}_{1}={\bf\hat{n}}_{2}=-{\bf\hat{n}}_{3}~{},$
(A.59)

i.e., one observed redshift $z$ and one angular vector ${\bf\hat{n}}$ for the
three points in the sky. In terms of three-dimensional position vectors, the
squeezed triangular configuration is represented by

$\mathbf{x}_{2}=\mathbf{x}_{1}~{},\qquad\qquad\mathbf{L}:=\mathbf{x}_{1}-\mathbf{x}_{3}~{},$
(A.60)

and our primary interest lies in the squeezed triangular configuration in the
limit $L\rightarrow\infty$. The squeezed correlation function can then be
further simplified as

$\xi_{\rm
sqz}:=\Big{\langle}{\mathfrak{D}}(\mathbf{x}_{1}){\mathfrak{D}}(\mathbf{x}_{2}){\mathfrak{D}}(\mathbf{x}_{3})\Big{\rangle}=\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}~{}e^{i(\mathbf{k}_{1}+\mathbf{k}_{2})\cdot\mathbf{L}}\bigg{[}B_{112}+B_{211}+B_{121}\bigg{]}(\mathbf{k}_{1},\mathbf{k}_{2},-\mathbf{k}_{12})~{},$
(A.61)

where $\mathbf{k}_{12}:=\mathbf{k}_{1}+\mathbf{k}_{2}$. The connected
bispectra with dimension $L^{6}$ in the integrand are defined as

$\displaystyle B_{112}$ $\displaystyle:=$ $\displaystyle
D_{1}^{2}k_{1}^{2}k_{2}^{2}P_{\cal R}(k_{1})P_{\cal R}(k_{2})\bigg{[}{\cal
F}(-\mathbf{k}_{1},-\mathbf{k}_{2};-{\bf\hat{n}})+{\cal
F}(-\mathbf{k}_{2},-\mathbf{k}_{1};-{\bf\hat{n}})\bigg{]}~{},$ (A.62)
$\displaystyle B_{211}$ $\displaystyle:=$ $\displaystyle
D_{1}^{2}k_{2}^{2}k_{3}^{2}P_{\cal R}(k_{2})P_{\cal R}(k_{3})\bigg{[}{\cal
F}(-\mathbf{k}_{2},-\mathbf{k}_{3};{\bf\hat{n}})+{\cal
F}(-\mathbf{k}_{3},-\mathbf{k}_{2};{\bf\hat{n}})\bigg{]}~{},$ (A.63)
$\displaystyle B_{121}$ $\displaystyle:=$ $\displaystyle
D_{1}^{2}k_{3}^{2}k_{1}^{2}P_{\cal R}(k_{3})P_{\cal R}(k_{1})\bigg{[}{\cal
F}(-\mathbf{k}_{3},-\mathbf{k}_{1};{\bf\hat{n}})+{\cal
F}(-\mathbf{k}_{1},-\mathbf{k}_{3};{\bf\hat{n}})\bigg{]}~{},$ (A.64)

where the subscripts represent the order of perturbations, e.g.,

$B_{112}:=\left\langle{\mathfrak{D}}^{(1)}_{{\bf
k}_{1}}{\mathfrak{D}}^{(1)}_{{\bf k}_{2}}{\mathfrak{D}}^{(2)}_{{\bf
k}_{3}}\right\rangle~{},$ (A.65)

and the dimension of the bispectra is $L^{6}$. The task for computing the
squeezed correlation boils down to computing the Fourier kernels for Eqs.
(A.62)$-$(A.64), which are presented in Section B. Despite the simplicity in
the notation, it is noted that the kernels ${\cal F}_{o}$ and ${\cal F}_{nl}$
for the observer and the non-local contributions include the exponential
factors or the line-of-sight integral. Moreover, though those two types of
Fourier kernels are not symmetric over the arguments, their contributions to
the squeezed correlation function are made symmetric in the bispectra.
However, mind the dependence of ${\bf\hat{n}}$ and its sign.

Treating the squeezed three-point correlation function as a two-point
correlation or just a function of separation ${\bf L}$, the Fourier counter
part or the power spectrum can be defined as

$P(k):=\int d^{3}L~{}e^{-i\mathbf{k}\cdot\mathbf{L}}~{}\xi_{\rm
sqz}\simeq\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\bigg{[}B_{112}+B_{211}+B_{121}\bigg{]}(\mathbf{k}_{1},\mathbf{k}-\mathbf{k}_{1},-\mathbf{k})~{},$
(A.66)

where $\mathbf{k}=\mathbf{k}_{1}+\mathbf{k}_{2}$. It is shown [60] that this
power spectrum is a useful quantity to measure the correction to the standard
galaxy power spectrum, arising from the non-Gaussian nature of the matter
density fluctuation, or the squeezed three-point correlation $\xi_{\rm sqz}$.
However, it is important to note that $\xi_{\rm sqz}$ here is the observed
three-point correlation function, which is not just a function of its
separation $\bm{L}$ alone. Furthermore, $\xi_{\rm sqz}$ is defined on the past
light cone, not in an infinite hypersurface of simultaneity at the observed
redshift, so that any integration over $\bm{L}$ should involves the variation
in time or the observed redshift. Therefore, the final equality for the power
spectrum in terms of three bispectrum should be taken as a simple theoretical
measure of non-Gaussianity, rather than a real observed power spectrum that
can be obtained in the light cone volume. While one can still define the
observed power spectrum as in Eq. (A.66) with the volume integration over the
light-cone volume, it is not equal to the expression in the RHS of Eq. (A.66)
in terms of the bispectra.

Before we proceed to derive the individual Fourier kernels, we discuss the
general scaling of the squeezed correlation function (see also [148]). Though
we have computed the three-point correlation function $\xi_{\rm sqz}$, this
squeezed correlation function contributes to the galaxy two-point correlation
function as a non-Gaussian correction discussed in Section 2.2. For such two-
point correlation function, in the absence of non-vanishing ensemble average,
we expect it to vanish in the limit $L\rightarrow\infty$. Assuming a power-law
relation, we expect the power-law index to be at least positive:

$\xi_{\rm sqz}\propto L^{-n}~{},\qquad\qquad\lim_{L\rightarrow\infty}\xi_{\rm
sqz}=0~{},\qquad\qquad n>0~{}.$ (A.67)

This expectation allows the possibility for the power spectrum $P(k)\propto
k^{n-3}$ to scale as $P(k)\propto P_{m}/k^{2}\propto P_{\delta\phi}$ or the
primordial non-Gaussianity signature discussed in Section 2.2.

In the limit $L\rightarrow\infty$, the squeezed triangular configuration takes
the form

$\mathbf{k}_{l}:=\mathbf{k}_{12}=\mathbf{k}_{1}+\mathbf{k}_{2}\rightarrow
0~{},\qquad\qquad\mathbf{k}_{1}\approx-\mathbf{k}_{2}~{},\qquad\qquad\mathbf{k}_{3}=-\mathbf{k}_{12}=-\mathbf{k}_{l}~{},$
(A.68)

and the connected bispectra in this limit are

$\displaystyle B_{112}$ $\displaystyle=$ $\displaystyle
D_{1}^{2}k_{1}^{4}P_{\cal R}(k_{1})P_{\cal R}(k_{1})\bigg{[}{\cal
F}(-\mathbf{k}_{1},\mathbf{k}_{1};-{\bf\hat{n}})+{\cal
F}(\mathbf{k}_{1},-\mathbf{k}_{1};-{\bf\hat{n}})\bigg{]}~{},$ (A.69)
$\displaystyle B_{211}$ $\displaystyle=$ $\displaystyle
D_{1}^{2}k_{1}^{2}k_{l}^{2}P_{\cal R}(k_{1})P_{\cal R}(k_{l})\bigg{[}{\cal
F}(\mathbf{k}_{1},\mathbf{k}_{l};{\bf\hat{n}})+{\cal
F}(\mathbf{k}_{l},\mathbf{k}_{1};{\bf\hat{n}})\bigg{]}~{},$ (A.70)
$\displaystyle B_{121}$ $\displaystyle=$ $\displaystyle
D_{1}^{2}k_{l}^{2}k_{1}^{2}P_{\cal R}(k_{l})P_{\cal R}(k_{1})\bigg{[}{\cal
F}(\mathbf{k}_{l},-\mathbf{k}_{1};{\bf\hat{n}})+{\cal
F}(-\mathbf{k}_{1},\mathbf{k}_{l};{\bf\hat{n}})\bigg{]}~{}.$ (A.71)

While $k_{l}\rightarrow 0$, the other wave vector $\mathbf{k}_{1}$ can take
any arbitrary value. As the exponential factor in Eq. (A.61) vanishes, the
squeezed correlation function is determined by the integration of these
connected bispectra. Note that the contribution of $B_{112}$ is a constant in
this limit or independent of $k_{l}$, and its ratio to the other contributions
is

${B_{112}\over B_{211}}\sim{B_{112}\over B_{121}}\sim{k_{1}^{2}P_{\cal
R}(k_{1})\over k_{l}^{2}P_{\cal R}(k_{l})}\sim\left({k_{l}\over
k_{1}}\right)^{2-n_{s}}~{},$ (A.72)

suppressed at least by $k_{l}/k_{1}$ before considering extra suppression
factors from the ratio of the Fourier kernels. The calculations in Appendix B
show that the ratio of two Fourier kernels vanish due to the symmetry of the
kernels in $B_{112}$ or scales with some power of $k_{l}$. Therefore, we will
ignore the contribution from $B_{112}$. The long-mode contributions in
$B_{211}$ and $B_{121}$ to the power spectrum are schematically

$\Delta P(k_{l})\propto\int d\ln k_{1}~{}k_{1}P_{m}(k_{1})P_{m}(k_{l})\int
d\mu_{1}~{}{[{\cal F}+{\cal F}]\over k_{l}^{2}}~{},\qquad\qquad
P_{m}(k)\propto k^{4}P_{\cal R}(k)~{},$ (A.73)

suggesting that the scaling with $k_{l}$ for the sum of two kernels in
$B_{211}$ or $B_{121}$ after the angular integration should be

$\Delta P\propto k_{l}^{n_{s}-2+n_{{\cal F}}}~{},\qquad\qquad n_{{\cal
F}}>-1-n_{s}\simeq-2~{},$ (A.74)

to satisfy the condition in Eq. (A.67), where we approximated

$\int d\mu\bigg{[}{\cal F}+{\cal F}\bigg{]}\propto k_{l}^{n_{{\cal F}}}~{}.$
(A.75)

As demonstrated in Appendix B, the power-law slope is non-negative: $n_{\cal
F}\geq 0$.

##  Appendix B Detailed calculations of Fourier kernels

Here we present the detailed calculations of the Fourier kernels for three
different types of contributions to the observed matter density fluctuation.
With the kernels in hands, the one-point ensemble average and three connected
bispectra can be readily computed. Each subsection specifies the calculations
according to three different types of Fourier kernels.

The squeezed triangular configuration, which is our primary interest, is
specified by only two angular directions ${\bf\hat{n}}_{1}$ and
${\bf\hat{n}}_{2}=-{\bf\hat{n}}_{1}$, opposite directions in the sky. To
facilitate the calculation, we set this direction to be aligned with
$z$-direction
(${\bf\hat{n}}={\bf\hat{n}}_{1}={\bf\hat{z}}=-{\bf\hat{n}}_{2}$), in which the
azimuthal integration is trivial. We will use the following relations for the
polar integration in terms of spherical Bessel functions $j_{n}(x)$

$\displaystyle\int\mathrm{d\mu\over 2}\cos(\mu
a)=j_{0}(a)\,,\qquad\qquad\int{\mathrm{d}\mu\over 2}\mu\sin(\mu
a)=j_{1}(a)\,,$ (B.1) $\displaystyle\int{\mathrm{d}\mu\over 2}\mu^{2}\cos(\mu
a)={j_{0}(a)-2j_{2}(a)\over 3}~{}.$

###  B.1 Contributions at the source position

The contributions at the source position arise from three different origins.
The first is the second-order matter density contributions in the standard
perturbation theory shown Eq. (2.1), and the second is the second-order
relativistic contribution to the matter density fluctuation in Eq. (2.4).
Finally, the last arises from the coupling terms of the linear-order matter
density fluctuation and the relativistic contribution at the source position
in Eq. (4.13). We present the detailed calculations according to this
classification.

Contributions from the standard perturbation theory: The second-order
contributions in the standard perturbation theory consist of two terms in Eq.
(2.1). The Fourier kernel for the first term is

$\bullet~{}~{}\frac{5}{7}D_{A}(z)\nabla_{\alpha}\left({\cal
R}^{\alpha}\Delta{\cal R}\right)~{}:\qquad\qquad
F_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=\frac{5}{7}D_{A}(z)\left(\mathbf{k}_{1}+\mathbf{k}_{2}\right)\cdot\mathbf{k}_{1}~{}k_{2}^{2}~{},$
(B.2)

and given the structure of the connected bispectra we symmetrize the Fourier
kernel

${\cal
F}_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=\frac{5}{14}D_{A}(z)\left(\mathbf{k}_{1}+\mathbf{k}_{2}\right)\cdot\left(k_{2}^{2}\mathbf{k}_{1}+k_{1}^{2}\mathbf{k}_{2}\right)~{}.$
(B.3)

Having derived the Fourier kernel, the one-point ensemble average in Eq.
(A.41) is readily shown to be vanishing

$\Xi=0~{},$ (B.4)

as the sum of two wave vectors is zero. The contribution to the connected
bispectra can be derived by using Eqs. (A.69)$-$(A.71). The first connected
bispectra $B_{112}$ in Eq. (A.69) is independent of $k_{l}$ and is simply zero

$B_{112}=0$ (B.5)

due to symmetry of the wave vectors. Since the pre-factors of three bispectra
are the same for different Fourier kernels, we present calculations of the sum
of Fourier kernels in the square bracket in Eqs. (A.69)$-$(A.71). The other
two connected bispectra are non-zero,

$\displaystyle
B_{211}\ni\frac{5}{7}D_{A}(z)\left[2k_{l}^{2}k_{1}^{2}+(k_{l}^{2}+k_{1}^{2})\mathbf{k}_{1}\cdot\mathbf{k}_{l}\right]\rightarrow
0~{},$ $\displaystyle
B_{121}\ni\frac{5}{7}D_{A}(z)\left[2k_{l}^{2}k_{1}^{2}-(k_{l}^{2}+k_{1}^{2})\mathbf{k}_{1}\cdot\mathbf{k}_{l}\right]\rightarrow
0~{}$ (B.6)

but they vanish in the limit $k_{l}\rightarrow 0$.

The Fourier kernel for the second term of the standard perturbation theory is
already symmetric:

$\bullet~{}~{}\frac{1}{7}D_{B}(z)\Delta\left({\cal R}^{\alpha}{\cal
R}_{\alpha}\right):\qquad\qquad{\cal
F}_{s}=F_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=\frac{1}{7}D_{B}(z)\left(\mathbf{k}_{1}+\mathbf{k}_{2}\right)^{2}\mathbf{k}_{1}\cdot\mathbf{k}_{2}~{},$
(B.7)

and the one-point ensemble average and the first connected bispectra vanish
due to symmetry of the wave vectors

$\Xi=0~{},\qquad\qquad B_{112}=0~{}.$ (B.8)

The two connected bispectra are obtained by Eqs. (A.70) and (A.71)

$B_{211}\ni\frac{2}{7}D_{B}(z)(\mathbf{k}_{1}+\mathbf{k}_{l})^{2}\mathbf{k}_{1}\cdot\mathbf{k}_{l}\rightarrow
0~{},\qquad\qquad
B_{121}\ni-\frac{2}{7}D_{B}(z)(\mathbf{k}_{1}-\mathbf{k}_{l})^{2}\mathbf{k}_{1}\cdot\mathbf{k}_{l}\rightarrow
0~{},$ (B.9)

and they also vanish in the limit $k_{l}\rightarrow 0$.

In the standard perturbation theory, two growth functions $D_{A}(t)$ and
$D_{B}(t)$ are assumed to be equal to $D_{1}^{2}(t)$, and the equality is only
valid in the Einstein-de Sitter universe, while the approximation yields only
small errors in the late universe [44]. The first kernel in Eq. (B.3) can then
be re-arranged as

${5\over 7}D_{1}^{2}(t)\left[1+{\mathbf{k}_{1}\cdot\mathbf{k}_{2}\over
2k_{1}k_{2}}\left({k_{1}\over k_{2}}+{k_{2}\over k_{1}}\right)\right]~{},$
(B.10)

and the second kernel in Eq. (B.7) is

$\frac{1}{7}D_{1}^{2}(t)\left[{\mathbf{k}_{1}\cdot\mathbf{k}_{2}\over
k_{1}k_{2}}\left({k_{1}\over k_{2}}+{k_{2}\over
k_{1}}\right)+2\left({\mathbf{k}_{1}\cdot\mathbf{k}_{2}\over
k_{1}k_{2}}\right)^{2}\right]~{},$ (B.11)

where we scaled each Fourier kernel by $k_{1}^{2}k_{2}^{2}$ in the
denominator, as the Fourier kernels in the standard perturbation theory are
expressed in terms of the matter density fluctuation $\delta({\bf
k})=-k^{2}{\cal R}({\bf k})$. They add up to yield the standard kernel

$F_{2}(\mathbf{k}_{1},\mathbf{k}_{2})=D_{1}^{2}(t)\left[\frac{5}{7}+\frac{1}{2}{\mathbf{k}_{1}\cdot\mathbf{k}_{2}\over
k_{1}k_{2}}\left({k_{1}\over k_{2}}+{k_{2}\over
k_{1}}\right)+\frac{2}{7}\left({\mathbf{k}_{1}\cdot\mathbf{k}_{2}\over
k_{1}k_{2}}\right)^{2}\right]~{},$ (B.12)

but note that the shift term in proportion to
$\mathbf{k}_{1}\cdot\mathbf{k}_{2}$ originates from two terms with different
time dependence.

Relativistic contributions in the matter density fluctuation: At the second
order in perturbations, there exist two extra relativistic contributions to
the matter density fluctuation in addition to the standard contributions
discussed in the context of the standard perturbation theory. They originate
from the constraint equation of general relativity. The Fourier kernels for
these two terms in Eq. (2.4) are computed here. The Fourier kernel for the
first relativistic contribution is

$\bullet~{}~{}\frac{3}{2}D_{1}(z){\cal R}^{,\alpha}{\cal
R}_{,\alpha}:\qquad\qquad{\cal
F}_{s}=F_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=-\frac{3}{2}D_{1}(z)\mathbf{k}_{1}\cdot\mathbf{k}_{2}~{},$
(B.13)

and its one-point ensemble average is non-vanishing

$\Xi=\frac{3}{2}D_{1}(z)\sigma_{2}\,,$ (B.14)

where the variance $\sigma_{n}$ is defined in Eq. (2.16). The contributions to
the connected bispectra are rather simple

$B_{112}\ni 3D_{1}k_{1}^{2}~{},\qquad\qquad
B_{211}=-B_{121}\ni-3D_{1}\mathbf{k}_{1}\cdot\mathbf{k}_{l}\rightarrow 0~{}.$
(B.15)

The Fourier kernel for the second relativistic contribution in Eq. (2.4) is

$\bullet~{}~{}4D_{1}(z){\cal R}\Delta{\cal R}:\qquad\qquad
F_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=-4D_{1}(z)k^{2}_{2}~{},$ (B.16)

and its symmetrized Fourier kernel is

${\cal
F}_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=-2D_{1}\left(k_{1}^{2}+k_{2}^{2}\right)~{}.$
(B.17)

For the bispectrum contributions, the kernels $F_{s}$ get symmetrized, so we
can use the symmetric kernel ${\cal F}_{s}$, but for the one-point ensemble
average in Eq. (A.41), the kernel $F_{s}$ should be used. The one-point
ensemble average of the second relativistic contribution is

$\Xi=-4D_{1}(z)\sigma_{2}~{}.$ (B.18)

Finally, the contributions to the connected bispectra are then obtained by
using Eqs. (A.69)$-$(A.71)

$B_{112}\ni-8D_{1}k_{1}^{2}~{},\qquad\qquad
B_{211}=B_{121}\ni-4D_{1}(k_{1}^{2}+k_{l}^{2})\rightarrow-4D_{1}k_{1}^{2}~{},$
(B.19)

and all of them are non-zero in the limit $k_{l}\rightarrow 0$.

Non-Gaussian contributions in the presence of $f_{\rm NL}$: The primordial
non-Gaussianity in the initial condition with non-vanishing $f_{\rm NL}$ in
Eq. (2.7) gives rise two extra terms in Eq. (2.25), and their Fourier kernels
can be readily read-off from Eqs. (B.13) and (B.16) as

$\displaystyle\bullet~{}~{}-\frac{6}{5}D_{1}(z)f_{\rm NL}~{}{\cal
R}^{,\alpha}{\cal R}_{\,a}:\qquad\quad{\cal
F}_{s}=F_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=\frac{6}{5}D_{1}(z)f_{\rm
NL}~{}\mathbf{k}_{1}\cdot\mathbf{k}_{2}~{},$ (B.20)
$\displaystyle\bullet~{}~{}-\frac{6}{5}D_{1}(z)f_{\rm NL}~{}{\cal
R}\Delta{\cal R}:\qquad\qquad{\cal
F}_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=\frac{3}{5}D_{1}(z)f_{\rm
NL}\left(k_{1}^{2}+k_{2}^{2}\right)~{}.$ (B.21)

As opposed to the two relativistic contributions in the matter density
fluctuations, the coefficients for these two contributions arrange in a way
that the sum of these two contributions in proportion to $f_{\rm NL}$ is

${\cal F}_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=\frac{3}{5}D_{1}(z)f_{\rm
NL}|\mathbf{k}_{1}+\mathbf{k}_{2}|^{2}~{},$ (B.22)

and the one-point ensemble average vanishes by cancellation of the two terms.
The contribution to the first bispectrum vanishes, and the other two bispectra

$B_{211}\ni\frac{6}{5}D_{1}f_{\rm nl}|\mathbf{k}_{1}+{\bf
k}_{l}|^{2}\quad\rightarrow~{}~{}\frac{6}{5}D_{1}f_{\rm
nl}k_{1}^{2}~{},\qquad\qquad B_{121}\ni\frac{6}{5}D_{1}f_{\rm
nl}|\mathbf{k}_{1}-{\bf
k}_{l}|^{2}\quad\rightarrow~{}~{}\frac{6}{5}D_{1}f_{\rm nl}k_{1}^{2}~{},$
(B.23)

survive in the squeezed limit.

Coupling contributions at the source position: The linear-order matter density
fluctuation is coupled with the relativistic effects in the light propagation,
which appears as $\Delta x^{\mu}~{}\partial_{\mu}\delta$ in Eq. (4.13). This
coupling results in numerous extra terms detailed in Eq. (A.1), as there are
many contribution terms in $\Delta x^{\mu}$. Among those, there are four
terms, contributing at the source position, and we compute the Fourier kernels
for those four terms.

The first two terms arise from the coupling
$\delta\eta~{}\partial_{\eta}\delta$, where $\delta\eta=\delta z/\mathcal{H}$
and $\delta^{\prime}\propto D_{1}^{\prime}=D_{V}$. The first of the two
contributions is the Sachs-Wolfe contribution at the source, and its Fourier
kernel is

$\bullet~{}~{}\frac{D_{V}}{\mathcal{H}}{\cal R}\Delta{\cal R}:\qquad\qquad
F_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=-{D_{V}(z)\over\mathcal{H}_{z}}k_{2}^{2}~{},\qquad\qquad{\cal
F}_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=-\frac{D_{V}(z)}{2\mathcal{H}_{z}}\left(k_{1}^{2}+k_{2}^{2}\right)~{}.$
(B.24)

The one-point ensemble average of this coupling term is

$\Xi=-\frac{D_{V}(z)}{\mathcal{H}_{z}}\sigma_{2}~{},$ (B.25)

and the connected bispectra are

$B_{112}\ni-{2D_{V}(z)\over\mathcal{H}_{z}}k_{1}^{2}~{},\qquad\qquad
B_{211}=B_{121}\ni-{D_{V}(z)\over\mathcal{H}_{z}}(k_{1}^{2}+k_{l}^{2})\rightarrow{D_{V}(z)\over\mathcal{H}_{z}}k_{1}^{2}~{}.$
(B.26)

None of the one-point ensemble average or the connected bispectra vanish in
the limit $k_{l}\rightarrow 0$.

The second contribution in the coupling $\delta\eta~{}\partial_{\eta}\delta$
arises from the line-of-sight velocity contribution at the source position.
The Fourier kernel for the contribution is

$\bullet~{}~{}-\frac{D_{V}^{2}}{\mathcal{H}_{z}}\partial_{r}{\cal
R}\Delta{\cal R}:\qquad\qquad
F_{s}(\mathbf{k}_{1},\mathbf{k}_{2})={D_{V}^{2}(z)\over\mathcal{H}_{z}}i\mu_{1}k_{1}k_{2}^{2}~{},$
(B.27)

and its symmetrized Fourier kernel is

${\cal
F}_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=\frac{D_{V}^{2}(z)}{2\mathcal{H}_{z}}ik_{1}k_{2}\left(\mu_{1}\,k_{2}+\mu_{2}k_{1}\right)~{}.$
(B.28)

The one-point ensemble average and the first connected bispectrum then
trivially vanish due to the symmetry of the wave vectors

$\Xi=0~{},\qquad\qquad B_{112}=0~{},$ (B.29)

and the remaining two connected bispectra also vanish

$B_{211}\ni{D_{V}^{2}\over\mathcal{H}_{z}}ik_{1}k_{l}\left(\mu_{1}k_{l}+\mu_{l}k_{1}\right)\rightarrow
0~{},\qquad\qquad
B_{121}\ni{D_{V}^{2}\over\mathcal{H}_{z}}ik_{1}k_{l}\left(\mu_{l}k_{1}-\mu_{1}k_{l}\right)\rightarrow
0~{},$ (B.30)

upon angle average in the limit.

The remaining two contributions at the source position come from the coupling
term $\delta r~{}\partial_{r}\delta$ (no contributions at the source position
from the angular distortion). The first term is the coupling of the matter
density fluctuation and the gravitational potential at the source position.
The Fourier kernel for this contribution is

$\bullet~{}~{}\frac{D_{1}D_{\Psi}}{\mathcal{H}}{\cal R}\partial_{r}\Delta{\cal
R}:\qquad\qquad
F_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=-{D_{1}(z)D_{\Psi}(z)\over\mathcal{H}_{z}}i\mu_{2}k_{2}^{3}~{},$
(B.31)

and its symmetrized kernel is

${\cal
F}_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=-\frac{D_{1}(z)D_{\Psi}(z)}{2\mathcal{H}_{z}}i\left(\mu_{1}\,k_{1}^{3}+\mu_{2}\,k_{2}^{3}\right)~{}.$
(B.32)

The one-point ensemble average and the first connected bispectrum vanish

$\Xi=0~{},\qquad\qquad B_{112}=0~{},$ (B.33)

upon angle average. Mind that the one-point ensemble average is computed by
using $F_{s}$ in Eq. (A.41), while the connected bispectra are computed by
using ${\cal F}_{s}$ in Eq. (A.69). The remaining two connected bispectra

$B_{211}\ni-{D_{1}D_{\Psi}\over\mathcal{H}_{z}}i\left(\mu_{1}k_{1}^{3}+\mu_{l}k_{l}^{3}\right)\rightarrow
0~{},\qquad\qquad
B_{121}\ni-{D_{1}D_{\Psi}\over\mathcal{H}_{z}}i\left(-\mu_{1}k_{1}^{3}+\mu_{l}k_{l}^{3}\right)\rightarrow
0~{},$ (B.34)

also vanish upon angle average in the limit $k_{l}\rightarrow 0$.

The last contribution is the coupling with the line-of-sight velocity at the
source position. The Fourier kernel for the fourth term is

$\bullet~{}~{}\frac{D_{1}D_{V}}{\mathcal{H}_{z}}\partial_{r}R\partial_{r}\left(\Delta
R\right):\qquad\qquad
F_{s}(\mathbf{k}_{1},\mathbf{k}_{2})={D_{1}(z)D_{V}(z)\over\mathcal{H}_{z}}\mu_{1}k_{1}\mu_{2}k_{2}^{3}~{},$
(B.35)

and its symmetrized kernel is

${\cal
F}_{s}(\mathbf{k}_{1},\mathbf{k}_{2})=\frac{D_{1}(z)D_{V}(z)}{2\mathcal{H}_{z}}\mu_{1}\mu_{2}\,k_{1}k_{2}\left(k_{1}^{2}+k_{2}^{2}\right)~{}.$
(B.36)

The one-point ensemble average and the first connected bispectrum are non-zero

$\Xi=-\frac{D_{1}(z)D_{V}(z)}{3\mathcal{H}_{z}}\sigma_{4}\,,\qquad\qquad
B_{112}=-{2D_{1}(z)D_{V}(z)\over\mathcal{H}_{z}}\mu_{1}^{2}k_{1}^{4}\rightarrow-{2D_{1}(z)D_{V}(z)\over
3\mathcal{H}_{z}}k_{1}^{4}~{},$ (B.37)

where the angular dependence is removed upon angular integration. The
remaining two connected bispectra, however,

$B_{211}=-B_{121}\ni{D_{1}(z)D_{V}(z)\over\mathcal{H}_{z}}\mu_{1}\mu_{l}k_{1}k_{l}(k_{1}^{2}+k_{l}^{2})\rightarrow
0~{},$ (B.38)

vanish, if we take the limit $k_{l}\rightarrow 0$.

###  B.2 Coupling contributions from the source and the observer positions

Since cosmological observables are measured by the observer, these observables
depend not only on the physical properties of the source, but also on the
state of the observer. As the observers perform measurements in the rest frame
different from the FRW frame, there exist various contributions in the
cosmological observables from the gravitational potential and the line-of-
sight velocity at the observer position, as shown in Eqs. (4.14)$-$(4.15).
These contributions contribute only through the coupling terms in $\Delta
x^{\mu}~{}\partial_{\mu}\delta$ in Eq. (4.13). As discussed in Appendix A.2,
the Fourier kernels ${\cal F}_{o}$ for such coupling contributions involve the
exponential factor in Eq. (A.45), in addition to the ordinary Fourier kernels
$F_{o}$ from the coupling terms.

There exist six different contributions in Eq. (A.1), involving those at the
observer position. The first two terms come from the coupling
$\delta\eta~{}\partial_{\eta}\delta$. The Sachs-Wolfe contribution at the
observer position makes the first of such terms, and its Fourier kernel is

$\bullet~{}~{}-\frac{D_{V}}{\mathcal{H}}{\cal R}_{o}\Delta{\cal
R}:\qquad\qquad
F_{o}(\mathbf{k}_{1},\mathbf{k}_{2})=\frac{D_{V}(z)}{\mathcal{H}_{z}}k_{2}^{2}~{},$
(B.39)

and its symmetrized kernel with the exponential factor is

${\cal F}_{o}={D_{V}(z)\over
2\mathcal{H}_{z}}\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{2}^{2}+e^{-i\mathbf{k}_{2}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{1}^{2}\right)~{},$
(B.40)

where the subscript $o$ indicates that ${\cal R}$ is evaluated at the observer
position and there is no time dependence of ${\cal R}_{o}$ on
$\bar{\eta}_{o}$, as the initial condition is constant in time. Furthermore,
we keep the convention that the first wave vector $\mathbf{k}_{1}$ describes
the quantity at the observer position, while the second wave vector
$\mathbf{k}_{2}$ describes one at the source position. The one-point ensemble
average and the first connected bispectrum are

$\Xi(z)={D_{V}(z)\over{\cal H}_{z}}\sigma_{2,0}(z)~{},\qquad\qquad
B_{112}\ni{2D_{V}(z)\over\mathcal{H}_{z}}k_{1}^{2}j_{0}(k_{1}\bar{r}_{z})~{},$
(B.41)

where the time-dependent variance $\sigma_{n,m}(z)$ is defined in Eq. (4.20).
Furthermore, in the derivation we used

$\int_{-1}^{1}\mathrm{d\mu\over 2}~{}e^{\pm ix\mu}=j_{0}(x)~{},$ (B.42)

and also for the following derivations we will use the useful relation
involving the integration of the exponential factor:

$\int_{-1}^{1}\mathrm{d\mu\over 2}~{}\mu e^{\pm ix\mu}=\pm
ij_{1}(x)~{},\qquad\qquad\int_{-1}^{1}\mathrm{d\mu\over 2}~{}\mu^{2}e^{\pm
ix\mu}={j_{0}(x)-2j_{2}(x)\over 3}~{}.$

Two connected bispectra are non-vanishing

$\displaystyle B_{211}$ $\displaystyle\ni$
$\displaystyle{D_{V}(z)\over\mathcal{H}_{z}}\left(e^{-i\mathbf{k}_{1}\cdot{\bar{r}_{z}}{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot{\bar{r}_{z}}{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow{D_{V}(z)\over\mathcal{H}_{z}}k_{1}^{2}~{},$
(B.43) $\displaystyle B_{121}$ $\displaystyle\ni$
$\displaystyle{D_{V}(z)\over\mathcal{H}_{z}}\left(e^{i\mathbf{k}_{1}\cdot{\bar{r}_{z}}{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot{\bar{r}_{z}}{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow{D_{V}(z)\over\mathcal{H}_{z}}k_{1}^{2}~{},$
(B.44)

even in the limit $k_{l}\rightarrow 0$.

The second term arises from the line-of-sight velocity at the observer
position, contributing to $\delta\eta$. The Fourier kernel for this term is

$\bullet~{}~{}\frac{D_{V}D_{V}^{o}}{\mathcal{H}}(\partial_{r}{\cal
R})_{o}\Delta{\cal R}:\qquad\qquad
F_{o}(\mathbf{k}_{1},\mathbf{k}_{2})=-\frac{D_{V}(z)D_{V}(\bar{\eta}_{o})}{\mathcal{H}_{z}}i\mu_{1}k_{1}k_{2}^{2}\,,$
(B.45)

and its symmetrized kernel is

${\cal
F}_{o}=-\frac{D_{V}(z)D_{V}(\bar{\eta}_{o})}{2\mathcal{H}_{z}}ik_{1}k_{2}\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{1}k_{2}+e^{-i\mathbf{k}_{2}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{2}k_{1}\right)~{}.$
(B.46)

The one-point ensemble average and the first connected bispectrum are

$\Xi=-{D_{V}(z)D_{V}(\bar{\eta}_{o})\over{\cal
H}_{z}}\sigma_{3,1}(z)~{},\qquad\qquad
B_{112}\ni-\frac{2D_{V}(z)D_{V}(\bar{\eta}_{o})}{\mathcal{H}_{z}}k_{1}^{3}j_{1}(k_{1}\bar{r}_{z})~{},$
(B.47)

non-vanishing, while two remaining connected bispectra

$\displaystyle B_{211}$ $\displaystyle\ni$
$\displaystyle-{D_{V}(z)D_{V}(\bar{\eta}_{o})\over\mathcal{H}_{z}}ik_{1}k_{l}\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{1}k_{l}+e^{-i\mathbf{k}_{l}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{l}k_{1}\right)\rightarrow
0~{},$ (B.48) $\displaystyle B_{121}$ $\displaystyle\ni$
$\displaystyle-{D_{V}(z)D_{V}(\bar{\eta}_{o})\over\mathcal{H}_{z}}ik_{1}k_{l}\left(-e^{i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{1}k_{l}+e^{-i\mathbf{k}_{l}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{l}k_{1}\right)\rightarrow
0~{},$ (B.49)

vanish in the limit $k_{l}\rightarrow 0$, upon angular integration.

The next two terms come from the coupling term $\delta
r~{}\partial_{r}\delta$, in which the gravitational potential and the line-of-
sight velocity terms couple to the derivative of the matter density
fluctuation. The Fourier kernel for the gravitational potential contribution
at the observer position is

$\bullet~{}~{}-D_{1}\left(D_{V}^{o}-\frac{1}{\mathcal{H}_{z}}\right){\cal
R}_{o}\partial_{r}\left(\Delta{\cal R}\right):\qquad\qquad
F_{o}(\mathbf{k}_{1},\mathbf{k}_{2})=D_{1}(z)\left[D_{V}(\bar{\eta}_{o})-\frac{1}{\mathcal{H}_{z}}\right]i\mu_{2}k_{2}^{3}\,,$
(B.50)

and its symmetrized kernel is

${\cal F}_{o}={D_{1}(z)\over
2}\left[D_{V}(\bar{\eta}_{o})-{1\over\mathcal{H}_{z}}\right]i\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{2}k_{2}^{3}+e^{-i\mathbf{k}_{2}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{1}k_{1}^{3}\right)~{}.$
(B.51)

Evident in the time-dependent pre-factors, there are multiple contributions
from the gravitational potential at the observer position with different pre-
factors. The one-point ensemble average and the first connected bispectrum are

$\Xi(z)=-D_{1}(z)\left[D_{V}(\bar{\eta}_{o})-{1\over\mathcal{H}_{z}}\right]\sigma_{3,1}(z)~{},\qquad\qquad
B_{112}\ni-2D_{1}(z)\left[D_{V}(\bar{\eta}_{o})-{1\over\mathcal{H}_{z}}\right]k_{1}^{3}j_{1}(k_{1}\bar{r}_{z})~{},$
(B.52)

and the remaining connected bispectra are

$\displaystyle B_{211}$ $\displaystyle\ni$ $\displaystyle
D_{1}(z)\left[D_{V}(\bar{\eta}_{o})-{1\over\mathcal{H}_{z}}\right]i\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{l}k_{l}^{3}+e^{-i\mathbf{k}_{l}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{1}k_{1}^{3}\right)\rightarrow
0~{},$ (B.53) $\displaystyle B_{121}$ $\displaystyle\ni$ $\displaystyle
D_{1}(z)\left[D_{V}(\bar{\eta}_{o})-{1\over\mathcal{H}_{z}}\right]i\left(e^{i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{l}k_{l}^{3}-e^{-i\mathbf{k}_{l}\cdot\bar{r}_{z}{\bf\hat{n}}}\mu_{1}k_{1}^{3}\right)\rightarrow
0~{},$ (B.54)

vanishing in the limit $k_{l}\rightarrow 0$ upon angular integration.

The next term is the line-of-sight velocity contribution to the coupling term
at the observer position, and its Fourier kernel is

$\bullet~{}~{}-\frac{D_{1}D_{V}^{o}}{\mathcal{H}_{z}}(\partial_{r}{\cal
R}_{o})\partial_{r}\left(\Delta{\cal R}\right):\qquad\qquad
F_{o}(\mathbf{k}_{1},\mathbf{k}_{2})=-\frac{D_{1}(z)D_{V}(\bar{\eta}_{o})}{\mathcal{H}_{z}}\mu_{1}\mu_{2}k_{1}k_{2}^{3}\,,$
(B.55)

and its symmetrized kernel is

${\cal F}_{o}=-{D_{1}(z)D_{V}(\bar{\eta}_{o})\over
2\mathcal{H}_{z}}\mu_{1}\mu_{2}k_{1}k_{2}\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{2}^{2}+e^{-i\mathbf{k}_{2}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{1}^{2}\right)~{},$
(B.56)

The one-point ensemble average and the first connected bispectrum are

$\Xi(z)={D_{1}(z)D_{V}(\bar{\eta}_{o})\over
3\mathcal{H}_{z}}\left[\sigma_{4,0}(z)-2\sigma_{4,2}(z)\right]~{},\qquad\qquad
B_{112}\ni{2D_{1}(z)D_{V}(\bar{\eta}_{o})\over
3\mathcal{H}_{z}}k_{1}^{4}\left[j_{0}(k_{1}\bar{r}_{z})-2j_{2}(k_{1}\bar{r}_{z})\right]~{},$
(B.57)

non-vanishing, while the remaining connected bispectra vanish

$\displaystyle B_{211}$ $\displaystyle\ni$
$\displaystyle-{D_{1}(z)D_{V}(\bar{\eta}_{o})\over\mathcal{H}_{z}}\mu_{1}\mu_{l}k_{1}k_{l}\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow
0~{},$ (B.58) $\displaystyle B_{121}$ $\displaystyle\ni$
$\displaystyle{D_{1}(z)D_{V}(\bar{\eta}_{o})\over\mathcal{H}_{z}}\mu_{1}\mu_{l}k_{1}k_{l}\left(e^{i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow
0~{}.$ (B.59)

The other coupling contribution
$(\delta\theta~{}\partial_{\theta}+\delta\phi~{}\partial_{\phi})\delta_{m}$
comes from the angular distortion, which includes only the line-of-sight
velocity term at the observer position. The Fourier kernel for this term is

$\bullet~{}~{}D_{1}D_{V}^{o}(\partial^{\alpha}{\cal
R})_{o}\hat{\nabla}_{\alpha}(\Delta{\cal R}):\qquad\qquad
F_{o}(\mathbf{k}_{1},\mathbf{k}_{2})=D_{1}(z)D_{V}(\bar{\eta}_{o}){\bar{r}_{z}}\left(-\mu_{1}\mu_{2}k_{1}k_{2}+\mathbf{k}_{1}\cdot\mathbf{k}_{2}\right)k_{2}^{2}\,,$
(B.60)

and its symmetrized kernel is

${\cal
F}_{o}=\frac{1}{2}D_{1}(z)D_{V}(\bar{\eta}_{o})\bar{r}_{z}\left(-\mu_{1}\mu_{2}k_{1}k_{2}+\mathbf{k}_{1}\cdot\mathbf{k}_{2}\right)\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{2}^{2}+e^{-i\mathbf{k}_{2}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{1}^{2}\right)~{}.$
(B.61)

Note that there exist two contributions in total, one from $\delta\theta$ and
one from $\delta\phi$, which are combined as the angular gradient
$\hat{\nabla}$. The one-point ensemble average and the first connected
bispectrum are

$\Xi=-2D_{1}(z)D_{V}(\bar{\eta}_{o})\sigma_{3,1}(z)~{},\qquad\qquad
B_{112}\ni-4D_{1}(z)D_{V}(\bar{\eta}_{o})k_{1}^{3}j_{1}(k_{1}\bar{r}_{z})~{},$
(B.62)

non-zero, while the remaining connected bispectra are

$\displaystyle B_{211}$ $\displaystyle\ni$ $\displaystyle
D_{1}D_{V}^{o}\bar{r}_{z}\left(-\mu_{1}\mu_{l}k_{1}k_{l}+\mathbf{k}_{1}\cdot\mathbf{k}_{l}\right)\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow
0~{},$ (B.63) $\displaystyle B_{121}$ $\displaystyle\ni$ $\displaystyle
D_{1}D_{V}^{o}\bar{r}_{z}\left(\mu_{1}\mu_{l}k_{1}k_{l}-\mathbf{k}_{1}\cdot\mathbf{k}_{l}\right)\left(e^{i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow
0~{},$ (B.64)

vanishing in the limit $k_{l}\rightarrow 0$, where we used
$j_{1}(x)=x(j_{0}+j_{2})/3$.

One last contribution is from the spatial shift $\delta x^{\alpha}_{o}$ of the
observer position, contained in spatial distortions $\delta r$,
$\delta\theta$, and $\delta\phi$. While it contains an integration over time,
its spatial position is fixed at the observer position at the perturbation
order of our interest, rendering it essentially the same as any other
contributions at the observer position in this section. The Fourier kernel for
the spatial shift contribution is

$\bullet~{}~{}-D_{1}\left(\int_{0}^{\bar{\eta}_{o}}\mathrm{d}\bar{\eta}\,D_{V}(\bar{\eta})\nabla_{\alpha}{\cal
R}\right)\nabla^{\alpha}(\Delta{\cal R}):\qquad\qquad
F_{o}(\mathbf{k}_{1},\mathbf{k}_{2})=-D_{1}(z)\left(\int_{0}^{\bar{\eta}_{o}}\mathrm{d}\bar{\eta}\,D_{V}(\bar{\eta})\right)\left(\mathbf{k}_{1}\cdot\mathbf{k}_{2}\right)k_{2}^{2}~{},$
(B.65)

and its symmetrized kernel is

${\cal F}_{o}=-{D_{1}(z)\over
2}\left(\int_{0}^{\bar{\eta}_{o}}\mathrm{d}\bar{\eta}\,D_{V}(\bar{\eta})\right)\left(\mathbf{k}_{1}\cdot\mathbf{k}_{2}\right)\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{2}^{2}+e^{-i\mathbf{k}_{2}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{1}^{2}\right)~{},$
(B.66)

where the integration over $\bar{\eta}$ is along the time coordinate with
spatial position fixed (not the line-of-sight integration). The one-point
ensemble average and the first connected bispectrum are

$\Xi(z)=D_{1}(z)\sigma_{4,0}(z)\int_{0}^{\bar{\eta}_{o}}\mathrm{d}\bar{\eta}~{}D_{V}(\bar{\eta})~{},\qquad\qquad
B_{112}\ni
2D_{1}(z)\left(\int_{0}^{\bar{\eta}_{o}}\mathrm{d}\bar{\eta}\,D_{V}(\bar{\eta})\right)k_{1}^{4}j_{0}(k_{1}\bar{r}_{z})~{},$
(B.67)

and the remaining connected bispectra are

$\displaystyle B_{211}$ $\displaystyle\ni$
$\displaystyle-{D_{1}(z)}\left(\int_{0}^{\bar{\eta}_{o}}\mathrm{d}\bar{\eta}\,D_{V}(\bar{\eta})\right)\mathbf{k}_{1}\cdot\mathbf{k}_{l}\left(e^{-i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow
0~{},$ (B.68) $\displaystyle B_{121}$ $\displaystyle\ni$
$\displaystyle{D_{1}(z)}\left(\int_{0}^{\bar{\eta}_{o}}\mathrm{d}\bar{\eta}\,D_{V}(\bar{\eta})\right)\mathbf{k}_{1}\cdot\mathbf{k}_{l}\left(e^{i\mathbf{k}_{1}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot\bar{r}_{z}{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow
0~{}.$ (B.69)

In the Einstein-de Sitter universe, the pre-factors are further simplified as

$-D_{1}(z)\left[D_{V}(\bar{\eta}_{o})-\frac{1}{\mathcal{H}_{z}}\right]={\bar{\eta}_{z}^{2}\over
10}\left({\bar{\eta}_{o}\over 5}-{\bar{\eta}_{z}\over
2}\right)~{},\qquad\qquad-
D_{1}\left(\int_{0}^{\bar{\eta}_{o}}\mathrm{d}\bar{\eta}\,D_{V}(\bar{\eta})\right)=-{\bar{\eta}_{z}^{2}\bar{\eta}_{o}^{2}\over
100}~{}.$ (B.70)

###  B.3 Coupling contributions involving the line-of-sight integral

The last type of contributions to the observed matter density fluctuation
involves the coupling with contributions from the line-of-sight direction. The
well-known components in this category are the contributions of the
gravitational lensing and the integrated Sachs-Wolfe effects. The Fourier
kernels for such contributions in Eq. (A.49) include not only the line-of-
sight integration, but also the exponential factor. There exist four coupling
terms arising from the line-of-sight integration.

The first term arises from the coupling term in
$\delta\eta~{}\partial_{\eta}\delta$, and the integrated Sachs-Wolfe
contribution in $\delta\eta=\delta z/\mathcal{H}$ results in the coupling
contribution involving the line-of-sight integral. Its Fourier kernel is

$\bullet~{}~{}\frac{2D_{V}}{\mathcal{H}}\left(\int_{0}^{\bar{r}_{z}}d\bar{r}~{}D_{\Psi}(\bar{r})\partial_{r}{\cal
R}\right)\Delta{\cal R}:\qquad\qquad
F_{nl}(\bar{r},\mathbf{k}_{1},\mathbf{k}_{2})=-\frac{2D_{V}(z)}{\mathcal{H}_{z}}D_{\Psi}(\bar{r})i\mu_{1}k_{1}k_{2}^{2}\,,$
(B.71)

where the first wave vector $\mathbf{k}_{1}$ belongs to the contribution along
the line-of-sight direction and the second wave vector $\mathbf{k}_{2}$
describes the matter fluctuation at the source position. The full and
symmetrized Fourier kernel is

${\cal
F}_{nl}=-\frac{D_{V}(z)}{\mathcal{H}_{z}}\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})ik_{1}k_{2}\left(e^{-i\mathbf{k}_{1}\Delta
r{\bf\hat{n}}}\mu_{1}k_{2}+e^{-i\mathbf{k}_{2}\Delta
r{\bf\hat{n}}}\mu_{2}k_{1}\right)\,,$ (B.72)

and its one-point ensemble average is

$\Xi(z)=-{2D_{V}(z)\over{\cal
H}_{z}}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}~{}k^{3}P_{\cal
R}(k)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})j_{1}(k\bar{r})~{},$
(B.73)

where $\Delta r:={\bar{r}_{z}}-\bar{r}$. The connected bispectra are then
obtained as

$\displaystyle B_{112}$ $\displaystyle\ni$
$\displaystyle-\frac{4D_{V}(z)}{\mathcal{H}_{z}}\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})k_{1}^{3}j_{1}(k_{1}\Delta
r)~{},$ (B.74) $\displaystyle B_{211}$ $\displaystyle\ni$
$\displaystyle-{2D_{V}(z)\over\mathcal{H}_{z}}\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\psi}(\bar{r})ik_{1}k_{l}\left(e^{-i\mathbf{k}_{1}\cdot\Delta
r{\bf\hat{n}}}\mu_{1}k_{l}+e^{-i\mathbf{k}_{l}\cdot\Delta
r{\bf\hat{n}}}\mu_{l}k_{1}\right)\rightarrow 0~{},$ (B.75) $\displaystyle
B_{121}$ $\displaystyle\ni$
$\displaystyle-{2D_{V}(z)\over\mathcal{H}_{z}}\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\psi}(\bar{r})ik_{1}k_{l}\left(-e^{i\mathbf{k}_{1}\cdot\Delta
r{\bf\hat{n}}}\mu_{1}k_{l}+e^{-i\mathbf{k}_{l}\cdot\Delta
r{\bf\hat{n}}}\mu_{l}k_{1}\right)\rightarrow 0~{}.$ (B.76)

In the Einstein-de Sitter universe, the gravitational potential is constant in
time:

$D_{\Psi}=-\frac{3}{5}~{},$ (B.77)

allowing us to analytically perform the line-of-sight integration by using

$k\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,j_{1}(k\Delta
r)=1-j_{0}(k{\bar{r}_{z}})\,,\qquad\qquad\frac{k}{3}\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,\Big{[}j_{0}(k\Delta
r)-2j_{2}(k\Delta r)\Big{]}=j_{1}(k{\bar{r}_{z}})\,.$ (B.78)

The one-point ensemble average and the first connected bispectrum becomes

$\Xi(z)=\frac{3\bar{\eta}_{z}^{2}}{25}\left[\sigma_{2}-\sigma_{2,0}(z)\right]\,,\qquad\qquad
B_{112}\ni{24\over 25}{k_{1}^{2}\left[1-j_{0}(k_{1}{\bar{r}_{z}})\right]\over
H_{0}^{2}(1+z)}~{}.$ (B.79)

The second term comes again from the integrated Sachs-Wolfe contribution in
$\delta z$, but through the coupling term $\delta r~{}\partial_{r}\delta$.
Hence the line-of-sight integration is identical to the first term, but with
different derivative coupling to the matter density fluctuation. The Fourier
kernel for this contribution is

$\bullet~{}~{}-\frac{2D_{1}}{\mathcal{H}}\left(\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,D_{\Psi}(\bar{r})\partial_{r}{\cal
R}\right)\partial_{r}(\Delta{\cal R}):\qquad\qquad
F_{nl}(\bar{r},\mathbf{k}_{1},\mathbf{k}_{2})=-\frac{2D_{1}(z)}{\mathcal{H}_{z}}D_{\Psi}(\bar{r})\mu_{1}\mu_{2}k_{1}k_{2}^{3}\,,$
(B.80)

and its symmetrized kernel is

${\cal
F}_{nl}=-\frac{D_{1}(z)}{\mathcal{H}_{z}}\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})\mu_{1}\mu_{2}k_{1}k_{2}\left(e^{-i\mathbf{k}_{1}\Delta
r{\bf\hat{n}}}k_{2}^{2}+e^{-i\mathbf{k}_{2}\Delta
r{\bf\hat{n}}}k_{1}^{2}\right)\,.$ (B.81)

The one-point ensemble average is

$\Xi(z)={2D_{1}(z)\over 3{\cal
H}_{z}}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}~{}k^{4}P_{\cal
R}(k)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})\left[j_{0}(k\bar{r})-2j_{2}(k\bar{r})\right]~{},$
(B.82)

and the first connected bispectrum is

$B_{112}\ni{4D_{1}(z)\over
3\mathcal{H}_{z}}\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})k_{1}^{4}\left[j_{0}(k_{1}\Delta
r)-2j_{2}(k_{1}\Delta r)\right]~{}.$ (B.83)

Two remaining connected bispectra

$\displaystyle B_{211}$ $\displaystyle\ni$
$\displaystyle-{2D_{1}(z)\over\mathcal{H}_{z}}\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})\mu_{1}\mu_{l}k_{1}k_{l}\left(e^{-i\mathbf{k}_{1}\cdot\Delta
r{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot\Delta
r{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow 0~{},$ (B.84) $\displaystyle
B_{121}$ $\displaystyle\ni$
$\displaystyle{2D_{1}(z)\over\mathcal{H}_{z}}\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})\mu_{1}\mu_{l}k_{1}k_{l}\left(e^{i\mathbf{k}_{1}\cdot\Delta
r{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot\Delta
r{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow 0~{},$ (B.85)

vanish in the limit $k_{l}\rightarrow 0$ upon angular integration. In the
Einstein-de Sitter universe, the one-point ensemble average and the first
connected bispectrum are further simplified as

$\Xi=-\frac{3\bar{\eta}_{z}^{3}}{50}\sigma_{3,1}(z)~{},\qquad\qquad
B_{112}\ni-{24\over 25}{k_{1}^{3}j_{1}(k_{1}{\bar{r}_{z}})\over
H_{0}^{3}(1+z)^{3/2}}~{}.$ (B.86)

The third term arises from the line-of-sight integration of the gravitational
potential contribution in the radial distortion $\delta r$. Its Fourier kernel
is

$\bullet~{}~{}-2D_{1}\left(\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,D_{\Psi}(\bar{r}){\cal
R}\right)\partial_{r}(\Delta{\cal R}):\qquad\qquad
F_{nl}(\bar{r},\mathbf{k}_{1},\mathbf{k}_{2})=2D_{1}(z)D_{\Psi}(\bar{r})i\mu_{2}k_{2}^{3}\,,$
(B.87)

and its symmetrized kernel is

${\cal
F}_{nl}=D_{1}(z)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})i\left(e^{-i\mathbf{k}_{1}\Delta
r{\bf\hat{n}}}\mu_{2}k_{2}^{3}+e^{-i\mathbf{k}_{2}\Delta
r{\bf\hat{n}}}\mu_{1}k_{1}^{3}\right)\,.$ (B.88)

The one-point ensemble average is

$\Xi(z)=-2D_{1}(z)\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}~{}k^{3}P_{\cal
R}(k)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})j_{1}(k\bar{r})~{},$
(B.89)

and the first connected bispectrum is

$B_{112}\ni-4D_{1}(z)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})k_{1}^{3}j_{1}(k_{1}\Delta
r)~{}.$ (B.90)

The remaining connected bispectra are

$\displaystyle B_{211}$ $\displaystyle\ni$ $\displaystyle
2D_{1}(z)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})i\left(e^{-i\mathbf{k}_{1}\cdot\Delta
r{\bf\hat{n}}}\mu_{l}k_{l}^{3}+e^{-i\mathbf{k}_{l}\cdot\Delta
r{\bf\hat{n}}}\mu_{1}k_{1}^{3}\right)\rightarrow 0~{},$ (B.91) $\displaystyle
B_{121}$ $\displaystyle\ni$ $\displaystyle
2D_{1}(z)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})i\left(e^{i\mathbf{k}_{1}\cdot\Delta
r{\bf\hat{n}}}\mu_{l}k_{l}^{3}-e^{-i\mathbf{k}_{l}\cdot\Delta
r{\bf\hat{n}}}\mu_{1}k_{1}^{3}\right)\rightarrow 0~{}.$ (B.92)

In the Einstein-de Sitter universe, we have

$\displaystyle\Xi(z)=\frac{3\bar{\eta}_{z}^{2}}{25}\left[\sigma_{2}-\sigma_{2,0}(z)\right]\,,\qquad\qquad
B_{112}\ni{24\over 25}{k_{1}^{2}\left[1-j_{0}(k_{1}{\bar{r}_{z}})\right]\over
H_{0}^{2}(1+z)}~{}.$ (B.93)

The last term comes from the gravitational lensing contribution, which couples
to the matter density fluctuation through the angular distortion
$(\delta\theta~{}\partial_{\theta}+\delta\phi~{}\partial_{\phi})\delta_{m}$.
The Fourier kernel for this contribution takes the form

$\displaystyle\bullet~{}~{}2D_{1}\left(\int_{0}^{\bar{r}_{z}}\mathrm{d}\bar{r}\,\left(\frac{{\bar{r}_{z}}-\bar{r}}{{\bar{r}_{z}}\bar{r}}\right)D_{\Psi}(\bar{r})\hat{\nabla}^{\alpha}{\cal
R}\right)\hat{\nabla}_{\alpha}\left(\Delta{\cal R}\right):$ (B.94)
$\displaystyle\qquad\qquad
F_{nl}(\bar{r},\mathbf{k}_{1},\mathbf{k}_{2})=2D_{1}(z)\Delta
rD_{\Psi}(\bar{r})k_{2}^{2}\left(-\mu_{1}\mu_{2}k_{1}k_{2}+\mathbf{k}_{1}\cdot\mathbf{k}_{2}\right)\,,$

and its symmetrized kernel is

${\cal
F}_{nl}=D_{1}(z)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})\Delta
r\left(e^{-i\mathbf{k}_{1}\Delta
r{\bf\hat{n}}}k_{2}^{2}+e^{-i\mathbf{k}_{2}\Delta
r{\bf\hat{n}}}k_{1}^{2}\right)\left(-\mu_{1}\mu_{2}k_{1}k_{2}+\mathbf{k}_{1}\cdot\mathbf{k}_{2}\right)~{}.$
(B.95)

The one-point ensemble average is

$\Xi(z)=-4D_{1}(z)\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}~{}k^{3}P_{\cal
R}(k)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})j_{1}(k\bar{r})~{},$
(B.96)

and the first connected bispectrum is

$B_{112}\ni-8D_{1}(z)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}D_{\Psi}(\bar{r})k_{1}^{3}j_{1}(k_{1}\Delta
r)~{},$ (B.97)

and the remaining connected bispectra are

$\displaystyle B_{211}$ $\displaystyle\ni$ $\displaystyle
2D_{1}(z)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}\Delta
rD_{\Psi}(\bar{r})\left(-\mu_{1}\mu_{l}k_{1}k_{l}+\mathbf{k}_{1}\cdot\mathbf{k}_{l}\right)\left(e^{-i\mathbf{k}_{1}\cdot\Delta
r{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot\Delta
r{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow 0~{},$ (B.98) $\displaystyle
B_{121}$ $\displaystyle\ni$ $\displaystyle
2D_{1}(z)\int_{0}^{{\bar{r}_{z}}}d\bar{r}~{}\Delta
rD_{\Psi}(\bar{r})\left(\mu_{1}\mu_{l}k_{1}k_{l}-\mathbf{k}_{1}\cdot\mathbf{k}_{l}\right)\left(e^{i\mathbf{k}_{1}\cdot\Delta
r{\bf\hat{n}}}k_{l}^{2}+e^{-i\mathbf{k}_{l}\cdot\Delta
r{\bf\hat{n}}}k_{1}^{2}\right)\rightarrow 0~{}.$ (B.99)

In the Einstein-de Sitter universe, the non-vanishing components are further
simplified as

$\Xi={6\bar{\eta}_{z}^{2}\over
25}\left[\sigma_{2}-\sigma_{2,0}(z)\right]~{},\qquad\qquad B_{112}\ni{48\over
25}{k_{1}^{2}\left[1-j_{0}(k_{1}{\bar{r}_{z}})\right]\over
H_{0}^{2}(1+z)}~{}.$ (B.100)

## References

  * [1] M. Tegmark et al., “Cosmological constraints from the SDSS luminous red galaxies,” Phys. Rev. D 74 no. 12, (Dec., 2006) 123507, arXiv:astro-ph/0608632.
  * [2] P. A. R. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont, C. Baccigalupi, A. J. Banday, and et al., “Planck 2013 results. XXIV. Constraints on primordial non-Gaussianity,” Astron. Astrophys. 571 (Nov., 2014) A24, arXiv:1303.5084.
  * [3] Planck Collaboration et al., “Planck 2018 results. VI. Cosmological parameters,” arXiv e-prints (July, 2018) , arXiv:1807.06209.
  * [4] D. H. D. H. Lyth and A. A. Riotto, “Particle physics models of inflation and the cosmological density perturbation,” Phys. Rep. 314 no. 1-2, (Jun, 1999) 1–146, arXiv:hep-ph/9807278.
  * [5] A. Linde, “Inflationary cosmology,” Phys. Rep. 333 (Aug, 2000) 575–591.
  * [6] J. Maldacena, “Non-gaussian features of primordial fluctuations in single field inflationary models,” J. High Energy Phys. 5 (May, 2003) 13, astro-ph/0210603.
  * [7] N. Bartolo, E. Komatsu, S. Matarrese, and A. Riotto, “Non-Gaussianity from inflation: theory and observations,” Phys. Rep. 402 (Nov., 2004) 103–266, arXiv:0406398.
  * [8] N. Dalal, O. Doré, D. Huterer, and A. Shirokov, “Imprints of primordial non-Gaussianities on large-scale structure: Scale-dependent bias and abundance of virialized objects,” Phys. Rev. D 77 no. 12, (June, 2008) 123514, arXiv:0710.4560.
  * [9] L. Verde and S. Matarrese, “Detectability of the Effect of Inflationary Non-Gaussianity on Halo Bias,” Astrophys. J. Lett. 706 (Nov., 2009) L91–L95, arXiv:0909.3224.
  * [10] S. Camera, C. Carbone, C. Fedeli, and L. Moscardini, “Neglecting primordial non-Gaussianity threatens future cosmological experiment accuracy,” Phys. Rev. D 91 no. 4, (Feb., 2015) 043533, arXiv:1412.5172 [astro-ph.CO].
  * [11] S. Camera, R. Maartens, and M. G. Santos, “Einstein’s legacy in galaxy surveys.,” Mon. Not. R. Astron. Soc. 451 (July, 2015) L80–L84, arXiv:1412.4781 [astro-ph.CO].
  * [12] A. Moradinezhad Dizgah and C. Dvorkin, “Scale-dependent galaxy bias from massive particles with spin during inflation,” J. Cosmol. Astropart. Phys. 2018 no. 1, (Jan., 2018) 010, arXiv:1708.06473 [astro-ph.CO].
  * [13] M. Levi, C. Bebek, T. Beers, R. Blum, R. Cahn, D. Eisenstein, B. Flaugher, K. Honscheid, R. Kron, O. Lahav, P. McDonald, N. Roe, et al., “The DESI Experiment, a whitepaper for Snowmass 2013,” ArXiv e-prints (Aug., 2013) , arXiv:1308.0847.
  * [14] C. W. Stubbs, D. Sweeney, J. A. Tyson, and LSST Collaboration, “An Overview of the Large Synoptic Survey Telescope (LSST) System,” in American Astronomical Society Meeting Abstracts, vol. 36 of Bulletin of the American Astronomical Society, p. 108.02.  Dec., 2004.
  * [15] R. Laureijs, J. Amiaux, S. Arduini, J. . Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, et al., “Euclid Definition Study Report,” ArXiv e-prints (Oct., 2011) , arXiv:1110.3193.
  * [16] J. Green, P. Schechter, C. Baltay, R. Bean, D. Bennett, R. Brown, C. Conselice, M. Donahue, et al., “Wide-Field InfraRed Survey Telescope (WFIRST) Final Report,” arXiv:1208.4012.
  * [17] M. Bruni, J. C. Hidalgo, N. Meures, and D. Wands, “Non-Gaussian Initial Conditions in $\Lambda$CDM: Newtonian, Relativistic, and Primordial Contributions,” Astrophys. J. 785 (Apr., 2014) 2, arXiv:1307.1478.
  * [18] M. Bruni, J. C. Hidalgo, and D. Wands, “Einstein’s Signature in Cosmological Large-scale Structure,” Astrophys. J. Lett. 794 (Oct., 2014) L11, arXiv:1405.7006.
  * [19] E. Villa, L. Verde, and S. Matarrese, “General relativistic corrections and non-Gaussianity in large-scale structure,” Classical and Quantum Gravity 31 no. 23, (Dec., 2014) 234005, arXiv:1409.4738 [astro-ph.CO].
  * [20] J. Yoo and J.-O. Gong, “Exact analytic solution for non-linear density fluctuation in a $\Lambda$CDM universe,” J. Cosmol. Astropart. Phys. 7 (July, 2016) 017, arXiv:1602.06300.
  * [21] S. Matarrese, L. Pilo, and R. Rollo, “Resilience of long modes in cosmological observables,” J. Cosmol. Astropart. Phys. 2021 no. 1, (Jan., 2021) 062, arXiv:2007.08877 [astro-ph.CO].
  * [22] E. Pajer, F. Schmidt, and M. Zaldarriaga, “The Observed squeezed limit of cosmological three-point functions,” Phys. Rev. D 88 no. 8, (Oct., 2013) 083502, arXiv:1305.0824.
  * [23] R. de Putter, O. Doré, and D. Green, “Is there scale-dependent bias in single-field inflation?,” J. Cosmol. Astropart. Phys. 2015 no. 10, (Oct., 2015) 024, arXiv:1504.05935 [astro-ph.CO].
  * [24] R. de Putter, O. Doré, D. Green, and J. Meyers, “Single-field inflation and the local ansatz: Distinguishability and consistency,” Phys. Rev. D 95 no. 6, (Mar., 2017) 063501, arXiv:1610.00785 [hep-th].
  * [25] J. Yoo, A. L. Fitzpatrick, and M. Zaldarriaga, “New perspective on galaxy clustering as a cosmological probe: General relativistic effects,” Phys. Rev. D 80 no. 8, (Oct., 2009) 083514, arXiv:0907.0707.
  * [26] J. Yoo, “General relativistic description of the observed galaxy power spectrum: Do we understand what we measure?,” Phys. Rev. D 82 no. 8, (Oct., 2010) 083508, arXiv:1009.3021.
  * [27] C. Bonvin and R. Durrer, “What galaxy surveys really measure,” Phys. Rev. D 84 no. 6, (Sept., 2011) 063505, arXiv:1105.5280.
  * [28] A. Challinor and A. Lewis, “Linear power spectrum of observed source number counts,” Phys. Rev. D 84 no. 4, (Aug., 2011) 043516, arXiv:1105.5292.
  * [29] D. Jeong, F. Schmidt, and C. M. Hirata, “Large-scale clustering of galaxies in general relativity,” Phys. Rev. D 85 (2012) 023504, arXiv:1107.5427.
  * [30] M. Bruni, R. Crittenden, K. Koyama, R. Maartens, C. Pitrou, and D. Wands, “Disentangling non-Gaussianity, bias and GR effects in the galaxy distribution,” Phys. Rev. D 85 (2012) 041301, arXiv:1106.3999.
  * [31] A. Kehagias, A. Moradinezhad Dizgah, J. Noreña, H. Perrier, and A. Riotto, “A consistency relation for the observed galaxy bispectrum and the local non-Gaussianity from relativistic corrections,” J. Cosmol. Astropart. Phys. 2015 no. 8, (Aug., 2015) 018, arXiv:1503.04467 [astro-ph.CO].
  * [32] E. Di Dio, H. Perrier, R. Durrer, G. Marozzi, A. Moradinezhad Dizgah, J. Noreña, and A. Riotto, “Non-Gaussianities due to relativistic corrections to the observed galaxy bispectrum,” J. Cosmol. Astropart. Phys. 2017 no. 3, (Mar., 2017) 006, arXiv:1611.03720 [astro-ph.CO].
  * [33] O. Umeh, S. Jolicoeur, R. Maartens, and C. Clarkson, “A general relativistic signature in the galaxy bispectrum: the local effects of observing on the lightcone,” J. Cosmol. Astropart. Phys. 2017 no. 3, (Mar., 2017) 034, arXiv:1610.03351 [astro-ph.CO].
  * [34] S. Jolicoeur, O. Umeh, R. Maartens, and C. Clarkson, “Imprints of local lightcone \projection effects on the galaxy bispectrum. Part II,” J. Cosmol. Astropart. Phys. 2017 no. 9, (Sept., 2017) 040, arXiv:1703.09630 [astro-ph.CO].
  * [35] K. Koyama, O. Umeh, R. Maartens, and D. Bertacca, “The observed galaxy bispectrum from single-field inflation in the squeezed limit,” J. Cosmol. Astropart. Phys. 2018 no. 7, (July, 2018) 050, arXiv:1805.09189 [astro-ph.CO].
  * [36] D. Bertacca, A. Raccanelli, N. Bartolo, M. Liguori, S. Matarrese, and L. Verde, “Relativistic wide-angle galaxy bispectrum on the light cone,” Phys. Rev. D 97 no. 2, (Jan., 2018) 023531, arXiv:1705.09306 [astro-ph.CO].
  * [37] O. Umeh, K. Koyama, R. Maartens, F. Schmidt, and C. Clarkson, “General relativistic effects in the galaxy bias at second order,” J. Cosmol. Astropart. Phys. 2019 no. 5, (May, 2019) 020, arXiv:1901.07460 [astro-ph.CO].
  * [38] C. Clarkson, E. M. de Weerd, S. Jolicoeur, R. Maartens, and O. Umeh, “The dipole of the galaxy bispectrum,” Mon. Not. R. Astron. Soc. 486 no. 1, (June, 2019) L101–L104, arXiv:1812.09512 [astro-ph.CO].
  * [39] O. Umeh and K. Koyama, “The galaxy bias at second order in general relativity with non-Gaussian initial conditions,” J. Cosmol. Astropart. Phys. 2019 no. 12, (Dec., 2019) 048, 1907.08094 [astro-ph.CO].
  * [40] R. Maartens, S. Jolicoeur, O. Umeh, E. M. De Weerd, and C. Clarkson, “Local primordial non-Gaussianity in the relativistic galaxy bispectrum,” J. Cosmol. Astropart. Phys. 2021 no. 4, (Apr., 2021) 013, arXiv:2011.13660 [astro-ph.CO].
  * [41] M. H. Goroff, B. Grinstein, S. Rey, and M. B. Wise, “Coupling of modes of cosmological mass density fluctuations,” Astrophys. J. 311 (Dec., 1986) 6–14.
  * [42] N. Makino, M. Sasaki, and Y. Suto, “Analytic approach to the perturbative expansion of nonlinear gravitational fluctuations in cosmological density and velocity fields,” Phys. Rev. D 46 (July, 1992) 585–602.
  * [43] B. Jain and E. Bertschinger, “Second-order power spectrum and nonlinear evolution at high redshift,” Astrophys. J. 431 (Aug., 1994) 495–505, arXiv:astro-ph/9311070.
  * [44] F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, “Large-scale structure of the Universe and cosmological perturbation theory,” Phys. Rep. 367 (Sept., 2002) 1–248, arXiv:astro-ph/0112551.
  * [45] J. M. Bardeen, “Gauge-invariant cosmological perturbations,” Phys. Rev. D 22 (Oct., 1980) 1882–1905.
  * [46] J. Yoo, “Proper-time hypersurface of nonrelativistic matter flows: Galaxy bias in general relativity,” Phys. Rev. D 90 no. 12, (Dec., 2014) 123507, arXiv:1408.5137.
  * [47] N. Kaiser, “On the spatial correlations of Abell clusters,” Astrophys. J. Lett. 284 (Sept., 1984) L9–L12.
  * [48] H. D. Politzer and M. B. Wise, “Relations between spatial correlations of rich clusters of galaxies,” Astrophys. J. Lett. 285 (Oct., 1984) L1–L3.
  * [49] J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser, “Excursion set mass functions for hierarchical Gaussian fluctuations,” Astrophys. J. 379 (Oct., 1991) 440–460.
  * [50] J. N. Fry, “The Evolution of Bias,” Astrophys. J. Lett. 461 (Apr., 1996) L65+.
  * [51] V. Desjacques, D. Jeong, and F. Schmidt, “Large-scale galaxy bias,” Phys. Rep. 733 (Feb., 2018) 1–193, arXiv:1611.09787 [astro-ph.CO].
  * [52] J. Yoo, N. Hamaus, U. Seljak, and M. Zaldarriaga, “Going beyond the Kaiser redshift-space distortion formula: A full general relativistic account of the effects and their detectability in galaxy clustering,” Phys. Rev. D 86 no. 6, (Sept., 2012) 063514, arXiv:1206.5809.
  * [53] P. Creminelli and M. Zaldarriaga, “A single-field consistency relation for the three-point function,” J. Cosmol. Astropart. Phys. 10 (Oct., 2004) 6, arXiv:astro-ph/0407059.
  * [54] P. Creminelli, J. Noreña, M. Simonović, and F. Vernizzi, “Single-field consistency relations of large scale structure,” J. Cosmol. Astropart. Phys. 2013 no. 12, (Dec., 2013) 025, arXiv:1309.3557 [astro-ph.CO].
  * [55] A. Kehagias and A. Riotto, “Symmetries and consistency relations in the large scale structure of the universe,” Nuclear Physics B 873 no. 3, (Aug., 2013) 514–529, arXiv:1302.0130 [astro-ph.CO].
  * [56] M. Peloso and M. Pietroni, “Galilean invariance and the consistency relation for the nonlinear squeezed bispectrum of large scale structure,” J. Cosmol. Astropart. Phys. 2013 no. 5, (May, 2013) 031, arXiv:1302.0223 [astro-ph.CO].
  * [57] L. Hui, A. Joyce, and S. S. C. Wong, “Inflationary soft theorems revisited: a generalized consistency relation,” J. Cosmol. Astropart. Phys. 2019 no. 2, (Feb., 2019) 060, arXiv:1811.05951 [hep-th].
  * [58] E. Mitsou and J. Yoo, “The spatial gauge-dependence of single-field inflationary bispectra,” Physics Letters B 828 (May, 2022) 137018, arXiv:2109.13154 [hep-th].
  * [59] S. Matarrese, F. Lucchin, and S. A. Bonometto, “A path-integral approach to large-scale matter distribution originated by non-Gaussian fluctuations,” Astrophys. J. Lett. 310 (Nov., 1986) L21–L26.
  * [60] S. Matarrese and L. Verde, “The Effect of Primordial Non-Gaussianity on Halo Bias,” Astrophys. J. Lett. 677 (Apr., 2008) L77–L80, arXiv:0801.4826.
  * [61] A. Slosar, C. Hirata, U. Seljak, S. Ho, and N. Padmanabhan, “Constraints on local primordial non-Gaussianity from large scale structure,” J. Cosmol. Astropart. Phys. 8 (Aug., 2008) 31, arXiv:0805.3580.
  * [62] W. H. Press and P. Schechter, “Formation of Galaxies and Clusters of Galaxies by Self-Similar Gravitational Condensation,” Astrophys. J. 187 (Feb., 1974) 425–438.
  * [63] J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, “The statistics of peaks of Gaussian random fields,” Astrophys. J. 304 (May, 1986) 15–61.
  * [64] J. Yoo and J.-O. Gong, “Relativistic effects and primordial non-Gaussianity in the matter density fluctuation,” Physics Letters B 754 (Mar., 2016) 94–98, arXiv:1509.08466 [astro-ph.CO].
  * [65] N. Grimm, F. Scaccabarozzi, J. Yoo, S. G. Biern, and J.-O. Gong, “Galaxy power spectrum in general relativity,” J. Cosmol. Astropart. Phys. 2020 no. 11, (Nov., 2020) 064, arXiv:2005.06484 [astro-ph.CO].
  * [66] H. A. Feldman, N. Kaiser, and J. A. Peacock, “Power-spectrum analysis of three-dimensional redshift surveys,” Astrophys. J. 426 (May, 1994) 23–37, arXiv:astro-ph/9304022.
  * [67] M. Tegmark, A. J. S. Hamilton, M. A. Strauss, M. S. Vogeley, and A. S. Szalay, “Measuring the Galaxy Power Spectrum with Future Redshift Surveys,” Astrophys. J. 499 (May, 1998) 555, arXiv:9708020.
  * [68] A. Raccanelli, L. Samushia, and W. J. Percival, “Simulating redshift-space distortions for galaxy pairs with wide angular separation,” Mon. Not. R. Astron. Soc. 409 (Dec., 2010) 1525–1533, arXiv:1006.1652.
  * [69] J. Yoo and V. Desjacques, “All-sky analysis of the general relativistic galaxy power spectrum,” Phys. Rev. D 88 no. 2, (July, 2013) 023502, arXiv:1301.4501.
  * [70] J. Yoo and U. Seljak, “Joint analysis of gravitational lensing, clustering, and abundance: Toward the unification of large-scale structure analysis,” Phys. Rev. D 86 no. 8, (Oct., 2012) 083504, arXiv:1207.2471.
  * [71] G. F. R. Ellis and M. Bruni, “Covariant and gauge-invariant approach to cosmological density fluctuations,” Phys. Rev. D 40 no. 6, (Sept., 1989) 1804–1818.
  * [72] M. Bruni and S. Sonego, “Observables and gauge invariance in the theory of nonlinear spacetime perturbations,” Classical and Quantum Gravity 16 (July, 1999) L29–L36, gr-qc/9906017.
  * [73] J. Yoo and R. Durrer, “Gauge-transformation properties of cosmological observables and its application to the light-cone average,” J. Cosmol. Astropart. Phys. 9 (Sept., 2017) 016, arXiv:1705.05839.
  * [74] T. Buchert and J. Ehlers, “Averaging inhomogeneous Newtonian cosmologies.,” Astron. Astrophys. 320 (Apr., 1997) 1–7, astro-ph/9510056.
  * [75] V. F. Mukhanov, L. R. W. Abramo, and R. H. Brandenberger, “Backreaction Problem for Cosmological Perturbations,” Physical Review Letters 78 (Mar., 1997) 1624–1627.
  * [76] C. M. Hirata and U. Seljak, “Can superhorizon cosmological perturbations explain the acceleration of the universe?,” Phys. Rev. D 72 no. 8, (Oct., 2005) 083501, astro-ph/0503582.
  * [77] J. Larena, “Spatially averaged cosmology in an arbitrary coordinate system,” Phys. Rev. D 79 no. 8, (Apr., 2009) 084006, arXiv:0902.3159.
  * [78] A. Ishibashi and R. M. Wald, “Can the acceleration of our universe be explained by the effects of inhomogeneities?,” Classical and Quantum Gravity 23 (Jan., 2006) 235–250, gr-qc/0509108.
  * [79] C. Clarkson, G. Ellis, J. Larena, and O. Umeh, “Does the growth of structure affect our dynamical models of the universe? The averaging, backreaction and fitting problems in cosmology,” Reports on Progress in Physics 74 (2011) 112901, arXiv:1109.2314.
  * [80] J. Yoo, E. Mitsou, N. Grimm, R. Durrer, and A. Refregier, “Cosmological information contents on the light-cone,” J. Cosmol. Astropart. Phys. 2019 no. 12, (Dec, 2019) , arXiv:1905.08262.
  * [81] J. Yoo, E. Mitsou, Y. Dirian, and R. Durrer, “$\bar{T}$: A New Cosmological Parameter?,” Phys. Rev. D 100 (Sept., 2019) 063510, arXiv:1905.09288.
  * [82] S. Baumgartner and J. Yoo, “Monopole Fluctuation of the CMB and its Gauge Invariance,” accepted for publication in PRD (Dec., 2021) , arXiv:2012.03968 [astro-ph.CO].
  * [83] E. Mitsou, J. Yoo, R. Durrer, F. Scaccabarozzi, and V. Tansella, “General and consistent statistics for cosmological observations,” Physical Review Research 2 no. 3, (July, 2020) 033004, arXiv:1905.01293 [astro-ph.CO].
  * [84] S. Weinberg, “Adiabatic modes in cosmology,” Phys. Rev. D 67 no. 12, (June, 2003) 123504, astro-ph/0302326.
  * [85] P. Creminelli, J. Noreña, and M. Simonović, “Conformal consistency relations for single-field inflation,” J. Cosmol. Astropart. Phys. 2012 no. 7, (July, 2012) 052, arXiv:1203.4595 [hep-th].
  * [86] K. Hinterbichler, L. Hui, and J. Khoury, “Conformal symmetries of adiabatic modes in cosmology,” J. Cosmol. Astropart. Phys. 2012 no. 8, (Aug., 2012) 017, arXiv:1203.6351 [hep-th].
  * [87] C. Cheung, A. L. Fitzpatrick, J. Kaplan, L. Senatore, and P. Creminelli, “The effective field theory of inflation,” Journal of High Energy Physics 2008 no. 3, (Mar., 2008) 014–014, arXiv:0709.0293 [hep-th].
  * [88] K. Hinterbichler, L. Hui, and J. Khoury, “An infinite set of Ward identities for adiabatic modes in cosmology,” J. Cosmol. Astropart. Phys. 2014 no. 1, (Jan., 2014) 039, arXiv:1304.5527 [hep-th].
  * [89] L. Dai, E. Pajer, and F. Schmidt, “Conformal Fermi Coordinates,” ArXiv e-prints (Feb., 2015) , arXiv:1502.02011.
  * [90] G. Cabass, E. Pajer, and F. Schmidt, “How Gaussian can our Universe be?,” J. Cosmol. Astropart. Phys. 2017 no. 1, (Jan., 2017) 003, arXiv:1612.00033 [hep-th].
  * [91] E. Mitsou and J. Yoo, “Tetrad formalism for exact cosmological observables,” Springer Briefs in Physics (Aug., 2020) 1–139, arXiv:1908.10757 [astro-ph.CO].
  * [92] E. Mitsou and J. Yoo, “Large gauge transformations, local coordinates and cosmological observables,” arXiv e-prints (May, 2022) , arXiv:2205.01544 [astro-ph.CO].
  * [93] J. Yoo, “Complete Treatment of Galaxy Two-Point Statistics: Gravitational Lensing Effects and Redshift-Space Distortions,” Phys. Rev. D 79 no. 2, (Jan., 2009) 023517, arXiv:0808.3138.
  * [94] J. N. Fry and E. Gaztanaga, “Biasing and hierarchical statistics in large-scale structure,” Astrophys. J. 413 (Aug., 1993) 447–452, astro-ph/9302009.
  * [95] R. J. Scherrer and D. H. Weinberg, “Constraints on the Effects of Locally Biased Galaxy Formation,” Astrophys. J. 504 (Sept., 1998) 607, astro-ph/9712192.
  * [96] R. K. Sheth and G. Tormen, “Large-scale bias and the peak background split,” Mon. Not. R. Astron. Soc. 308 (Sept., 1999) 119–126, astro-ph/9901122.
  * [97] T. Baldauf, U. Seljak, V. Desjacques, and P. McDonald, “Evidence for quadratic tidal tensor bias from the halo bispectrum,” Phys. Rev. D 86 no. 8, (Oct., 2012) 083540, arXiv:1201.4827.
  * [98] N. Dalal, U.-L. Pen, and U. Seljak, “Large-scale BAO signatures of the smallest galaxies,” J. Cosmol. Astropart. Phys. 11 (Nov., 2010) 7, arXiv:1009.4704.
  * [99] J. Yoo, N. Dalal, and U. Seljak, “Supersonic relative velocity effect on the baryonic acoustic oscillation measurements,” J. Cosmol. Astropart. Phys. 7 (July, 2011) 18, arXiv:1105.3732.
  * [100] J. Yoo and U. Seljak, “Signatures of first stars in galaxy surveys: Multitracer analysis of the supersonic relative velocity effect and the constraints from the BOSS power spectrum measurements,” Phys. Rev. D 88 no. 10, (Nov., 2013) 103520, arXiv:1308.1401.
  * [101] F. Schmidt, “Effect of relative velocity and density perturbations between baryons and dark matter on the clustering of galaxies,” Phys. Rev. D 94 no. 6, (Sept., 2016) 063508, arXiv:1602.09059 [astro-ph.CO].
  * [102] P. McDonald, “Primordial non-Gaussianity: Large-scale structure signature in the perturbative bias model,” Phys. Rev. D 78 no. 12, (Dec., 2008) 123519, arXiv:0806.1061.
  * [103] P. McDonald and A. Roy, “Clustering of dark matter tracers: generalizing bias for the coming era of precision LSS,” J. Cosmol. Astropart. Phys. 8 (Aug., 2009) 20, arXiv:0902.0991.
  * [104] R. A. Porto, L. Senatore, and M. Zaldarriaga, “The Lagrangian-space Effective Field Theory of large scale structures,” J. Cosmol. Astropart. Phys. 2014 no. 5, (May, 2014) 022, arXiv:1311.2168 [astro-ph.CO].
  * [105] L. Senatore, “Bias in the effective field theory of large scale structures,” J. Cosmol. Astropart. Phys. 2015 no. 11, (Nov., 2015) 007, arXiv:1406.7843 [astro-ph.CO].
  * [106] Z. Vlah, M. White, and A. Aviles, “A Lagrangian effective field theory,” J. Cosmol. Astropart. Phys. 2015 no. 9, (Sept., 2015) 014, arXiv:1506.05264 [astro-ph.CO].
  * [107] F. Schmidt, F. Elsner, J. Jasche, N. M. Nguyen, and G. Lavaux, “A rigorous EFT-based forward model for large-scale structure,” J. Cosmol. Astropart. Phys. 2019 no. 1, (Jan., 2019) 042, arXiv:1808.02002 [astro-ph.CO].
  * [108] Z. Vlah, N. E. Chisari, and F. Schmidt, “An EFT description of galaxy intrinsic alignments,” J. Cosmol. Astropart. Phys. 2020 no. 1, (Jan., 2020) 025, arXiv:1910.08085 [astro-ph.CO].
  * [109] D. J. Eisenstein, J. Annis, J. E. Gunn, A. S. Szalay, A. J. Connolly, R. C. Nichol, N. A. Bahcall, M. Bernardi, S. Burles, F. J. Castander, M. Fukugita, et al., “Spectroscopic Target Selection for the Sloan Digital Sky Survey: The Luminous Red Galaxy Sample,” Astron. J. 122 (Nov., 2001) 2267–2280, arXiv:astro-ph/0108153.
  * [110] M. White, M. Blanton, A. Bolton, D. Schlegel, J. Tinker, A. Berlind, L. da Costa, E. Kazin, Y.-T. Lin, M. Maia, C. K. McBride, N. Padmanabhan, J. Parejko, W. Percival, F. Prada, B. Ramos, E. Sheldon, F. de Simoni, R. Skibba, D. Thomas, D. Wake, I. Zehavi, Z. Zheng, R. Nichol, D. P. Schneider, M. A. Strauss, B. A. Weaver, and D. H. Weinberg, “The Clustering of Massive Galaxies at z ~ 0.5 from the First Semester of BOSS Data,” Astrophys. J. 728 (Feb., 2011) 126, arXiv:1010.4915.
  * [111] J. Yoo and M. Zaldarriaga, “Beyond the Linear-Order Relativistic Effect in Galaxy Clustering: Second-Order Gauge-Invariant Formalism,” Phys. Rev. D 90 no. 10, (July, 2014) 023513, arXiv:1406.4140.
  * [112] E. Di Dio, R. Durrer, G. Marozzi, and F. Montanari, “Galaxy number counts to second order and their bispectrum,” J. Cosmol. Astropart. Phys. 12 (Dec., 2014) 17, arXiv:1407.0376.
  * [113] D. Bertacca, R. Maartens, and C. Clarkson, “Observed galaxy number counts on the lightcone up to second order: I. Main result,” J. Cosmol. Astropart. Phys. 9 (Sept., 2014) 037, arXiv:1405.4403.
  * [114] D. Bertacca, R. Maartens, and C. Clarkson, “Observed galaxy number counts on the lightcone up to second order: II. Derivation,” J. Cosmol. Astropart. Phys. 11 (Nov., 2014) 13, arXiv:1406.0319.
  * [115] J.-C. Hwang and H. Noh, “Second-order perturbations of a zero-pressure cosmological medium: Proofs of the relativistic-Newtonian correspondence,” Phys. Rev. D 72 no. 4, (Aug., 2005) 044011, arXiv/0412128.
  * [116] J.-C. Hwang and H. Noh, “Third-order perturbations of a zero-pressure cosmological medium: Pure general relativistic nonlinear effects,” Phys. Rev. D 72 no. 4, (Aug., 2005) 044012, gr-qc/0412129.
  * [117] J. Hwang and H. Noh, “Why Newtonian gravity is reliable in large-scale cosmological simulations,” Mon. Not. R. Astron. Soc. 367 (Apr., 2006) 1515–1520, arXiv:0507159.
  * [118] F. R. Bouchet, S. Colombi, E. Hivon, and R. Juszkiewicz, “Perturbative Lagrangian approach to gravitational instability.,” Astron. Astrophys. 296 (Apr., 1995) 575, arXiv:astro-ph/9406013 [astro-ph].
  * [119] D. Bertacca, N. Bartolo, M. Bruni, K. Koyama, R. Maartens, S. Matarrese, M. Sasaki, and D. Wands, “Galaxy bias and gauges at second order in general relativity,” Classical and Quantum Gravity 32 no. 17, (Sept., 2015) 175019, arXiv:1501.03163 [astro-ph.CO].
  * [120] J. E. Gunn, “On the Propagation of Light in Inhomogeneous Cosmologies. I. Mean Effects,” Astrophys. J. 150 (Dec., 1967) 737.
  * [121] R. K. Sachs and A. M. Wolfe, “Perturbations of a Cosmological Model and Angular Variations of the Microwave Background,” Astrophys. J. 147 (Jan., 1967) 73+.
  * [122] N. Kaiser, “Clustering in real space and in redshift space,” Mon. Not. R. Astron. Soc. 227 (July, 1987) 1–21.
  * [123] J. Miralda-Escudé, “The correlation function of galaxy ellipticities produced by gravitational lensing,” Astrophys. J. 380 (Oct., 1991) 1–8.
  * [124] N. Kaiser, “Weak gravitational lensing of distant galaxies,” Astrophys. J. 388 (Apr., 1992) 272–286.
  * [125] R. Narayan, “Gravitational lensing and quasar-galaxy correlations,” Astrophys. J. Lett. 339 (Apr., 1989) L53–L56.
  * [126] M. Bartelmann, “Cosmological parameters from angular correlations between QSOs and galaxies.,” Astron. Astrophys. 298 (June, 1995) 661, arXiv:astro-ph/9407091.
  * [127] B. Jain, “Magnification Effects as Measures of Large-Scale Structure,” Astrophys. J. Lett. 580 (Nov., 2002) L3–L6, arXiv:astro-ph/0208515.
  * [128] J. Yoo, “Relativistic effect in galaxy clustering,” Class. Quant. Grav. 31 no. 23, (Sept., 2014) 234001, arXiv:1409.3223.
  * [129] H. J. Mo and S. D. M. White, “An analytic model for the spatial clustering of dark matter haloes,” Mon. Not. R. Astron. Soc. 282 (Sept., 1996) 347–361, arXiv:9512127.
  * [130] P. Catelan, F. Lucchin, S. Matarrese, and C. Porciani, “The bias field of dark matter haloes,” Mon. Not. R. Astron. Soc. 297 (July, 1998) 692–712, astro-ph/9708067.
  * [131] P. Catelan, C. Porciani, and M. Kamionkowski, “Two ways of biasing galaxy formation,” Mon. Not. R. Astron. Soc. 318 (Nov., 2000) L39–L44, astro-ph/0005544.
  * [132] R. K. Sheth, H. J. Mo, and G. Tormen, “Ellipsoidal collapse and an improved model for the number and spatial distribution of dark matter haloes,” Mon. Not. R. Astron. Soc. 323 (May, 2001) 1–12, astro-ph/9907024.
  * [133] F. Schmidt, D. Jeong, and V. Desjacques, “Peak-background split, renormalization, and galaxy clustering,” Phys. Rev. D 88 no. 2, (July, 2013) 023515, arXiv:1212.0868.
  * [134] R. Scoccimarro, E. Sefusatti, and M. Zaldarriaga, “Probing primordial non-Gaussianity with large-scale structure,” Phys. Rev. D 69 no. 10, (May, 2004) 103513, arXiv:astro-ph/0312286.
  * [135] E. Sefusatti and E. Komatsu, “Bispectrum of galaxies from high-redshift galaxy surveys: Primordial non-Gaussianity and nonlinear galaxy bias,” Phys. Rev. D 76 no. 8, (Oct., 2007) 083004, arXiv:0705.0343.
  * [136] J. Yoo, N. Grimm, E. Mitsou, A. Amara, and A. Refregier, “Gauge-invariant formalism of cosmological weak lensing,” J. Cosmol. Astropart. Phys. 4 (Apr., 2018) 029, arXiv:1802.03403.
  * [137] E. Komatsu and D. N. Spergel, “Acoustic signatures in the primary microwave background bispectrum,” Phys. Rev. D 63 no. 6, (Mar., 2001) 063002, arXiv:astro-ph/0005036.
  * [138] J. Yoo and F. Scaccabarozzi, “Unified treatment of the luminosity distance in cosmology,” J. Cosmol. Astropart. Phys. 9 (Sept., 2016) 046, 1606.08453.
  * [139] S. G. Biern and J. Yoo, “Correlation function of the luminosity distances,” J. Cosmol. Astropart. Phys. (Sept., 2017) 026, 1704.07380.
  * [140] F. Scaccabarozzi, J. Yoo, and S. G. Biern, “Galaxy two-point correlation function in general relativity,” J. Cosmol. Astropart. Phys. 10 (Oct., 2018) 024, arXiv:1807.09796.
  * [141] M. Magi and J. Yoo, “Second-Order Gauge-Invariant Formalism for the Cosmological Observables: Complete Verification of their Gauge-Invariance,” submitted to J. Cosmol. Astropart. Phys. (Apr., 2022) , arXiv:2204.01751.
  * [142] J. Binney and T. Quinn, “Gaussian random fields in spherical coordinates,” Mon. Not. R. Astron. Soc. 249 (Apr., 1991) 678–683.
  * [143] K. B. Fisher, O. Lahav, Y. Hoffman, D. Lynden-Bell, and S. Zaroubi, “Wiener reconstruction of density, velocity and potential fields from all-sky galaxy redshift surveys,” Mon. Not. R. Astron. Soc. 272 (Feb., 1995) 885–908, arXiv:9406009.
  * [144] A. F. Heavens and A. N. Taylor, “A spherical harmonic analysis of redshift space,” Mon. Not. R. Astron. Soc. 275 (July, 1995) 483–497, arXiv:9409027.
  * [145] A. Rassat and A. Refregier, “3D spherical analysis of baryon acoustic oscillations,” Astron. Astrophys. 540 (Apr., 2012) A115, arXiv:1112.3100.
  * [146] H. Noh and J.-C. Hwang, “Second-order perturbations of the Friedmann world model,” Phys. Rev. D 69 no. 10, (May, 2004) 104011, arXiv:0305123.
  * [147] C.-P. Ma and E. Bertschinger, “Cosmological Perturbation Theory in the Synchronous and Conformal Newtonian Gauges,” Astrophys. J. 455 (Dec., 1995) 7, arXiv:astro-ph/9401007.
  * [148] P. J. E. Peebles, The large-scale structure of the universe.  Princeton University Press, Princeton, 1980.