#  Gravitational wave fossils in nonlinear regime: halo tidal bias and
intrinsic alignments from gravitational wave separate universe simulations

Kazuyuki Akitsu<EMAIL_ADDRESS>School of Natural Sciences, Institute for
Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA  Yin Li  Center for
Computational Astrophysics & Center for Computational Mathematics, Flatiron
Institute, 162 5th Avenue, New York, NY 10010, USA  Teppei Okumura  Institute
of Astronomy and Astrophysics, Academia Sinica, No. 1, Sec. 4, Roosevelt Road,
Taipei 10617, Taiwan  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), UTIAS, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa,
Chiba 277-8583, Japan

###### Abstract

We investigate impacts of long-wavelength gravitational waves (GWs) on
nonlinear structure formation by utilizing the tidal separate universe
simulations. Based on the equivalence of a long-wavelength GW to a uniform
tidal field in a local frame, we provide a way to incorporate a long-
wavelength GW into the tidal separate universe simulation as an effective
anisotropic expansion. This methodology enables us to study effects of GWs on
large-scale structure efficiently. We measure the anisotropic imprint in the
local power spectrum from the tidal separate universe simulations with GWs,
which corresponds to the scalar-scalar-tensor bispectrum in squeezed limit or
the so-called power spectrum response to GWs. We also detect the halo tidal
bias induced by GWs from the response of the halo-matter cross-power spectrum
to GWs, as well as the linear shape bias (or the linear alignment coefficient)
induced by GWs from the one-point function of the halo ellipticity. In
contrast to the case of the tidal field induced by scalar perturbations, we
discover that the wavenumber dependence of the temporal evolution of GWs
naturally causes these biases to be scale-dependent. We also find that this
scale dependence is well approximated by the second-order density induced by
the coupling between scalar and tensor perturbation. This highlights that the
structure formation, especially the process to determine the halo shape, is
nonlocal in time. Our findings lay the foundation for predicting the impact of
GWs on large-scale structure.

##  I Introduction

Gravitational waves (GWs) serve as vital means to observe the universe. In
particular, since long-wavelength GWs in the Gpc-Mpc range are unlikely to
originate from astrophysical events, they are thought to have a cosmological
origin (e.g., inflation Starobinsky (1979); Sato (1981); Guth (1981); Linde
(1982); Maggiore (2000)), making them an important probe in cosmology. In
spite of the various experiments trying to detect the primordial $B$-mode
signal in the polarization of the Cosmic Microwave Background (CMB) Komatsu
_et al._ (2011); Akrami  _et al._ (2020); Ade  _et al._ (2018); Adachi  _et
al._ (2022); Abazajian  _et al._ (2016), which is one of the most powerful
methods to hunt for long-wavelength cosmological GWs Zaldarriaga and Seljak
(1997); Kamionkowski  _et al._ (1997), such GWs have not been observed yet.

Compared to the numerous studies for the effect of GWs on the perturbations of
the CMB, that on large-scale structure (LSS) of the universe has not received
as much attention. As large-scale structure is dominantly sourced by scalar
perturbations, there is a long history for the study on how the scalar
perturbations have shaped LSS, including the linear and nonlinear perturbation
theory Dodelson (2003); Suto and Sasaki (1991); Makino  _et al._ (1992);
Bernardeau  _et al._ (2002); Baumann  _et al._ (2012) and the $N$-body
simulation Springel (2005); Potter  _et al._ (2016); Springel  _et al._
(2021); Angulo and Hahn (2021). On the other hand, there are fewer studies of
the effect of the tensor perturbations (GWs) on nonlinear structure formation,
as we list below.

There are two types of effects of GWs on LSS observables: the dynamical effect
and the projection effect. The former refers to the effect of GWs on nonlinear
structure formation itself, whereas the latter refers to the effect of GWs on
the light path emitted from distant galaxies to us. In other words, the
dynamical effect would be observed by a comoving observer in a local frame,
while the projection effect comes from the mapping of observables from the
galaxy’s local frame to our local frame (at the earth), which includes the
Sachs-Wolfe effect and the gravitational lensing effect caused by GWs. Both
effects have been studied by means of the perturbation theory; Refs. Kaiser
and Jaffe (1997); Dodelson  _et al._ (2003); Yoo  _et al._ (2009); Jeong  _et
al._ (2012); Schmidt and Jeong (2012a); Jeong and Schmidt (2012) formulated
the projection effect on the galaxy clustering and the galaxy shape (shear) by
GWs, while Ref Schmidt and Jeong (2012a) pointed out that long-wavelength GWs
can also contribute the intrinsic alignments of galaxy shapes and Refs. Masui
and Pen (2010); Dai  _et al._ (2013); Schmidt  _et al._ (2014) computed the
second order matter density contrast induced by the coupling between scalar
perturbations and long-wavelength tensor perturbations (GWs). For the
projection effect the perturbative treatment would be adequate. For the
dynamical effect, however, the perturbation theory breaks down in nonlinear
scales and thus the nonlinear nature of LSS requires $N$-body simulations with
GWs to capture the fully nonlinear impact of GWs on structure formation.
Furthermore, given that the biased tracers of LSS such as halos are themselves
nonlinear objects, such simulations are necessary to understand their biases
to GWs even in the linear regime.

An $N$-body simulation with GWs is generally challenging because the usual
$N$-body simulations are based on Newtonian gravitational dynamics in an
expanding background whereas GWs are a purely general-relativistic effect. The
most straightforward way to introduce GWs into $N$-body simulations is to
develop a general-relativistic cosmological simulation Adamek  _et al._ (2013,
2016a, 2016b), though Refs. Adamek  _et al._ (2016a, b) only considered the
second-order (induced) tensor perturbations by scalar perturbations.111
Technically the induced tensor perturbations computed in Refs. Adamek  _et
al._ (2016a, b) involve non-GWs contributions, which are tensor modes but not
propagating waves. See e.g., Ref. Domènech and Sasaki (2021) for detailed
discussion on this issue. Also, since GWs are much smaller than scalar
perturbations, it is difficult to single out the effect of GWs on structure
formation from that of the scalar perturbations in this sort of simulations.

In this paper, we utilize a separate universe approach to circumvent this
issue on $N$-body simulations with GWs. In the separate universe approach, the
influence of a long-wavelength perturbation is absorbed into the cosmic
expansion observed in the local frame, thereby the local expansion becomes
different from the global one. Accordingly, nonlinear structure formation in
the local region responds to this difference in the background expansions.
Using this technique, the response to the long-wavelength perturbation can be
accurately measured in $N$-body simulations Sirko (2005); Gnedin  _et al._
(2011); Baldauf  _et al._ (2011); Li  _et al._ (2014); Wagner  _et al._
(2015); Lazeyras  _et al._ (2016); Li  _et al._ (2016); Baldauf  _et al._
(2016). Recently, the anisotropic extension of the separate universe
simulation was developed in Refs. Stücker  _et al._ (2021); Masaki  _et al._
(2020); Akitsu  _et al._ (2021a). They considered a long-wavelength tidal
perturbation sourced by the long-wavelength scalar tides
$\propto\left(\partial_{i}\partial_{j}-\frac{1}{3}\delta^{\textsc{k}}_{ij}\right)\Phi$
Dai  _et al._ (2015a); Akitsu  _et al._ (2017); Akitsu and Takada (2018);
Akitsu  _et al._ (2019); Li  _et al._ (2018); Taruya and Akitsu (2021), and
the perturbation was absorbed into the local background, making the local
cosmic expansion anisotropic.

Taking advantage of the equivalence of a long-wavelength GW to a long-
wavelength tidal field in a local region, we apply this tidal separate
universe simulation to measure the impact of GWs on structure formation. In
contrast to the scalar case, the time evolution of GWs is scale-dependent (or
wavenumber-dependent). As a result, when mimicking the long-wavelength GW as
the local anisotropic expansion, the anisotropic expansion rate would be
different depending on the wavenumber of GWs, and so is the response of large-
scale structure, as considered for isotropic scale-dependent long-wavelength
perturbation in Refs. Hu  _et al._ (2016); Chiang  _et al._ (2016, 2018). The
purpose of this paper is to generalize the tidal separate universe simulation
for GWs in order to study the scale-dependent responses for long-wavelength
GWs of different wavelengths.

The remainder of this paper provides a way to implement long-wavelength GWs
into the tidal separate universe simulation and presents imprints of GWs on
large-scale structure measured from newly developed simulations. After giving
a brief review about the local coordinates in the presence of GWs and
perturbative results in Sec. II, we construct the tidal separate universe with
GWs in Sec. III. In Sec. IV, we measure the power spectrum responses for
matter auto-, halo-matter cross-, and halo auto-power spectra, which is
related to the scalar-scalar-tensor bispectrum in squeezed limit. Sec. V and
Sec. VI are devoted to the measurements of the halo tidal bias and linear
alignment coefficient (or the linear shape bias) from GWs, respectively. We
discuss possible observables in large-scale structure to probe GWs in Sec.
VII. Throughout this paper we adopt cosmological parameters consistent with
Planck result: $\Omega_{\rm r0}=4.1577\times 10^{-5}$, $\Omega_{\rm
m0}=0.3089$, $\Omega_{\Lambda 0}=0.6911$, $H_{0}=67.74$ Ade  _et al._ (2016).

##  II Local frame and the perturbative results

Here we first introduce the local coordinates in the presence of long-
wavelength GWs. We then briefly summarize the derivation of the second-order
density perturbations induced by the coupling between scalar and tensor
perturbations (GWs), following Ref. Schmidt  _et al._ (2014). We employ the
Lagrangian perturbation formalism, which can be straightforwardly used to
construct the tidal separate universe in the next section.

###  II.1 Long-wavelength gravitational waves in the conformal Fermi
coordinates

Figure 1:  Left panels from top to bottom: Transfer function of GWs
$\mathcal{T}(a;k_{L})$, growth coefficient $\alpha(a;k_{L})$ and dilation
coefficient $\beta(a;k_{L})$ in Eq. (29) as a function of scale factor $a$ for
various wavenumbers $k_{L}$. Right panels from top to bottom: Same as the left
panels but as a function of wavenumber of GWs, $k_{L}$, at various redshifts.
Note that the functional shape of $\beta(a;k_{L})$ is equivalent to that of
$\mathcal{T}(a;k_{L})$ because
$\beta(a;k_{L})=-\frac{1}{2}\left[\mathcal{T}(a;k_{L})-1\right]$.

In the cosmological context, GWs are defined as the trace-free transverse
components in the perturbed FLRW metric,

$\displaystyle
ds^{2}=a^{2}[-d\eta^{2}+(\delta^{\textsc{k}}_{ij}+h_{ij})dx^{i}dx^{j}],$ (1)

where $a$ is the scale factor, $\eta$ is the conformal time,
$\delta^{\textsc{k}}_{ij}$ is Kronecker’s delta, and $h_{ij}$ is GWs
satisfying $\tensor{h}{{}^{i}_{i}}=0$ and $\partial^{i}h_{ij}=0$. We use
$k_{L}$ to denote the wavenumber of GWs in what follows since we focus on the
long-wavelength GWs in this paper. We introduce the transfer function of GWs,
$\mathcal{T}(\eta;k_{L})$, through

$\displaystyle h_{ij}(\eta;{\bf k}_{L})=\mathcal{T}(\eta;k_{L})h^{\rm
ini}_{ij}({\bf k}_{L}),$ (2)

where $h_{ij}^{\rm ini}(k_{L})$ denotes the primordial value of GWs, i.e.,
$h_{ij}^{\rm ini}(k_{L})=h_{ij}(0;k_{L})$. The transfer function
$\mathcal{T}(\eta;k_{L})$ obeys

$\displaystyle\mathcal{T}^{\prime\prime}(\eta;k_{L})+2\mathcal{H}\mathcal{T}^{\prime}(\eta;k_{L})+k_{L}^{2}\mathcal{T}(\eta;k_{L})=0,$
(3)

with ${}^{\prime}=\differential/\differential\eta$, and $\mathcal{H}=aH$ being
the conformal Hubble paramter. The top-left panel of Fig. 1 shows
$\mathcal{T}(\eta;k_{L})$ as a function of the scale factor for various
wavenumbers of GWs. Similarly, the top-right panel shows that as a function of
the wavenumbers of GWs at various redshifts. For any $k_{L}$,
$\mathcal{T}(\eta;k_{L})$ remains unity when $k_{L}\eta\ll 1$, which means
that GWs are frozen before they enter the horizon.

For later convenience, we expand GWs by the helicity basis:

$\displaystyle h_{ij}(\eta;{\bf k}_{L})=\sum_{\lambda=\pm
2}h_{(\lambda)}(\eta;{\bf k}_{L})e^{(\lambda)}_{ij}(\hat{k}_{L}),$ (4)

where $e_{ij}^{(\pm 2)}\equiv e_{i}^{(\pm)}e_{j}^{(\pm)}$ and ${\bf
e}^{(\pm)}\equiv({\bf e}_{1}\mp i{\bf e}_{2})/\sqrt{2}$ with $\\{{\bf
k}_{L},{\bf e}_{1},{\bf e}_{2}\\}$ being an orthonormal set. The helicity
basis satisfies $e^{(\lambda)}_{ij}e^{ij}_{(-\lambda)}=1$,
$e^{(\lambda)}_{ij}e^{ij}_{(\lambda)}=0$, and
$h_{(\lambda)}=e^{ij}_{(-\lambda)}h_{ij}$. The power spectrum of GWs for each
helicity mode is defined as

$\displaystyle\langle h_{(\lambda)}({\bf k}_{L};\eta)h^{*}_{(\lambda)}({\bf
k}^{\prime}_{L};\eta^{\prime})\rangle=(2\pi)^{3}\delta_{\rm D}^{(3)}({\bf
k}_{L}-{\bf k}^{\prime}_{L})P_{h_{(\lambda)}}(k_{L};\eta,\eta^{\prime}),$ (5)

which is related to the primordial power spectrum $P_{h_{(\lambda)}}(k_{L})$
through

$\displaystyle
P_{h_{(\lambda)}}(k_{L};\eta,\eta^{\prime})=\mathcal{T}(\eta;k_{L})\mathcal{T}(\eta^{\prime};k_{L})P_{h_{(\lambda)}}(k_{L}).$
(6)

For unpolarized GWs, the two power spectra has equal power, i.e.,
$P_{h_{(+2)}}(k)=P_{h_{(-2)}}(k)\equiv P_{h}(k)/2$ with $P_{h}(k)$ being the
total power spectrum of GWs. For chiral GWs, the chiral parameter $\chi(k)$
defined as

$\displaystyle\chi(k)\equiv\frac{P_{(+2)}(k)-P_{(-2)}(k)}{P_{h}(k)},$ (7)

measures the degree of the parity-breaking in GWs. The total power spectrum of
GWs is often characterized via

$\displaystyle\frac{k^{3}P_{h}(k)}{2\pi^{2}}\equiv
rA_{s}\left(\frac{k}{k_{*}}\right)^{n_{T}}$ (8)

with $r$ being the scalar-tensor ratio, $A_{s}$ being the amplitude for the
primordial curvature perturbations at the pivot scale $k_{*}$, and $n_{T}$
being the tensor tilt.

To investigate the physical effects of long-wavelength GWs on scalar
perturbations at smaller scales, consider a local region centered at the
timelike geodesic of a comoving observer. In this local patch we can construct
the so-called conformal Fermi coordinates (CFC), which is an extension of the
Fermi normal coordinates Manasse and Misner (1963), developed in Refs. Schmidt
_et al._ (2014); Dai  _et al._ (2015b). The metric of this coordinates,
$g_{\mu\nu}^{F}$, takes the form of the FLRW metric along the central geodesic
with leading-order corrections of ${\cal O}(x^{2}_{F})$. We are interested in
the interaction between long-wavelength GWs and non-relativistic matter. In
this case $g_{00}^{F}$ encodes all the relevant impact from the long-
wavelength perturbations because the dynamics of non-relativistic matter is
solely determined by the usual Newtonian potential in $g_{00}^{F}$. One can
show that given the global coordinates Eq. (1), $g_{00}^{F}$ is computed as

$\displaystyle g_{00}^{F}=-a^{2}\left[1+\tau_{ij}x_{F}^{i}x_{F}^{j}\right],$
(9)

with

$\displaystyle\tau_{ij}(\eta;k_{L})=-\frac{1}{2}\left[a^{-1}\left(ah_{ij}^{\prime}\right)^{\prime}\right]=-\frac{1}{2}\left(h_{ij}^{\prime\prime}+\mathcal{H}h_{ij}^{\prime}\right)\equiv
T(\eta;k_{L})h^{\rm ini}_{ij}$ (10)

representing the effective tidal field induced by the long-wavelength GWs in
the local region. The derivation is summarized in App. A. It is worth noting
that $\tau_{ij}$ takes effect only after GWs cross the horizon as implied by
the time derivative on $h_{ij}$. Note that $x_{F}^{i}$ corresponds to the
comoving distance, unlike the physical (proper) distance in the Fermi normal
coordinates. In what follows, we work in this frame and drop the subscript $F$
in $x_{F}^{i}$, namely, denote $x^{i}$ as $x_{F}^{i}$.

###  II.2 Second order density induced by the interaction between GWs and
scalar perturbations

The equation of motion of a matter particle in the local frame is 222Strictly
speaking, this is justified by considering the geodesic equation in the local
frame.

$\displaystyle\frac{\differential^{2}r_{i}}{\differential
t^{2}}=\frac{1}{a^{2}}\left[\frac{\differential^{2}}{\differential\eta^{2}}-\mathcal{H}\frac{\differential}{\differential\eta}\right]r_{i}=-\frac{\partial}{\partial
r_{i}}\left(\Phi_{\rm iso}+\phi\right),$ (11)

where $r_{i}=ax_{i}$ and we split the gravitational potential into a
background potential $\Phi_{\rm iso}$ and a peculiar potential $\phi$. The
subscript “iso” in $\Phi$ stands for the potential sourced only by the usual
isotropic background,

$\displaystyle\Phi_{\rm iso}\equiv\frac{2}{3}\pi G\bar{\rho}_{\rm
m}r^{2}-\frac{\Lambda}{6}r^{2},$ (12)

where $\bar{\rho}_{\rm m}$ is the mean density of matter. The peculiar
potential $\phi$ includes the potential sourced by local inhomogeneities
$\phi_{s}$ as well as the effective tidal field sourced by the long-wave
gravitational wave, $\phi=\phi_{s}+\frac{1}{2}\tau_{ij}x^{i}x^{j}$. In the
next sections, the effective tidal potential is absorbed into the back ground
potential, $\Phi=\Phi_{\rm iso}+\frac{1}{2}\tau_{ij}x^{i}x^{j}$ (see Eq. (33)
below). The potential sourced by local inhomogeneities, $\phi_{s}$, satisfies
the Poisson equation

$\displaystyle\nabla_{r}^{2}\phi_{s}=4\pi G\bar{\rho}_{\rm
m}\delta=\frac{3}{2}\frac{\Omega_{\rm m}(\eta)\mathcal{H}^{2}}{a^{2}}\delta,$
(13)

where $\delta=\rho/\bar{\rho}_{\mathrm{m}}-1$ denotes the overdensity field.
We can also decompose the left-hand side of Eq. (11) into the background and
peculiar parts, resulting in

$\displaystyle\mathcal{H}^{\prime}x_{i}$
$\displaystyle=-\frac{\partial}{\partial x_{i}}\Phi_{\rm iso},$ (14)
$\displaystyle x_{i}^{\prime\prime}+\mathcal{H}x_{i}^{\prime}$
$\displaystyle=-\frac{\partial}{\partial
x_{i}}\left(\phi_{s}+\frac{1}{2}\tau_{\ell m}x^{\ell}x^{m}\right).$ (15)

One can verify that Eq. (14) with Eq. (12) gives the usual Friedmann equation.
Eq. (15) can be regarded as the evolution equation for the displacement field
$\Psi_{i}$, which relates Lagrangian position $q_{i}$ to Eulerian position
$x_{i}$ via $x_{i}=q_{i}+\Psi_{i}$. Before shell-crossing the Jacobian
determinant of this mapping gives the overdensity as

$\displaystyle\delta=\left|\frac{\partial x_{i}}{\partial
q_{j}}\right|^{-1}-1=\left|\delta^{\textsc{k}}_{ij}+\frac{\partial\Psi_{i}}{\partial
q_{j}}\right|^{-1}-1.$ (16)

At linear order this reduces to $\delta^{(1)}=-\partial\Psi^{(1)}_{i}/\partial
q_{i}$. We split the linear displacement into the one sourced by pure scalar
contributions and that sourced by $\tau_{ij}$, each of which satisfies

$\displaystyle\Psi^{(1)\prime\prime}_{s,i}+\mathcal{H}\Psi^{(1)\prime}_{s,i}$
$\displaystyle=-\frac{3}{2}\Omega_{m}(\eta)\mathcal{H}^{2}\frac{\partial_{q}^{i}}{\partial^{2}_{q}}\delta^{(1)},$
(17)
$\displaystyle\Psi^{(1)\prime\prime}_{t,i}+\mathcal{H}\Psi^{(1)\prime}_{t,i}$
$\displaystyle=-\frac{1}{2}\partial_{q}^{i}\left[\tau_{kl}q^{k}q^{l}\right],$
(18)

where $\Psi^{(1)}_{s,i}$ and $\Psi^{(1)}_{t,i}$ denote the linear displacement
caused by the scalar and tensor perturbations, respectively, and we have used
Eq. (13) and the fact that $x_{i}$ can be replaced by $q_{i}$ at linear order.
The non-decaying solutions for $\Psi^{(1)}_{s,i}$ and $\Psi^{(1)}_{t,i}$ are

$\displaystyle\Psi^{(1)}_{s,i}(\eta)$
$\displaystyle=-\frac{D(\eta)}{D(\eta_{0})}\frac{\partial_{q}^{i}}{\partial_{q}^{2}}\delta^{(1)}(\eta_{0}),$
(19) $\displaystyle\Psi^{(1)}_{t,i}(\eta;k_{L})$
$\displaystyle=\frac{1}{2}\left[\mathcal{T}(\eta;k_{L})-1\right]h^{\rm
ini}_{ij}(k_{L})q^{j}\equiv-\beta(\eta;k_{L})h_{ij}^{\rm ini}(k_{L})q^{j},$
(20)

where we have used the boundary condition $\lim_{\eta\to
0}\mathcal{T}(\eta;k_{L})=1$ and have introduced the quantity
$\beta(\eta;k_{L})$ to characterize the linear displacement due to the tensor
perturbations, and $D(\eta)$ represents the linear growth function, which
follows

$\displaystyle
D^{\prime\prime}(\eta)+\mathcal{H}D^{\prime}(\eta)-\frac{3}{2}\Omega_{\rm
m}(\eta)\mathcal{H}^{2}D(\eta)=0.$ (21)

From this result one can confirm that GWs do not induce the linear density:
$\delta^{(1)}_{t}=-\partial\Psi^{(1)}_{t,i}/\partial q_{i}\propto h_{i}^{\
i}=0$, as expected.

At second order, however, GWs do affect the density field as well as the
displacement field. Expanding Eq. (16) up to second order leads to

$\displaystyle\delta^{(1)}({\bf x}({\bf q}))+\delta^{(2)}_{st}({\bf x}({\bf
q}))=$ $\displaystyle\delta^{(1)}({\bf
q})-\left.\frac{\partial\Psi^{(2)}_{st,i}}{\partial q_{i}}\right|_{{\bf
q}}+\left.\frac{\partial\Psi^{(1)}_{s,i}}{\partial
q_{j}}\frac{\partial\Psi^{(1)}_{t,j}}{\partial q_{i}}\right|_{{\bf q}}$
$\displaystyle=$ $\displaystyle\delta^{(1)}({\bf x})-\Psi_{t,i}^{(1)}({\bf
x})\partial_{x}^{i}\delta^{(1)}({\bf
x})-\left.\frac{\partial\Psi^{(2)}_{st,i}}{\partial q_{i}}\right|_{{\bf
q}={\bf x}}+\left.\frac{\partial\Psi^{(1)}_{s,i}}{\partial
q_{j}}\frac{\partial\Psi^{(1)}_{t,j}}{\partial q_{i}}\right|_{{\bf q}={\bf
x}},$ (22)

where we focus on the second-order density arising from the coupling between
the scalar and tensor modes, and $\delta_{st,i}^{(2)}$ and $\Psi_{st,i}^{(2)}$
respectively denote the second-order density and displacement induced by the
coupling. Namely, we omit the second-order contributions from the auto-
coupling between scalars (irrelevant) or tensors (subdominant). In deriving
the first equality above, we have used the fact that $\Psi^{(1)}_{t,i}$ is
divergence free. Taking the divergence of Eq. (15) with respect to ${\bf q}$
we obtain

$\displaystyle\psi_{st}^{(2)\prime\prime}+{\cal
H}\psi_{st}^{(2)\prime}=-\nabla_{x}^{2}\phi_{s}-\frac{\partial\Psi^{(1)}_{t,j}}{\partial
q_{i}}\frac{\partial^{2}\phi_{s}}{\partial x_{i}\partial
x_{j}}-\frac{\partial\Psi^{(1)}_{s,j}}{\partial q_{i}}\tau_{ij},$ (23)

where we have defined $\psi^{(2)}_{st}\equiv\partial\Psi_{st,i}^{(2)}/\partial
q_{i}$ and used
$\partial_{q}^{i}=\partial_{x}^{i}+\partial_{q}^{i}\Psi_{j}\partial_{q}^{j}\simeq\partial_{x}^{i}+\partial_{q}^{i}\Psi_{j}\partial_{x}^{j}$.
Rewriting $\phi_{s}$ in terms of $\delta$ by using Eq. (13) and using Eq. (22)
yield

$\displaystyle\psi_{st}^{(2)\prime\prime}+\mathcal{H}\psi_{st}^{(2)\prime}=\frac{3}{2}\Omega_{\rm
m}\mathcal{H}^{2}\left[-\frac{\partial\Psi^{(2)}_{st,i}}{\partial
q_{i}}+\frac{\partial\Psi^{(1)}_{s,i}}{\partial
q_{j}}\frac{\partial\Psi^{(1)}_{t,j}}{\partial
q_{i}}\right]-\frac{3}{2}\Omega_{\rm
m}\mathcal{H}^{2}\frac{\partial\Psi^{(1)}_{t,j}}{\partial
q_{i}}\frac{\partial\Psi^{(1)}_{s,i}}{\partial
q_{j}}-\frac{\partial\Psi^{(1)}_{s,j}}{\partial q_{i}}\tau_{ij}.$ (24)

Finally the equation for $\psi_{st}^{(2)}$ is found to be

$\displaystyle\psi_{st}^{(2)\prime\prime}+\mathcal{H}\psi_{st}^{(2)\prime}-\frac{3}{2}\Omega_{\rm
m}\mathcal{H}^{2}\psi_{st}^{(2)}=\left(\frac{\partial^{i}_{q}\partial^{j}_{q}}{\partial_{q}^{2}}\delta^{(1)}(\eta)\right)\tau_{ij}.$
(25)

We can write the solution as

$\displaystyle\psi_{st}^{(2)}({\bf
q},\eta;k_{L})=D_{st}^{(2)}(\eta;k_{L})\left(\frac{\partial^{i}_{q}\partial^{j}_{q}}{\partial_{q}^{2}}\delta^{(1)}({\bf
q},\eta_{0})\right)h_{ij}^{\rm ini}(k_{L}),$ (26)

where the time-dependent part of $\psi_{st}^{(2)}({\bf q},\eta;k_{L})$, which
we write $D_{st}^{(2)}(\eta;k_{L})$, satisfies

$\displaystyle
D_{st}^{(2)\prime\prime}(\eta;k_{L})+\mathcal{H}D_{st}^{(2)\prime}(\eta;k_{L})-\frac{3}{2}\Omega_{\rm
m}(\eta)\mathcal{H}^{2}D_{st}^{(2)}(\eta;k_{L})=\frac{D(\eta)}{D(\eta_{0})}\cdot\left[-\frac{1}{2a}\left(a\mathcal{T}^{\prime}(\eta;k_{L})\right)^{\prime}\right].$
(27)

We can solve this equation numerically given the cosmological parameters but
here we derive the analytic solution assuming the matter domination for the
later convenience. In the matter-dominated era where we have $\Omega_{\rm
m}=1$, $\mathcal{H}=2/\eta$ and $D(\eta)\propto a\propto\eta^{2}$, we can
write the solution using Green’s function:

$\displaystyle D_{st}^{(2)}(\eta;k_{L})=$
$\displaystyle\frac{1}{D(\eta_{0})}\int_{0}^{\eta}\differential\tilde{\eta}~{}\frac{1}{5}\left[\frac{\eta^{2}}{\tilde{\eta}}-\frac{\tilde{\eta}^{4}}{\eta^{3}}\right]D(\tilde{\eta})\left[-\frac{1}{2a}\left[a\mathcal{T}^{\prime}(\tilde{\eta};k_{L})\right]^{\prime}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{5}\frac{D(\eta)}{D(\eta_{0})}\left[\beta(\eta)+4\int_{0}^{\eta}\differential\tilde{\eta}\left(\frac{\tilde{\eta}}{\eta}\right)^{5}\beta^{\prime}(\tilde{\eta};k_{L})\right].$
(28)

Finally we find the second order density induced by the interaction between
the scalar and tensor perturbation as

$\displaystyle\delta^{(2)}_{st}({\bf x},\eta;k_{L})=$ $\displaystyle
h_{ij}^{\rm ini}\left[\beta(\eta;k_{L})x^{j}\partial_{x}^{i}\delta^{(1)}({\bf
x},\eta)-D_{st}^{(2)}(\eta;k_{L})\left(\frac{\partial^{i}_{x}\partial^{j}_{x}}{\partial_{x}^{2}}\delta^{(1)}({\bf
x},\eta_{0})\right)+\beta(\eta;k_{L})\left(\frac{\partial^{i}_{x}\partial^{j}_{x}}{\partial_{x}^{2}}\delta^{(1)}({\bf
x},\eta)\right)\right]$ $\displaystyle=$ $\displaystyle h_{ij}^{\rm
ini}\left[\alpha(\eta;k_{L})\frac{\partial^{i}_{x}\partial^{j}_{x}}{\partial_{x}^{2}}+\beta(\eta;k_{L})x^{j}\partial_{x}^{i}\right]\delta^{(1)}({\bf
x},\eta),$ (29)

where $\beta$ is defined in Eq. (20) and $\alpha$ is defined as

$\displaystyle\alpha(\eta;k_{L})\equiv$
$\displaystyle-\frac{D(\eta_{0})}{D(\eta)}D_{st}^{(2)}(\eta;k_{L})+\beta(\eta;k_{L})$
$\displaystyle=$
$\displaystyle\frac{4}{5}\left(\beta(\eta;k_{L})-\int_{0}^{\eta}\differential\tilde{\eta}\left(\frac{\tilde{\eta}}{\eta}\right)^{5}\beta^{\prime}(\tilde{\eta};k_{L})\right).$
(30)

Eq. (29) holds not only during the matter domination but also the entire
history of the universe as long as we solve Eq. (27) numerically, whereas the
second equality of Eq. (30) holds only for the matter domination. The first
term in the last line $(\propto\alpha)$ represents the changes in the short-
mode amplitude, called the growth effect, by the coupling between the scalar
tidal field and long-wavelength GWs. The second term in the last line
$(\propto\beta)$ represents the coordinate shift induced by the long-
wavelength GWs, known as the dilation effect. The power spectrum of the
density field in the presence of the long-wavelength GWs then becomes

$\displaystyle P_{\rm mm}({\bf k_{S}},\eta|h_{ij}(k_{L}))=P_{\rm
lin}(k_{S})\left[1+\hat{k}_{S}^{i}\hat{k}_{S}^{j}h_{ij}^{\rm
ini}(k_{L})\left(2\alpha(\eta;k_{L})-\beta(\eta;k_{L})\frac{\partial\ln{P_{\rm
lin}(k_{S},\eta)}}{\partial\ln{k_{S}}}\right)\right].$ (31)

We show $\mathcal{T}(k_{L};\eta)$, $\alpha(k_{L};\eta)$ and
$\beta(k_{L};\eta)$ as a function of the scale factor for various wavenumbers
of GWs (the left panels) and as a function of the wavemunber of GWs at various
redshifts at the right and left panels of Fig.1, respectively. To plot these
functions, we numerically integrated Eqs. (3) and (27) without assuming the
matter domination. In the limit of $k_{L}\eta\to 0$, it is obvious that there
is no physical effect from GWs, i.e., $\alpha(\eta;k_{L})\to 0$ and
$\beta(\eta;k_{L})\to 0$, since GWs are frozen on the super-horizon scales. On
the other hand, taking the limit $k_{L}\eta\gg 1$ is more interesting because
in this limit $\mathcal{T}(\eta;k_{L})\to 0$ but $\alpha(\eta;k_{L})$ and
$\beta(\eta;k_{L})$ do not vanish. In other words, even long after GWs have
decayed away, their impact on the growth and displacement remains, which is
sometimes called the “fossil” effect Masui and Pen (2010).

##  III Tidal separate universe with gravitational waves

In this section, we describe how to incorporate long-wavelength GWs into the
simulation background with the help of the tidal separate universe simulation
technique developed in Ref. Akitsu  _et al._ (2021a) (see also Refs. Stücker
_et al._ (2021); Masaki  _et al._ (2020)). We focus on differences arising
from GWs and refer the readers to Ref. Akitsu  _et al._ (2021a) about the
details of the implementation. In this section we do not employ Einstein’s
summation convention to avoid confusions.

###  III.1 Anisotropic background

In the tidal separate universe simulation we introduce anisotropic scale
factors by absorbing the long-wavelength tidal perturbations into the
background. In general, reflecting that the tidal perturbations are expressed
by the $3\times 3$ symmetric matrix $\tau_{ij}$ the anisotropic scale factors
are also written by the $3\times 3$ symmetric matrix $a_{ij}$, which relates
the physical coordinate $r_{i}$ to the comoving coordinate $x_{i}$ as
$r_{i}=\sum_{j}a_{ij}x_{j}$. However, we can always rotate the simulation
coordinates to align with the eigenvectors of $\tau_{ij}$, leaving only the
diagonal components non-zero so that
$\tau_{ij}=\tau_{i}\delta^{\textsc{k}}_{ij}$ and
$a_{ij}=a_{i}\delta^{\textsc{k}}_{ij}$. Now we have different scale factors
for each axis and characterize these differences by $\Delta_{i}$ defined via

$\displaystyle a_{i}=a(1+\Delta_{i}),$ (32)

with $a$ being the global scale factor. While the equation of motion remains
the same as in Sec. II.2, now we want to absorb the long-wavelength effective
tidal potential induced by GWs into the background potential so that

$\displaystyle\Phi=\Phi_{\rm
iso}+\frac{1}{2}\sum_{i}\tau_{i}x_{i}^{2}=\frac{2}{3}\pi G\bar{\rho}_{\rm
m}r^{2}-\frac{\Lambda}{6}r^{2}+\frac{1}{2}\sum_{i}\tau_{i}x_{i}^{2}.$ (33)

The background and peculiar equations in an anisotropic background become

$\displaystyle\frac{1}{a^{2}}\left[a_{i}^{\prime\prime}-\mathcal{H}a_{i}^{\prime}\right]x_{i}=-\frac{1}{a_{i}}\frac{\partial}{\partial
x_{i}}\Phi,$ (34)
$\displaystyle\frac{1}{a^{2}}\left[a_{i}x_{i}^{\prime\prime}+2a_{i}^{\prime}x_{i}^{\prime}-\mathcal{H}a_{i}x_{i}\right]=-\frac{1}{a_{i}}\frac{\partial}{\partial
x_{i}}\phi_{s}.$ (35)

Subtracting the isotropic background, which is determined by Eq. (14), from
Eq. (34) and linearizing it in $\Delta_{i}$ we find

$\displaystyle\Delta^{\prime\prime}_{i}(\eta;k_{L})+\mathcal{H}\Delta_{i}^{\prime}(\eta;k_{L})=\frac{1}{2}a^{-1}(\eta)\left[a(\eta)h_{i}^{\prime}(\eta;k_{L})\right]^{\prime}.$
(36)

Note that anisotropic scale factors depend on the wavenumber of GWs as a
consequence of the wavenumber dependence of the transfer function of GWs,
$\mathcal{T}(\eta;k_{L})$, while they do not for the case of scalar large-
scale tidal field where the linear growth of scalar perturbations is
independent of their wavenumbers, as characterized by $D(\eta)$. Note also
that the source term for the anisotrpic scale factors on the right-hand side
is non-zero only when $h_{i}^{\prime}(\eta;k_{L})$ does not vanish; in other
words long-wavelength GWs induces the anisotropic scale factors only after GWs
enters the horizon as expected. In fact, integrating Eq. (36) twice yields

$\displaystyle\Delta_{i}(\eta;k_{L})=\frac{h^{\rm
ini}_{i}(k_{L})}{2}\left[\mathcal{T}(\eta;k_{L})-1\right]=-h^{\rm
ini}_{i}(k_{L})\beta(\eta;k_{L}),$ (37)

which goes to zero when $\eta\to 0$ because $\lim_{\eta\to
0}\mathcal{T}(\eta;k_{L})=1$.

The appearance of the function $\beta(\eta;k_{L})$ defined in Eq. (20) is
expected and this result can be understood in a more intuitive way. What is
specifically done in the separate universe construction is to absorb the
displacement caused by the long-wavelength perturbation into the background
expansion while keeping the physical distance unchanged. In other words, we
introduce the local scale factor $a_{i}$ to satisfy

$\displaystyle a_{i}x_{i}=a(x_{i}+\Psi_{i}^{\rm long}),$ (38)

which implies $\Delta_{i}x_{i}=\Psi_{i}^{\rm long}$. Given the displacement
caused by long-wavelength GWs in Eq. (20), this immediately leads to
$\Delta_{i}(\eta;k_{L})=-h_{i}^{\rm ini}(k_{L})\beta(\eta;k_{L})$.

We note that this matching only works for the non-relativistic matter. In
other words, the effect of long-wavelength GWs can be captured by the
anisotropic expansion only when we focus on non-relativistic particles that do
not care about $g_{0i}$ and $g_{ij}$ components in the metric. For example,
the method presented here is not useful in order to study the impact of long-
wavelentgh GWs on the radiation perturbations. However, this treatment is
sufficient to study the influence on dark matter particles that we are
interested in and consistent with the usual Newtonian $N$-body method.

###  III.2 Initial conditions

The background anisotropy induces a correction to the 2LPT solution in the
isotropic background as discussed in Ref.Akitsu  _et al._ (2021a). The
correction depends on $\Delta_{i}(\eta)$, which is different for the scalar
tidal field and GWs. Here we derive this correction induced by the background
anisotropy governed by long-wavelength GWs.

The equation for the displacement in the anisotropic background can be
obtained by combining Eq. (35) with Eq. (13) as

$\displaystyle\sum_{ij}\left|\frac{\partial{\bf x}}{\partial{\bf
q}}\right|\left[\delta_{ij}+\Psi_{i,j}\right]^{-1}\left[\Psi_{i,j}^{\prime\prime}+(\mathcal{H}+2\Delta_{i}^{\prime})\Psi_{i,j}^{\prime}\right]=\frac{3}{2}\Omega_{\rm
m}\mathcal{H}^{2}\left(\left|\frac{\partial{\bf x}}{\partial{\bf
q}}\right|-1\right),$ (39)

where we adopt the notation $\Psi_{i,j}\equiv\partial_{q}^{j}\Psi_{i}$ in this
subsection. Taking $\Delta_{i}=0$ results in the usual master equation for the
LPT. We introduce the correction $\epsilon^{(1)}_{i}$ of order ${\cal
O}(\delta^{(1)}\Delta_{i})$ with $\delta^{(1)}$ being short-wavelength modes
in the simulations:

$\displaystyle\Psi_{i}=\Psi_{i}^{(1)}+\Psi_{i}^{(2)}+\epsilon^{(1)}_{i},$ (40)

where $\Psi_{i}^{(1)}+\Psi_{i}^{(2)}$ is the usual 2LPT solution in the
isotropic background. In the following we also introduce the potentials such
that $\Psi_{i}^{(1)}=-\psi^{(1)}_{,i}$ and
$\epsilon_{i}^{(1)}=-\epsilon^{(1)}_{,i}$ for the convenience. The equation
for $\epsilon^{(1)}$ can be found as

$\displaystyle\sum_{i}\epsilon_{,ii}^{(1)\prime\prime}+\mathcal{H}\sum_{i}\epsilon_{,ii}^{(1)\prime}-\frac{3}{2}\Omega_{\rm
m}\mathcal{H}^{2}\sum_{i}\epsilon_{,ii}^{(1)}=-2\sum_{i}\Delta_{i}^{\prime}\psi_{,ii}^{(1)\prime}.$
(41)

Going to Fourier space and decomposing $\epsilon^{(1)}({\bf k})$ as
$\epsilon^{(1)}({\bf k})=\sum_{i}\hat{k}_{i}^{2}\varepsilon_{i}^{(1)}({\bf
k})$, this equation can be rewritten as

$\displaystyle\varepsilon_{i}^{(1)\prime\prime}+\mathcal{H}\varepsilon_{i}^{(1)\prime}-\frac{3}{2}\Omega_{\rm
m}\mathcal{H}^{2}\varepsilon_{i}^{(1)}=-2\Delta_{i}^{\prime}\psi^{(1)\prime}.$
(42)

Notice that $\varepsilon_{i}^{(1)}$ is different from $\epsilon_{i}^{(1)}$. We
can derive the matter-dominated solution for $\varepsilon_{i}^{(1)}$ for
initial-condition generation. Since Green’s function for this equation is the
same as Eq. (27) the solution is

$\displaystyle\varepsilon_{i}^{(1)}(\eta)=$
$\displaystyle\int_{0}^{\eta}\differential\tilde{\eta}~{}\frac{1}{5}\left[\frac{\eta^{2}}{\tilde{\eta}}-\frac{\tilde{\eta}^{4}}{\eta^{3}}\right]\cdot\left(-2\Delta_{i}^{\prime}(\tilde{\eta})\psi^{(1)\prime}(\tilde{\eta})\right)$
$\displaystyle=$
$\displaystyle-\frac{2}{5}\frac{\psi^{(1)}(\eta)}{D(\eta)}\frac{\Delta_{i}(\eta)}{\beta(\eta)}\int_{0}^{\eta}\differential\tilde{\eta}~{}\left[\frac{\eta^{2}}{\tilde{\eta}}-\frac{\tilde{\eta}^{4}}{\eta^{3}}\right]\cdot
D^{\prime}(\tilde{\eta})\beta^{\prime}(\tilde{\eta})$ $\displaystyle=$
$\displaystyle-\frac{4}{5}\psi^{(1)}(\eta)h_{i}^{\rm
ini}\left[\beta(\eta)-\int_{0}^{\eta}\differential\tilde{\eta}~{}\left(\frac{\tilde{\eta}}{\eta}\right)^{5}\beta^{\prime}(\tilde{\eta})\right]$
$\displaystyle=$ $\displaystyle\psi^{(1)}(\eta)h_{i}^{\rm ini}\alpha(\eta),$
(43)

where we have used $D(\eta)\propto\eta^{2}$ and $\Delta_{i}(\eta)=h_{i}^{\rm
ini}\beta(\eta)$ and $\alpha(\eta)$ is introduced in Eq. (29). Thus

$\displaystyle\epsilon^{(1)}(\eta)=\psi^{(1)}(\eta)\alpha(\eta)\sum_{i}h_{i}^{\rm
ini}\hat{k}_{i}^{2}.$ (44)

This implies that the linear growth function has a direction-dependent
modulation:

$\displaystyle D(\eta,{\bf
k})=D(\eta)\left[1+\alpha(\eta;k_{L})\sum_{i}h_{i}^{\rm
ini}\hat{k}_{i}^{2}\right],$ (45)

which is consistent with Eq. (29). Using Eq. (44), the correction to the
velocity can be computed as

$\displaystyle\epsilon_{i}^{(1)\prime}=$
$\displaystyle\left(f_{\alpha}+f_{1}\right)\mathcal{H}\epsilon_{i}^{(1)},$
(46)

with $f_{\alpha}\equiv\differential\ln\alpha/\differential\ln a$ and
$f_{1}\equiv\differential\ln D/\differential\ln a$. We implement these
modifications in 2LPTIC Crocce  _et al._ (2006) and generate the initial
conditions at $z_{i}=99$.

###  III.3 Simulations

Figure 2:  Fractional anisotropic scale factor in the $x$-axis
$\Delta_{x}(a;k_{L})=a_{x}(a;k_{L})/a-1$ for various wavenumbers of GWs
$k_{L}$ as a function of the scale factor $a$ in the case of $h^{\rm
ini}_{ij}={\rm diag}(0.1,-0.1,0)$. The vertical dashed line represents the
starting redshift $z_{\rm ini}=99$.

We perform $N$-body simulations in the tidal backgrounds with $1024^{3}$
particles in $500~{}\mathrm{Mpc}/h$ boxes. The details of modifications of the
$N$-body code based on Gadget-2 Springel (2005) in the tidal background is
described in Ref. Akitsu  _et al._ (2021a). One important additional
modification to the code was made in the drift operator, which is discussed in
Appendix C.

After rotating $\tau_{ij}$ to align its eigenvectors with the simulation axis,
the remaining degrees of freedom can be completely characterized by two
parameters, which we can parametrize as $\tau_{\rm
e}\equiv-(\tau_{1}-\tau_{2})/2$ and $\tau_{\rm
p}=-\tau_{3}+(\tau_{1}+\tau_{2})/2$. Taking into account the transverse
condition of GWs, we cannot consider $\tau_{\rm p}$-type tides for the tidal
separate universe with GWs unlike the case of the scalar tidal field in Ref.
Akitsu  _et al._ (2021a). Hence in this paper we only consider $\tau_{\rm
e}$-type tides for the background anisotropy. In other words, we consider GWs
propagating along $z$ direction with $+$ mode polarization:

$\displaystyle h^{\rm ini}_{ij}=\begin{pmatrix}\pm\epsilon&0&0\\\
0&\mp\epsilon&0\\\ 0&0&0\end{pmatrix},$ (47)

where we choose $\epsilon=0.1$. Since we aim to measure the response of large-
scale structure to the long-wavelength GWs, the results should not be
dependent on the choice of the direction of the propagation and the
polarization basis.

In order to investigate effects of GWs over a wide range of wavenumbers, we
run tidal separate universe simulations with various wavenumbers of GWs:

$\displaystyle
k_{L}=\left\\{0.0001,\,0.0002,\,0.0005,\,0.001,\,0.002,\,0.005,\,0.01,\,0.02,\,0.05,\,0.1,\,0.2\right\\}\,[h/{\rm
Mpc}].$ (48)

These different wavenumbers of GWs give rise to the different time evolution
of the local anisotropic scale factors, through which long-wavelength GWs
affect the simulated small-scale structure formation. Some examples of the
time evolution of $\Delta_{x}$ are shown in Fig. 2. Given the box size of
$500~{}{\rm Mpc}/h$, some wavenumbers are larger than the fundamental mode in
the simulation: $k_{\rm F}=0.013~{}h/{\rm Mpc}$. For such larger wavemubers,
we cannot neglect the curvature of the long-wavelength modes and treat GW as a
unifiom tidal field over the whole simulation box and the approximation is
violated. Still, these GWs are longer-mode than the halo formation scale so GW
can be seen as a uniform tide in that local region. Therefore we can study the
impact of GWs on halos, in particular, the response of the halo shape as
discussed in Sec. VI, where we come back to this issue again.

For references, we also run fiducial simulations with the isotropic
background, and tidal separate universe simulations induced by the scalar
$\tau_{\rm e}$-type tidal field. These reference simulations share the same
parameters (including random seeds) as those used in the GW tidal separate
universe simulations. For each type of simulations (fiducial and tidal
separate universe with GWs and scalar tides), we run four realizations,
amounting to 100 simulations in total.

We use the AHF code Knollmann and Knebe (2009) to identify dark matter halos
in the simulations by spherical overdensity (SO) regions 200 times as dense as
the mean matter density. We need to identify SO halos in the _global_
coordinates on the _isotropic_ background, i.e., $a_{i}x_{i}/a$, while AHF by
default uses the local simulation coordinates $x_{i}$ on the anisotropic
background. To this end, we modify the AHF code to rescale the coordinate
correspondingly in distance computations.

##  IV Anisotropic power spectrum response: tensor fossils in nonlinear regime

In this section, as a first example of imprints of GWs on large-scale
structure, we present the anisotropic impact on the matter auto-, matter-halo
cross-, and halo auto-power spectra by long-wavelength GWs measured from our
$N$-body simulations.

###  IV.1 Growth-dilation decomposition

As we derived in Sec. II.2, long-wavelength GWs leave the anisotropic imprint
in the matter power spectrum in a given realization of GWs through the
nonlinear interaction of tidal fields. This tidal response consists of two
different contributions as in Eq. (31); the term proportional to $\alpha$ that
modulates the amplitude of the power spectrum and the term proportional to
$\beta$ that modulates the scales. In terms of the tidal separate universe
introduced in Sec. III.1, the former known as the growth effect describes the
changes in the amplitude of short-mode fluctuations measured in the local
anisotropic background, while the latter known as the dilation effect stems
from the anisotropic expansion of the local background with respect to the
global one. Specifically, the response of the power spectrum to long-
wavelength GWs can be decomposed as

$\displaystyle\left.\frac{\differential{\ln P_{\cal G}}}{\differential{h^{\rm
ini}_{ij}}}\right|_{{\bf k}_{\cal G}}=\left.\frac{\differential{\ln P_{\cal
L}}}{\differential{h^{\rm ini}_{ij}}}\right|_{{\bf k}_{\cal G}}=$
$\displaystyle\left.\frac{\partial\ln P_{\cal L}}{\partial h^{\rm
ini}_{ij}}\right|_{{\bf k}_{\cal L}}+\left.\frac{\partial\ln P_{\cal
L}}{\partial\ln k_{{\cal L},i^{\prime}}}\right|_{h_{ij}^{\rm
ini}}\left.\frac{\differential{\ln k_{{\cal L},i^{\prime}}}}{\partial h^{\rm
ini}_{ij}}\right|_{{\bf k}_{\cal G}}$ $\displaystyle\equiv$
$\displaystyle\,\hat{k}_{i}\hat{k}_{j}\left[R^{\rm GW}_{\rm
growth}(k;k_{L})+R^{\rm GW}_{\rm dilation}(k;k_{L})\right],$ (49)

where in the first line the power spectra defined with respect to the global
isotropic background and the local anisotropic background are denoted by
$P_{\cal G}$ and $P_{\cal L}$, respectively, and correspondingly the
wavenumbers by ${\bf k}_{\cal G}$ and ${\bf k}_{\cal L}$. The physical scale
should be unchanged in these two coordinates, $a_{i}x_{{\cal L},i}=ax_{{\cal
G},i}$, which implies $k_{{\cal L},i}=k_{{\cal G},i}(1+\Delta_{i})$. In the
first equality we have used that the variances must be conserved in the
coordinate transformation: $P_{\cal G}\,\differential^{3}{\bf k}_{\cal
G}=P_{\cal L}\,\differential^{3}{\bf k}_{\cal L}$, together with
$|\differential^{3}{\bf k}_{\cal L}/\differential^{3}{\bf k}_{\cal
G}|=\prod_{i=1}^{3}(1+\Delta_{i})=1$ at leading order. In the second line, We
relabel both $k_{\cal G}$ and $k_{\cal L}$ as $k$, since the responses $R^{\rm
GW}_{\rm growth}$ and $R^{\rm GW}_{\rm dilation}$ are already first order in
$h^{\rm ini}_{ij}$ (or $\Delta_{i}$) and hence here we do not need to
distinguish $k_{\cal G}$ and $k_{\cal L}$. Notice that we distinguish $k_{L}$
from $k_{\cal L}$; the former represents the wavenumber of GWs. Notice also
that we have defined the response with respect to the $h_{\rm ini}$, not
$h(z)={\cal T}(z)h_{\rm ini}$, which allows for simple computations of
observables in terms of the primordial tensor mode amplitude. Namely the
cosmology dependence other than the initial amplitude of GWs factorizes out in
the response function.

The power spectrum of biased tracers in the presence of long-wavelength GWs
then has a generic form,

$\displaystyle P_{XY}({\bf
k}|h_{ij}(k_{L}))=P_{XY}(k)\left[1+\hat{k}^{i}\hat{k}^{j}h_{ij}^{\rm
ini}(k_{L})\left(R^{\rm GW}_{{\rm growth};XY}(k;k_{L})+R^{\rm GW}_{{\rm
dilation};XY}(k;k_{L})\right)\right],$ (50)

where $X$ and $Y$ represent tracers being considered. For the matter auto-
power spectrum ($X=Y={\rm m}$), this corresponds to the nonlinear extension of
Eq. (31). We will also consider the matter-halo cross-power spectrum ($X={\rm
m}$ and $Y={\rm h}$) and the halo auto-power spectrum ($X=Y={\rm h}$) below.
Moreover, using these GW power spectrum responses, we can write down the
$X$-$Y$-GWs (scalar-scalar-tensor) bispectrum in squeezed limit as

$\displaystyle\lim_{k_{L}\to
0}B_{XYh_{(\lambda)}}(k,k^{\prime},k_{L})=\hat{k}^{i}\hat{k}^{j}e^{(\lambda)}_{ij}\left[R^{\rm
GW}_{{\rm growth};XY}(k;k_{L})+R^{\rm GW}_{{\rm
dilation};XY}(k;k_{L})\right]P_{XY}(k)P_{h_{(\lambda)}}(k_{L}),$ (51)

where we have defined the $X$-$Y$-GW bispectrum via $\langle X({\bf k})Y({\bf
k}^{\prime})h_{(\lambda)}({\bf k}_{L})\rangle=(2\pi)^{3}\delta^{(3)}_{\rm
D}({\bf k}+{\bf k}^{\prime}+{\bf
k}_{L})B_{XYh_{(\lambda)}}(k,k^{\prime},k_{L})$ and neglected the primordial
contribution.333 Although here we have neglected the primordial scalar-scalar-
tensor bispectrum, its contribution to the local observables appears only at
the order of ${\mathcal{O}}(k_{L}/k)^{2}$ and thus is subdominant for the
single-field inflation (see Ref. Pajer  _et al._ (2013)).

Making use of the growth-dilation decomposition, we can go a little further
than the perturbative result. Because the dilation effect purely captures the
coordinate transformation, using Eq. (37), we can obtain the non-perturbative
result as

$\displaystyle R^{\rm GW}_{{\rm
dilation};XY}(k;k_{L})=\beta(k_{L})\frac{\partial\ln P_{XY}(k)}{\partial\ln
k}.$ (52)

Eqs. (31) and (52) are different because former perturbative result involves
the linear power spectrum while the latter involves the nonlinear power
spectrum and can be applied to a non-linear regime. This means that we can
compute the dilation piece without running simulations even in the nonlinear
regime given the slope of the nonlinear power spectrum. On the other hand, the
growth piece comes from the dynamical effect where we cannot extend the
perturbative result in the nonlinear regime. We thus need to rely on
simulations in order to calibrate the growth response in the nonlinear regime.
Hence, in this paper we focus on measuring the growth term from the tidal
separate universe simulations with GWs. By comparing the measurement of the
growth term of the matter auto-power spectrum in simulations with the
perturbation theory prediction at quasi-nonlinear scales, we can validate our
methodology using the relation:

$\displaystyle\lim_{k\ll k_{\rm NL}}R^{\rm GW}_{\rm
growth;mm}(k;k_{L})=2\alpha(k_{L}).$ (53)

Note that since we assume GWs are long-wavelength modes compared with scalar
perturbations the perturbative result is valid only for the range of $k_{L}\ll
k\ll k_{\rm NL}$.

###  IV.2 Growth response of the matter auto-power spectrum from simulations

In the presence of GWs of $\tau_{\rm e}$-type configuration (Eq. (47)), the
matter auto-power spectrum takes a form of

$\displaystyle P_{\rm mm}({\bf k}|h_{\rm e}(k_{L}))=P_{\rm
mm}(k)\left[1+\frac{2}{3}R^{\rm GW}_{\rm growth;mm}(k;k_{L})\left({\cal
L}_{2}(\hat{k}_{1})-{\cal L}_{2}(\hat{k}_{2})\right)h^{\rm ini}_{\rm
e}\right],$ (54)

where $P_{\rm mm}(k)$ is the nonlinear matter power spectrum in the isotropic
background, ${\cal L}_{2}(x)$ is Legendre polynomial of order two,
$\hat{k}_{1}$ and $\hat{k}_{2}$ represent the $x$ and $y$ components of
$\hat{k}$ respectively, and $h_{\rm e}^{\rm ini}\equiv-\left(h_{11}^{\rm
ini}-h_{22}^{\rm ini}\right)/2=\pm\epsilon$. We can estimate the growth
response by taking the quadrupoles of the power spectrum along both $x$ and
$y$ axes:

$\displaystyle P^{\ell_{\rm e}=2}_{\rm mm}(k|h_{\rm e}(k_{L}))$
$\displaystyle\equiv P^{\ell_{x}=2}_{\rm mm}(k)-P^{\ell_{y}=2}_{\rm mm}(k)$
(55) $\displaystyle=2P_{\rm mm}(k)R^{\rm GW}_{\rm growth;mm}(k;k_{L})h^{\rm
ini}_{\rm e},$ (56)

where $P^{\ell_{i}=2}_{\rm mm}(k)$ ($i=\\{x,y\\}$) are defined as

$\displaystyle P^{\ell_{x}=2}_{XY}(k)$ $\displaystyle\equiv
5\int\frac{\differential^{2}\hat{{\bf k}}}{4\pi}P_{XY}({\bf k}){\cal
L}_{2}(\hat{k}_{1}),$ (57) $\displaystyle P^{\ell_{y}=2}_{XY}(k)$
$\displaystyle\equiv 5\int\frac{\differential^{2}\hat{{\bf
k}}}{4\pi}P_{XY}({\bf k}){\cal L}_{2}(\hat{k}_{2}).$ (58)

This leads to the estimator for $R_{\rm growth;mm}(k;k_{L})$ as

$\displaystyle R^{\rm GW}_{\rm growth;mm}(k;k_{L})=\frac{P^{\ell_{\rm
e}=2}_{\rm mm}(k|h^{\rm ini}_{\rm e}=+\epsilon)-P^{\ell_{\rm e}=2}_{\rm
mm}(k|h^{\rm ini}_{\rm e}=-\epsilon)}{4\epsilon P_{\rm mm}(k)}.$ (59)

Figure 3:  Growth response of matter auto-power spectrum to GWs, $R^{\rm
GW}_{\rm growth}(k_{L};k)$, as a function of the short-wavenmer $k$ at various
wavenumvers of GWs (left) and various redshifts (right), measured from the
simulations. The top-left and bottom-left panels show the responses at $z=0.5$
$z=2$, respectively. The top-right and bottom-right panels show the responses
at $k_{L}=0.0002~{}h/{\rm Mpc}$ and $k_{L}=0.001~{}h/{\rm Mpc}$, respectively.
The dashed lines with each different color represent the predictions from the
perturbation theory at each wavenumber and redshift.

Figure 4:  Growth response of matter auto-power spectrum to GWs, $R^{\rm
GW}_{\rm growth}(k_{L};k)$, as a function of the wavenumber of GWs $k_{L}$ at
$z=2,1,0.5,$ and $0$ (from the left to the right), measured from the
simulations. The upper and lower rows show the results for different $k$,
$k=0.376~{}h/{\rm Mpc}$ and $k=0.881~{}h/{\rm Mpc}$, respectively. The orange
line depicts the “rescaled growth response” to the scalar tides,
$\frac{7}{4}\alpha(k_{L})R_{\rm growth}^{\rm scalar}(k)$ with $R_{\rm
growth}^{\rm scalar}(k)$ measured from the tidal separate universe simulations
with the scalar tides.

We measure the growth response of the matter auto-power spectrum to GWs,
$R^{\rm GW}_{\rm growth;mm}(k;k_{L})$, from the simulations which share the
same initial random phase to reduce the sample variance. Since $R^{\rm
GW}_{\rm growth;mm}$ depends on the wavenumbers of the shote-mode, $k$, and of
GWs, $k_{L}$, as well as redshift $z$, it is hard to show all the dependences
in one figure. Thus, let us start by showing the measured growth response as a
function of $k$ for various $k_{L}$ at $z=0.5$ and $z=2$ in the upper- and
lower-left panels of Fig. 3, respectively. We then show it for
$k_{L}=0.0002h/{\rm Mpc}$ and $k_{L}=0.001h/{\rm Mpc}$ at various redshifts in
the upper- and lower-right panels, respectively. In each panel, the
corresponding perturbation theory predictions, Eq. (53), are shown in the
dashed lines. First, onc can see the excellent agreement between the measured
and predicted growth responses on large scales for all $k_{L}$, which verifies
that our methodology to incorporate GWs in simulations works correctly.
Second, the measured response deviates from the perturbation prediction on
nonlinear scales, in particular for larger $k_{L}$. Although the way it
deviates depends on $k_{L}$, the overall trend is similar among different
$k_{L}$ values; (i) at the smallest scales ($k\gtrsim 2\ h/{\rm Mpc}$) the
growth response decreases for all $k_{L}$ and redshifts. (ii) At earlier
redshift, the growth response is slightly enhanced compared to the
perturbation theory while at lower redshift it is largely suppressed. Third,
the redshift dependence of the growth response is clearer as shown in the
right panels. The agreement between the perturbation theory and the simulation
becomes worse at lower redshifts. At $z=2$ these two are in good agreement up
to $k\sim 0.2\ h/{\rm Mpc}$, while at $z=0$ the simulation results start to
differ from the perturbation theory around $k\sim 0.04\ h/{\rm Mpc}$.

These tendencies we found for the tensor mode above are actually very similar
to those for the case of the scalar tidal field which had been extensively
studied (see App.B and Refs. Akitsu  _et al._ (2021a); Stücker  _et al._
(2021); Masaki  _et al._ (2020) for details). Given the similarity of the
behavior of the growth response on nonlinear scales between the scalar and
tensor cases, it is natural to ask how similar are these two quantitatively.
Let us finish this subsection by answering this question. For this purpose, we
consider a “rescaled growth response” to the scalar tides,
$\frac{7}{4}\alpha(\eta;k_{L})R^{\rm scalar}_{\rm growth;mm}(k)$, where the
scalar tidal response $R^{\rm scalar}_{\rm growth;mm}(k)$ is introduced in Eq.
(96) in App. B. Since $R^{\rm scalar}_{\rm growth;mm}(k)$ approaches $8/7$ at
the large-scale limit (see Eq. (95)), this scalar tidal response matches the
tensor tidal response in the large scale limit (Eq. 53),

$\displaystyle\lim_{k\to 0}\frac{7}{4}\alpha(\eta;k_{L})R^{\rm scalar}_{\rm
growth;mm}(k)=2\alpha(\eta;k_{L}).$ (60)

In Fig.4 we compare this “rescaled growth response” to the scalar tides with
the measured tensor tidal response on nonlinear scales as a function of
$k_{L}$ for various redshifts. The upper and lower panels show the results for
$k=0.376\ h/{\rm Mpc}$ and $k=0.881\ h/{\rm Mpc}$, respectively. Overall, the
rescaled response to the scalar tides captures the general feature of the
growth response to GWs on nonlinear scales. Such agreement can be seen for
both the two wavenumbers at all the redshifts when $k_{L}\lesssim
10^{-3}~{}h/{\rm Mpc}$. On the other hand, for $k_{L}\gtrsim 10^{-3}~{}h/{\rm
Mpc}$ the response to GWs has greater values than the rescaled response,
except for some $k_{L}$ at $z=0$. This difference gets larger at higher
redshifts and larger $k$. These can be attributed to the different time
evolution of the anisotropic scale factors in the scalar tide and GWs cases.
The anisotropic scalar factor induced by the scalar tidal field grows
monotonically in time, following the linear growth rate $D(\eta)$ at leading
order, $\Delta^{\rm scalar}_{i}(\eta)\propto D(\eta)$ (Eq. (94)), whereas that
induced by GWs has the time dependence described by $\beta(\eta;k_{L})$, which
is generally monotonous in time for $k_{L}\lesssim 10^{-3}~{}h/{\rm Mpc}$
while not so for $k_{L}\gtrsim 10^{-3}~{}h/{\rm Mpc}$ (see Fig. 2).
Specifically, for $k_{L}\gtrsim 10^{-3}~{}h/{\rm Mpc}$ the anisotropic scale
factors from GWs reach their asymptotic value at early redshift, when the
anisotropic scale factors from scalar tides are still tiny, which causes the
tidal response to be stronger in these $k_{L}$.

Figure 5:  Growth response of the halo-matter cross-power spectrum to GWs,
$R^{\rm GW;hm}_{\rm growth}(k_{L};k)$, as a function of $k$ for various
$k_{L}$ (left), various redshifts (center), and various halo masses (right),
measured from the simulations. The blue dotted line in the right panel shows
the growth response of the matter auto-power spectrum.  Figure 6:  Similar to
Fig. 5 but the growth response of the halo auto-power spectrum to GWs, $R^{\rm
GW;hh}_{\rm growth}(k_{L};k)$ as a function of $k$ for various $k_{L}$ (left),
various redshifts (center), and various halo masses (right), measured from the
simulations. The upper and lower sets show the results normalized by the halo
auto-power spectrum with and without a shot noise, respectively. For the
latter, the shot noise contribution is subtracted assuming the Poisson
distribution. In the lower-right panel, we do not show the result for the mass
bin, $10^{14}M_{\odot}/h<M_{\rm vir}<10^{14.5}M_{\odot}/h$ because it is too
noisy.

###  IV.3 Growth responses of the halo-matter cross- and halo auto-power
spectra from simulations

The halo-matter cross- and halo auto-power spectra in a local region are also
affected by GWs. The growth responses of these power spectra to GWs can be
estimated from the simulations in the same way as the matter auto-power
spectrum,

$\displaystyle R^{\rm GW}_{\rm growth;XY}(k;k_{L})=\frac{P^{\ell_{\rm
e}=2}_{\rm XY}(k|h^{\rm ini}_{\rm e}=+\epsilon)-P^{\ell_{\rm e}=2}_{\rm
XY}(k|h^{\rm ini}_{\rm e}=-\epsilon)}{4\epsilon P_{\rm XY}(k)},$ (61)

with $XY\in\\{{\rm hm},{\rm hh}\\}$.

In Fig. 5 we show the growth response of the halo-matter cross-power spectrum
to GWs as a function of the wavenumber of the short modes. The left and middle
panels focus on the fixed halo mass ($10^{12}M_{\odot}/h<M_{\rm
vir}<10^{12.5}M_{\odot}/h$), and compares the growth responses for various
wavenumbers of GWs at $z=0.5$ in the left panel and various redshfits at
$k_{L}=0.001~{}h/{\rm Mpc}$ in the middle panel. As for the
$k_{L}$-dependence, the response tends to be greater for larger $k_{L}$,
similar to the case of the matter auto response. On the other hand, the
redshift dependence is slightly different from the matter auto case. For the
halo-matter cross-power spectrum, the growth response at late redshifts ($z=0$
and $z=0.5$) gets slightly enhanced compared to the linear regime and persists
up to $k\lesssim 1~{}h/{\rm Mpc}$, whereas the matter auto response starts to
decrease around $k\simeq 0.1~{}h/{\rm Mpc}$. This can be clearly seen in the
right panel of Fig. 5, which directly compares the growth response of the
matter auto power-spectrum to that of the halo-matter cross-power spectrum
with various halo masses at $z=0.5$ and $k_{L}=0.001~{}h/{\rm Mpc}$. First, it
turns out that the scale-dependence ($k$-dependence) of the growth response
varies with the halo mass. The enhancement of the growth response around
$k\simeq 1~{}h/{\rm Mpc}$ is greater for less massive halos. This is mainly
due to the normalization. Here we define the response with respect to the
halo-matter cross-power spectrum in the fiducial simulation. However if we
define the response with respect to the matter auto-power spectrum, the
response of more massive halos gets more amplified in the nonlinear regime.
Second, the difference of the responses between the matter auto and halo-
matter cross cases stems from the halo biases. In particular at largest scales
it should be explained by the halo tidal bias induced by GWs, which we
investigate in the next section.

Next, let us focus on the growth response of the halo auto-power spectrum. To
obtain it, for the denominator of Eq. (61) we use measurements of the halo
auto-power spectrum with and without the shot noise contribution. For the
latter, the shot noise contribution is subtracted assuming the Poisson
distribution, $1/\bar{n}_{\rm h}$. The upper and lower rows of Fig. 6 show the
resulting response of the halo auto-power spectrum as a function of $k$ with
and without the shot noise, respectively. As with Fig. 5, the left and middle
columns show the responses with the fixed mass range
$10^{12}M_{\odot}/h<M_{\rm vir}<10^{12.5}M_{\odot}/h$ for various wavenumber
of GWs at $z=0.5$ and for $k_{L}=0.001~{}h/{\rm Mpc}$ at various redshfits,
respectively, while the right columns show the result with various halo masses
for $k_{L}=0.001~{}h/{\rm Mpc}$ at $z=0.5$. Overall, the $k_{L}$, redshift,
and halo mass dependencies of the responses for the halo auto spectra both
with and without the shot noise are similar to that for the halo-matter cross
spectrum. In the upper set of Fig. 6, the responses approach zero on small
scales where the halo auto-power spectra are dominated by shot noises. On the
contrary, the peaky features around $k\sim 2~{}h/{\rm Mpc}$ in the lower set
reflects the non-Poissonian behavior of the halo shot noise due to, e.g., the
exclusion effect Smith  _et al._ (2007); Hamaus  _et al._ (2010); Baldauf  _et
al._ (2013). Therefore, the simple Poissonian shot-noise removal leads to
zero-crossing or negative, unphysical halo auto-power spectrum on small
scales, making the responses peaky. Nonetheless, up to $k\sim 0.5~{}h/{\rm
Mpc}$, where the shot noise contribution is still small, the responses are not
suppressed unlike the matter auto case. These trends seen in the response of
the halo-matter cross- and halo auto-power spectra are the same as those in
the scalar tidal field case (see App. B).

##  V Halo tidal bias induced by GWs

Figure 7:  Halo tidal bias induced by GWs as a function of halo mass for
various wavenumber of GWs at $z=0.5$ (the left panel) and at various redsfhits
for $k_{L}=0.001~{}h/{\rm Mpc}$ (the right panel).

As discussed in the previous section, the long-wavelength GWs can affect the
halo density field. In the perturbative regime, this effect should be
characterized in terms of the halo bias. When only scalar perturbations are
considered, which is the standard setup, the halo density field can be
expanded up to the second order as

$\displaystyle\delta_{\rm
h}=b_{1}\delta+\frac{1}{2}b_{2}\delta^{2}+b_{s^{2}}s^{2},$ (62)

where $s^{2}=s_{ij}s^{ij}$ with
$s_{ij}\equiv(\partial_{i}\partial_{j}/\partial^{2}-\delta^{\textsc{k}}_{ij}/3)\delta$.
Here $b_{s^{2}}$ is called the tidal bias that captures the effect of the
tidal fields on the halo density field. Since long-wavelength GWs are locally
indistinguishable from tidal fields induced by the scalar perturbations, it is
natural to expect that there is a tidal bias induced by GWs as well. Then, at
the linear order of GWs, the halo density field should acquire the following
contribution from GWs,

$\displaystyle\delta_{\rm h}^{(2)}(\eta)\supset b^{\rm
GW}_{s^{2}}(\eta;k_{L})s^{ij}(\eta)h_{ij}^{\rm ini}(k_{L}),$ (63)

where $b_{s^{2}}^{\rm GW}$ is the tidal bias coefficient induced by the
coupling between tidal fields induced by GWs and scalar density perturbations.
Note that we define $b_{s^{2}}^{\rm GW}$ with respect to the GWs at the
_initial_ epoch instead of the _same_ time, unlike the biases defined with
respect to the scalar perturbations. Because the GWs are additional degrees of
freedom to the scalar adiabatic perturbations, we expect the tidal bias to
depend on the wavenumber of GWs. Collecting the second order pieces in the
matter density field and halo density field, we find the tree-level halo-
matter-GWs bispectrum in the squeezed limit as

$\displaystyle\lim_{k_{L}\to 0}B_{{\rm
hm}h_{(\lambda)}}(k,k^{\prime},k_{L})=\hat{k}^{i}\hat{k}^{j}e^{(\lambda)}_{ij}\left[2b_{1}\alpha(\eta;k_{L})+2b_{s^{2}}^{\rm
GW}(\eta;k_{L})\right]P_{\rm lin}(k)P_{h_{(\lambda)}}(k_{L}),$ (64)

where we have omitted the dilation piece. This implies that the local halo-
matter power spectrum in the presence of long-wavelength GWs is

$\displaystyle P_{\rm hm}({\bf
k}|h_{ij}(k_{L}))=\left[b_{1}+\left[2b_{1}\alpha(\eta;k_{L})+2b^{\rm
GW}_{s^{2}}(\eta;k_{L})\right]k^{i}k^{j}h_{ij}^{\rm ini}\right]P_{\rm
lin}(k).$ (65)

Thus, we can estimate $b_{s^{2}}^{\rm GW}$ from the growth response of the
halo-matter power spectrum to GWs, which is presented in the previous section.
Although there are several ways to estimate $b_{s^{2}}^{\rm GW}$ from the
growth response of the halo-matter power spectrum, we use the following way to
reduce the uncertainty of $b_{1}$ and the sample variance. We first define the
local linear bias estimator as a ratio of the halo-matter and matter auto-
power spectra including the quadrupoles measured in the simulations with GWs,

$\displaystyle\hat{b}^{\ell_{x,y}}_{1}(k;\pm\epsilon)\equiv\frac{P^{\ell=0}_{\rm
hm}(k;\pm\epsilon)+P^{\ell_{x,y}=2}_{\rm hm}(k;\pm\epsilon)}{P^{\ell=0}_{\rm
mm}(k;\pm\epsilon)+P^{\ell_{x,y}=2}_{\rm
mm}(k;\pm\epsilon)}=\frac{b_{1}\pm(2b_{1}\alpha(\eta;k_{L})+2b_{s^{2}}^{\rm
GW})\epsilon}{1\pm 2\alpha(\eta;k_{L})\epsilon}\simeq\tilde{b}_{1}(k)\pm
2\tilde{b}_{s^{2}}^{\rm GW}(k)\epsilon,$ (66)

Therefore using this local linear bias, the estimator for $b_{s^{2}}^{\rm GW}$
is now

$\displaystyle\tilde{b}^{\rm
GW}_{s^{2}}(k)=\frac{1}{2}\left[\frac{\hat{b}_{1}^{\ell_{x}}(k;+\epsilon)-\hat{b}_{1}^{\ell_{x}}(k;-\epsilon)}{4\epsilon}+\frac{\hat{b}_{1}^{\ell_{y}}(k;+\epsilon)-\hat{b}_{1}^{\ell_{y}}(k;-\epsilon)}{4\epsilon}\right].$
(67)

We estimate $b_{s^{2}}^{\rm GW}$ utilizing the $\chi^{2}$ statistic, defined
as $\chi^{2}=\sum_{k=k_{\rm min}}^{k_{\rm max}}[b_{s^{2}}^{\rm
GW}-\tilde{b}_{s^{2}}^{\rm GW}(k)]^{2}/\sigma^{2}_{b_{s^{2}}^{\rm GW}}(k)$,
where $\sigma^{2}_{b_{s^{2}}^{\rm GW}}(k)$ is the variance of
$\tilde{b}_{s^{2}}^{\rm GW}(k)$ at each $k$-bin measured from simulations.
Adopting $k_{\rm max}=0.08~{}h/{\rm Mpc}$, we obtain the best-fitting value of
$b_{s^{2}}^{\rm GW}$ and its uncertainty by minimizing $\chi^{2}$. We restrict
the measurement of $b_{s^{2}}^{\rm GW}$ to a range of $k_{L}\leq
0.002~{}h/{\rm Mpc}$ because we rely on the perturbative results, which are
valid only when $k_{L}\ll k$.

Figure 8:  Upper panels: Halo tidal bias induced by GWs as a function of
wavenumber of GWs with halo mass of $10^{13}~{}M_{\odot}/h<M_{\rm
vir}<10^{13.5}~{}M_{\odot}/h$. Lower panels: Same as upper panels but with
halo mass of $10^{14}~{}M_{\odot}/h<M_{\rm vir}<10^{14.5}~{}M_{\odot}/h$. From
the left to right , we show the results at $z=1$, $0.5$ and $0$. The blue
points with error bars are the measurement from the simulations. The orange
line represents the ansatz $b_{s^{2}}^{\rm
GW}(k_{L};M)=\frac{7}{4}b_{s^{2}}^{\rm scalar}(M)\alpha(k_{L})$ with the
shaded region being the error coming from $b_{s^{2}}^{\rm scalar}(M)$.

Fig. 7 shows the tidal bias from GWs determined in this way as a function of
halo mass. In the left panel we plot the result for various wavenumbers of GWs
at $z=0.5$, and in the right panel we plot the result for various redshifts at
$k_{L}=0.001~{}h/{\rm Mpc}$. Regardless of the wavenumber of GWs and redshift,
more massive halos have greater absolute values of $b_{s^{2}}^{\rm GW}$, which
has the same trend as the usual tidal bias induced by the scalar tidal fields.
In addition to the mass dependence, the redshift dependence is also similar to
the scalar tidal bias case; the absolute value of the tidal bias is larger at
higher redshift. On the other hand, in contrast to the case of the scalar
tidal bias, the tidal bias from GWs has a particular wavenumber dependence as
shown in the left panel of Fig. 7. This wavenumber-dependence (or scale-
dependence) inherits from the wavenumber-dependent transfer function of GWs
while the growth function for the scalar density fluctuations is independent
of the wavenumber (Eq. (21)). To look at this scale dependence in more detail,
Fig. 8 displays $b_{s^{2}}^{\rm GW}$ as a function of $k_{L}$ for several
redshifts and halo masses. We compare the measurements with the following
ansatz:

$\displaystyle b_{s^{2}}^{\rm
GW}(\eta;k_{L})=\frac{7}{4}\alpha(\eta;k_{L})b^{\rm scalar}_{s^{2}}(\eta),$
(68)

where $b_{s^{2}}^{\rm scalar}$ is the tidal bias induced by the scalar tides,
introduced in Eq. (62) (see also App. B). This ansatz is motivated by the fact
that (i) the tidal bias term stems from the coupling of tidal fields, which is
given by $s^{ij}(z)s_{ij}(z)$ for the scalar tide case ($b_{s^{2}}^{\rm
scalar}$) while $s^{ij}(z)h^{\rm ini}_{ij}$ for the tensor tide case
($b_{s^{2}}^{\rm GW}$), and (ii) the second-order matter density induced by
tidal fields is $\frac{4}{7}s^{ij}(z)s_{ij}(z)$ for the scalar tides case
while $\alpha(z;k_{L})s^{ij}(z)h^{\rm ini}_{ij}(k_{L})$ for the tensor tides
case. Despite the large error bars, overall the measurements are well
explained by this ansatz. This result suggests that we do not have to
introduce a new free bias parameter for GWs at leading order once
$b_{s^{2}}^{\rm scalar}$ is known.

##  VI Intrinsic alignments induced by GWs

Figure 9:  Linear shape bias induced by GWs as a function of halo mass. The
left and right panels show $b_{K}^{\rm GW}$ for various wavenumber of GWs at
$z=0.5$ and for various redshifts at $k_{L}=0.001~{}h/{\rm Mpc}$,
respectively.

As we discussed so far, the tidal fields contribute to the density field only
at the second order because tidal fields are a tensor while density fields are
a scalar. Conversely, the tidal fields should contribute to tensor quantities
at linear order. One well-known observable of tensor quantities in large-scale
structure is the intrinsic alignments of halo (or galaxy) shapes Catelan  _et
al._ (2001); Hirata and Seljak (2004). The deviation of the halo intrinsic
shape from sphere is characterized by ellipticity at the lowest order. In
other words, the halo shape can be described by the trace-free rank-2 tensor
in three-dimensional space, $\gamma_{ij}$.

The linear alignment model of the intrinsic alignment relates the halo shape
with the tidal field as follows:

$\displaystyle\gamma_{ij}(\eta)$ $\displaystyle=b^{\rm
scalar}_{K}(\eta)s_{ij}(\eta),$ (69)

where $b^{\rm scalar}_{K}$ is the linear alignment coefficient or the linear
shape bias, which captures the sensitivity or the response of halo or galaxy
shape to the tidal field. Because long-wavelength GWs are locally equivalent
to the tidal field, we naturally expect that the halo shapes are also aligned
by GWs, namely

$\displaystyle\gamma_{ij}(\eta;k_{L})=b_{K}^{\rm GW}(\eta;k_{L})h^{\rm
ini}_{ij}(k_{L}),$ (70)

where we have introduced the linear shape bias $b_{K}^{\rm GW}$ for GWs, in
analogy with the tidal bias, $b^{\rm scalar}_{K}$. While $b^{\rm scalar}_{K}$
is defined as a response of the halo shape with respect to the scalar tidal
field at the same time, $b^{\rm scalar}_{K}$ is defined as that with respect
to GWs at the _initial_ time when the GWs are frozen, as in the tidal bias
case.

Figure 10:  Linear shape bias induced by GWs as a function of the wavenumber
of GWs for the halos with $10^{13.5}~{}M_{\odot}/h<M_{\rm
vir}<10^{14}~{}M_{\odot}/h$. From the left to the right panels we show the
results at $z=2,1,0.5,$ and $0$, respectively. The orange line depicts the
ansatz introduced in Eq. (75), $b_{K}^{\rm
GW}(k_{L};M,z)=\frac{7}{4}\alpha(k_{L},z)b_{K}^{\rm scalar}(M,z)$, and is not
a fitting. The orange dashed region corresponds to the 1-$\sigma$ error of
$b_{K}^{\rm scalar}$.

The linear shape bias introduced above can be efficiently measured in the
tidal separate universe simulation Akitsu  _et al._ (2021a); Stücker  _et al._
(2021), as the linear halo bias can be obtained precisely in the isotropic
separate universe simulation Li  _et al._ (2016); Baldauf  _et al._ (2016);
Lazeyras  _et al._ (2016). In this paper, we define the halo shape by its
reduced inertia tensor,

$\displaystyle J_{ij}=\sum_{n=1}^{N_{\rm p}}\frac{x_{n,i}x_{n,j}}{x_{n}^{2}},$
(71)

where $N_{\rm p}$ is the number of particles in the halo and $x_{n,i}$ is the
$i$-th component of the particle location with respect to the halo center.
Note that here $x$ represents the distance measured in the isotropic
background, i.e., the physical coordinates. As $\gamma_{ij}$ is regarded as
the trace-free part of $J_{ij}$, the linear alignment model states that this
halo shape can be written as

$\displaystyle
J_{ij}=J_{0}\left[\frac{1}{3}\delta^{\textsc{k}}_{ij}+\gamma_{ij}\right]=J_{0}\left[\frac{1}{3}\delta^{\textsc{k}}_{ij}+b_{K}K_{ij}\right],$
(72)

where $J_{0}$ is the normalization, for which we use the trace part of
$J_{ij}$: $J_{0}={\rm Tr}[J_{ij}]$, and $K_{ij}$ is the tidal field of either
$s_{ij}$ or $h_{ij}^{\rm ini}$. Therefore, we can measure $b_{K}$ as a
response of the one-point function of the ellipticity to the tidal field.
Specifically, introducing the following quantity

$\displaystyle J_{\rm e}\equiv\frac{J_{11}-J_{22}}{2},$ (73)

$b_{K}$ can be obtained as

$\displaystyle b_{K}(M,z)=\frac{J_{\rm e}(M,z|K_{\rm e}^{\rm
ini}=+\epsilon)-J_{\rm e}(M,z|K_{\rm e}^{\rm ini}=-\epsilon)}{2\epsilon
J_{0}(M,z)},$ (74)

where $K_{\rm e}^{\rm ini}=h_{\rm e}^{\rm ini}=\left(h_{11}^{\rm
ini}-h_{22}^{\rm ini}\right)/2$ for the tensor tides and $K_{\rm e}^{\rm
ini}=\left(s_{11}(z)-s_{22}(z)\right)/2$ for the scalar tides.

Fig. 9 shows $b_{K}^{\rm GW}$ as a function of halo mass for various
wavenumbers of GWs at $z=0.5$ (the left panel) and for $k_{L}=0.001~{}h/{\rm
Mpc}$ at various redshifts (the right panel). First, we find that GWs
influence the halo shape, implying that indeed intrinsic alignments can be
induced by GWs. Compared with the measurement of $b_{s^{2}}^{\rm GW}$, the
measurement of $b_{K}^{\rm GW}$ has much greater signal-to-noise ratio as we
use the one-point function to measure $b_{K}^{\rm GW}$ while $b_{s^{2}}^{\rm
GW}$ is measured from the large-scale limit of the power spectrum responses.
Second, the halo-mass and redshift dependences of $b_{K}^{\rm GW}$ are similar
to those of $b_{K}^{\rm scalar}$; that is, massive halos tend to more strongly
align with GWs and the strength of the alignment at fixed halo mass decreases
as the redshift gets smaller (see Fig. 16 and Ref. Akitsu  _et al._ (2021a)
for details).444 Note, however, that for the tensor case $b_{K}^{\rm GW}(z)$
is defined with respect to the _initial_ amplitude of GWs and it thus
represents the strength of the alignment with respect to $h_{ij}^{\rm ini}$,
while for the scalar tides case, $b_{K}^{\rm scalar}(z)$ represents the
strength of the alignment with respect to $s_{ij}(z)$. Taking this difference
into account, however, does not change the trend of the redshift dependence
though, since $s_{ij}\propto D$ that increases monotonically with time.
Third, the unique feature in the case of GWs is that $b_{K}^{\rm GW}$ is
wavenumber-dependent (or scale-dependent), as in the cases with the power
spectrum response and the halo tidal bias. Fig. 10 plots $b_{K}^{\rm GW}$ as a
function of the wavenumber of GWs for the halos with
$10^{13.5}~{}M_{\odot}/h<M_{\rm vir}<10^{14}~{}M_{\odot}/h$ at $z=2,1,0.5,0$
to directly demonstrate this scale dependence. The scale dependence is clearly
seen because of the high $S/N$.

Figure 11:  Ratio of the linear shape bias induced by GWs to the one induced
by scalar tidal fields as a function of the wavenumber of GWs for various halo
masses. From the left to the right, we show the results at $z=2,1,0.5,0$,
respectively. The orange line depicts the ansatz introduced in Eq. (75),
$b_{K}^{\rm GW}(k_{L};M,z)/b_{K}^{\rm
scalar}(M,z)=\frac{7}{4}\alpha(k_{L},z)$, and is not a fitting.

In Fig. 10 we compare the measurements with the following ansatz for the
intrinsic alignment from GWs Schmidt  _et al._ (2014),555 The original ansatz
considered in Ref. Schmidt  _et al._ (2014) is $b_{K}^{\rm
GW}(\eta;k_{L})=\frac{7}{2}\alpha(\eta;k_{L})b_{K}^{\rm scalar}(\eta)$ (see
their Eq. (74)). However, the prefactor should be $7/4$ given the factor of
two that comes from exchanging long-mode and short-mode in the second-order
coupling
($\delta^{(2)}=\frac{2}{7}s_{ij}s_{ij}\to\delta^{(2)}=\frac{4}{7}s^{\rm
long}_{ij}s^{\rm short}_{ij}$).

$\displaystyle b_{K}^{\rm
GW}(\eta;k_{L})=\frac{7}{4}\alpha(\eta;k_{L})b_{K}^{\rm scalar}(\eta).$ (75)

It is natural to assume that the halo shape is determined by the local tidal
environment around the halo, which is affected by long-wavelength tidal field
via nonlinear mode-coupling. The influence of the long-wavelength tides on the
small-scalar tides is captured by the response function that is discussed in
Sec. IV. In particular, in Fig. 4 we examine the relation between the response
function to the scalar tides and GWs and find that $R_{\rm growth;mm}^{\rm
GW}(k;k_{L})=\frac{7}{4}\alpha(\eta;k_{L})R_{\rm growth;mm}^{\rm scalar}(k)$
is a good approximation even in the nonlinear regime. As these responses serve
as the amplitude of the small-scale tides induced by the large-scale tides, we
expect that $\gamma^{\rm scalar}_{ij}\propto R_{\rm growth;mm}^{\rm
scalar}s_{ij}$ and $\gamma^{\rm GW}_{ij}\propto R_{\rm growth;mm}^{\rm
GW}h^{\rm ini}_{ij}$, leading to the above ansatz. The ansatz of Eq. (75) is
shown by the orange lines in Fig. 10. Remarkably, the measurements are well
explained by this prediction for all the redshifts. Fig. 11 investigates if
this agreement holds for all halo mass by displaying the ratio of the linear
shape biases, $b_{K}^{\rm GW}/b_{K}^{\rm scalar}$, which is equal to
$\frac{7}{4}\alpha$ in the ansatz regardless of halo mass. It turns out that
the trend seen in Fig. 10 remains the same for all halo masses. One important
point that follows from this agreement is that the process to determine the
halo shape is not local in time. For, if the halo shape responds to the tidal
field locally in time, meaning that the halo shape is related to the
instantaneous tides: $\gamma^{\rm scalar}_{ij}(\eta)\propto s_{ij}(\eta)$ and
$\gamma^{\rm GW}_{ij}(\eta)\propto\tau_{ij}(\eta)$, we expect $b_{K}^{\rm
GW}(\eta)=b_{K}^{\rm scalar}(\eta)T(\eta;k_{L})/a(\eta)$ with $T$ defined in
Eq. (10); however this is not the case.

Let us finish this section by considering possible reasons that cause the
difference between the measurements and the ansatz. First, the deviation is
also observed around $k_{L}\sim 10^{-3}~{}h/{\rm Mpc}$ at $z=0$. This might be
inherited from the large difference in the response functions to the scalar
and tensor tides shown in the upper row of Fig. 4, which implies that
relatively large-scale tides are more responsible for halo shapes than halo-
scale tides. Second, the deviation of the measurements from the ansatz gets
larger as $k_{L}$ increases. This trend is also consistent with the responses
of the matter power spectrum in Fig. 4. At the same time, this could be partly
because the approximation we employ in this paper is no longer valid for these
large $k_{L}$. In other words, GWs with larger-$k_{L}$ cannot be seen as a
uniform tidal field even in the halo formation region, given that the
Lagrangian halo radius $R_{M}=\left(4\pi\bar{\rho}_{\rm m}/3M_{\rm
vir}\right)^{1/3}$ is close to the wavelength of GWs. Specifically, ignoring
the curvature of GWs can leads the corrections to $b_{K}^{\rm GW}$, which
scales as ${\cal O}(k_{L}^{2}R_{M}^{2})$. In the case of $M_{\rm
vir}=1.0\times 10^{13}M_{\odot}/h$ and $k_{L}=0.1~{}h/{\rm Mpc}$, for example,
the correction is about $(0.1\cdot 3)^{2}\sim 0.1$, which cannot be
negligible. In order to accurately study the impact of such GWs, a different
technique is necessary.

We also note that the results for high-$k_{L}$ GWs can be verified by
analyzing the standard $N$-body simulation with staring an anisotropic initial
power spectrum because the effect of high-$k_{L}$ GWs is almost encoded in the
initial conditions. As shown in Fig. 2, the anisotropic scale factors induced
by high-$k_{L}$ GWs ($k_{L}\gtrsim 10^{-2}~{}h/{\rm Mpc}$) reach their
asymptote already at $z_{\rm ini}$ and thus their effect in the late time can
be absorbed into the overall time-independent rescaling of scale factors,
implying that it does not have a physical effect on structure formation. In
other words, for high-$k_{L}$ GWs halos are formed in absence of a large-scale
tidal field but with anisotropic small-scale modes, while for low-$k_{L}$ GWs
halos are formed in an evolving tidal field, in analogy to the difference in
$b_{1}$ (the response to the large-scale density field) and $b_{\phi}$ (the
response to the change of $\sigma_{8}$) in the density case Baldauf  _et al._
(2011); Desjacques  _et al._ (2018). This study is beyond our paper and we
leave it for future work.

##  VII Discussion

In this paper, we have quantified the impact of GWs on large-scale structure
by means of the tidal separate universe simulation. To the best of our
knowledge, this is the first study for the effect of GWs on nonlinear
structure formation using $N$-body simulations. We found that GWs indeed
influence nonlinear structure formation in both the clustering statistics and
the intrinsic alignment of halo shapes. Our main finding is that the impact of
GWs on large-scale structure can be well described by combining the impact of
the scalar tides with the perturbation theory. Specifically, the halo tidal
bias and linear shape bias induced by GWs, $b_{s^{2}}^{\rm GW}$ and
$b_{K}^{\rm GW}$ respectively, can be approximated by $b_{X}^{\rm
GW}(k_{L})=\frac{7}{4}\alpha(k_{L})b_{X}^{\rm scalar}$ where
$X=\\{s^{2},K\\}$.

Figure 12:  Shape auto-power spectra from GWs at $z=1$ assuming $r=0.1$. The
blue, orange, and green lines show the monopole of the $E$-mode auto-, the
monopole of the $B$-mode auto-, and the dipole of the $EB$ cross-power
spectra, respectively. For the $EB$ cross power spectrum we also assume a
maximally parity-violating case, i.e., $\chi(k)=1$. The blue dashed line
depicts the shape noise.

Let us discuss possible observables that could be used to probe GWs from LSS.
As GWs affect halo shapes at linear order shown in Sec. VI, the simplest probe
would be the shape correlation, namely the intrinsic alignment from GWs
Schmidt and Jeong (2012a); Schmidt  _et al._ (2014); Biagetti and Orlando
(2020). Although the shape correlation is primarily affected by the scalar
perturbations, we can use the $E/B$-decomposition to distinguish the GWs
contribution from the scalar contribution at linear order. Under the flat-sky
approximation and assuming the line-of-sight direction, $\hat{n}$, is parallel
to the $z$-axis, we can define $E$-mode and $B$-mode via

$\displaystyle E({\bf k},\hat{n})\pm iB({\bf k},\hat{n})\equiv{{}_{\pm
2}\gamma}({\bf k},\hat{n})e^{\mp 2i\phi_{k}},$ (76)

where

$\displaystyle{{}_{\pm 2}\gamma}({\bf k},\hat{n})\equiv
m_{\mp}^{i}(\hat{n})m_{\mp}^{j}(\hat{n})\gamma_{ij}({\bf k}),$ (77)

with ${\bf m}_{\pm}\equiv(1,\mp i,0)/\sqrt{2}$. At linear order, the scalar
perturbations only induce $E$-mode with vanishing $B$-mode Hirata and Seljak
(2004). On the other hand, GWs generate both $E$-mode and $B$-mode as

$\displaystyle E({\bf k}_{L},\hat{n})$ $\displaystyle=b_{K}^{\rm
GW}(k_{L})\frac{1}{8}(1+\mu_{L}^{2})\sum_{\lambda}h_{(\lambda)}({\bf k}_{L}),$
(78) $\displaystyle B({\bf k}_{L},\hat{n})$ $\displaystyle=-b_{K}^{\rm
GW}(k_{L})\frac{i}{2}\mu_{L}\sum_{\lambda}\frac{\lambda}{2}h_{(\lambda)}({\bf
k}_{L}),$ (79)

where we have used Eq. (70) and $\mu_{L}\equiv\hat{k}_{L}\cdot\hat{n}$. Note
that here we project the halo shapes onto two-dimensional plane but do not
project their position; we can combine photometric and spectroscopic surveys
to get the projected shapes and their three-dimensional positions (e.g.,
Okumura and Taruya, 2020; Kurita  _et al._ , 2021; Akitsu  _et al._ , 2021b;
Kurita and Takada, 2022). Defining the power spectra for $E$\- and $B$-modes
via

$\displaystyle\langle X({\bf k})Y^{*}({\bf
k}^{\prime})\rangle\equiv(2\pi)^{3}\delta^{(3)}_{\rm D}({\bf k}-{\bf
k}^{\prime})P_{XY}({\bf k}),$ (80)

we obtain

$\displaystyle P_{EE}(k_{L},\mu_{L};z)$
$\displaystyle=\frac{1}{64}(1+\mu_{L}^{2})^{2}(b_{K}^{\rm
GW}(z))^{2}P_{h}(k_{L})\simeq\frac{49}{1024}(1+\mu_{L}^{2})^{2}\alpha^{2}(k_{L};z)(b_{K}^{\rm
scalar}(z))^{2}P_{h}(k_{L}),$ (81) $\displaystyle P_{BB}(k_{L},\mu_{L};z)$
$\displaystyle=\frac{1}{4}\mu_{L}^{2}(b_{K}^{\rm
GW}(z))^{2}P_{h}(k_{L})\simeq\frac{49}{64}\mu_{L}^{2}\alpha^{2}(k_{L};z)(b_{K}^{\rm
scalar}(z))^{2}P_{h}(k_{L}),$ (82)

and $P_{EB}=0$ for unpolarized GWs, where we have used Eq. (75). Considering
the chiral GWs, there appears a non-vanishing $EB$ correlation as

$\displaystyle P_{EB}(k_{L},\mu;z)$
$\displaystyle=\frac{i}{16}\mu_{L}(1+\mu_{L}^{2})(b_{K}^{\rm
GW}(z))^{2}\chi(k_{L})P_{h}(k_{L})\simeq
i\frac{49}{256}\mu_{L}(1+\mu_{L}^{2})\alpha^{2}(k_{L};z)(b_{K}^{\rm
scalar}(z))^{2}\chi(k_{L})P_{h}(k_{L}),$ (83)

where $\chi$ is defined in Eq. (7). Notice that $P_{h}(k_{L})$ is the
primordial power spectrum of the tensor mode.

In Fig. 12, we plot the lowest order moment of the multipoles of these spectra
at $z=1$, i.e., the monopole of the $EE$ and $BB$ spectra and the dipole of
the $EB$ spectrum, assuming $b_{K}^{\rm scalar}=0.1$, $r=0.1$, and $\chi(k)=1$
as a demonstration. On large scales, the suppression comes from the behaviour
of $\alpha(k_{L})$ while on small scales the shape obeys a power-law with
$P(k_{L})\propto k_{L}^{-3}$ because $\alpha(k_{L})$ does not change much. We
also display the shape noise as the dashed line, assuming
$\sigma^{2}_{\gamma}=0.2$ and $\bar{n}_{g}=5.0\times 10^{-4}~{}(h/{\rm
Mpc})^{3}$. As is evident, for all cases the shape noise contribution is much
greater than the expected signals, meaning that the detection of GWs using
these spectra is very challenging. There should also be the contributions of
the scalar-mode both in $EE$ (at linear order) and $BB$ (at one-loop order),
though they are absent in the $EB$ spectrum Biagetti and Orlando (2020).
666Ref. Biagetti and Orlando (2020) argued that the shape noise is absent in
the $EB$ power spectrum, making it a cleaner probe of the chiral GWs. Although
this is true at the signal level, the covariance of the $EB$ spectrum includes
the $EE$ and $BB$ spectra with the shape noise. Thus even for the $EB$ power
spectrum the detectability is limited by the shape noise.  Still, this sort of
probes allows us to put the upper limit on the amplitude of GWs that are
generated after the recombination or the reionization, which is not
constrained by the CMB. Specifically, we could obtain the 1-$\sigma$ error on
$r$ of order $\sigma(r)\sim 10^{3}$ at $k\sim 10^{-3}~{}h/{\rm Mpc}$, implying
the total energy spectrum of GWs, $\Omega_{\rm GW}(k)$, could be constrained
as $\Omega_{\rm GW}(k)\lesssim 10^{-8}$ at $k\sim 10^{-3}~{}h/{\rm Mpc}$ from
the current galaxy surveys.

The density-density-shape bispectrum also involves the GWs contribution. For
instance, the tree-level bispectrum of $\delta_{\rm h}$-$\delta_{\rm h}$-$B$
in squeezed limit ($k_{L}\to 0$) is found to be

$\displaystyle\lim_{k_{L}\to 0}iB^{\rm grav.}_{{\rm hh}B}(k,k_{L};z)=$
$\displaystyle\frac{1}{2}b_{1}(z)b_{K}^{\rm
GW}(k_{L};z)\left[\left(2b_{1}(z)\alpha(k_{L};z)+2b^{\rm
GW}_{s^{2}}(k_{L};z)+b_{1}(z)\beta(k_{L};z)\frac{\partial}{\partial\ln{k}}\right)P_{\rm
lin}(k;z)\right]$ $\displaystyle\hskip
241.84842pt\times\mu_{L}\hat{k}^{i}\hat{k}^{j}\sum_{\lambda}\frac{\lambda}{2}e^{(\lambda)}_{ij}(\hat{k}_{L})P_{h_{(\lambda)}}(k_{L})$
(84) $\displaystyle\simeq$ $\displaystyle\frac{7}{8}b_{1}(z)b_{K}^{\rm
scalar}(z)\alpha(k_{L};z)\left[\left(2(b_{1}(z)+\frac{7}{4}b^{\rm
scalar}_{s^{2}}(z))\alpha(k_{L};z)+b_{1}(z)\beta(k_{L};z)\frac{\partial}{\partial\ln{k}}\right)P_{\rm
lin}(k;z)\right]$ $\displaystyle\hskip
241.84842pt\times\mu_{L}\hat{k}^{i}\hat{k}^{j}\sum_{\lambda}\frac{\lambda}{2}e^{(\lambda)}_{ij}(\hat{k}_{L})P_{h_{(\lambda)}}(k_{L}),$
(85)

where we neglect the contributions from the projection effect Schmidt and
Jeong (2012b) and the redshift-space distortion Kaiser (1987). Notice again
that here $P_{h_{(\lambda)}}(k_{L})$ is the primordial power spectrum and all
the redshift dependence is encoded in the bias coefficients ($b_{K}^{\rm
GW}(k_{L};z)$ and $b^{\rm GW}_{s^{2}}(k_{L};z)$), $\alpha(k_{L};z)$, and
$\beta(k_{L};z)$. The superscript “grav.” emphasises that this bispectrum is
induced by the gravitational interaction in the late time universe. In other
words, there could be an additional contribution from the primordial universe,
e.g., the scalar-scalar-tensor non-Gaussianity Endlich  _et al._ (2013);
Domènech  _et al._ (2017). In principle we can directly observe the scalar-
scalar-tensor non-Gaussianity by looking at this density-density-shape
bispectrum.

Another promising observable to probe GWs from LSS is the quadrupolar
anisotropic imprint in the local density power spectrum discussed in Sec. IV
and Sec. V. One can construct the optimal quadratic estimator for GWs by using
the anisotropic imprint in the local power spectrum Masui and Pen (2010);
Jeong and Kamionkowski (2012); Dimastrogiovanni  _et al._ (2014); Dai  _et
al._ (2013); Masui  _et al._ (2017). Our result can be used to increase
$k_{\rm max}$ of the estimator in this method, allowing for the improved
detectability. As we demonstrated for the first time how GWs affect the biased
tracer, a more realistic estimator for the biased tracer can be available.
Taking the cross correlation with the CMB would also be valuable to improve
the detectability Dodelson (2010); Alizadeh and Hirata (2012); Chisari  _et
al._ (2014); Philcox and Johnson (2022). We leave the detailed investigations
on these possibilities for future work.

Let us conclude by mentioning one potentially significant effect of GWs on LSS
observables. In this paper, we focused on the effect of GWs on the cold dark
matter (CDM) perturbations and ignored the effect on the photon-baryon fluid
in the early universe. However, given that the biased tracers such as galaxies
and halos trace the CDM-baryon perturbation, this could leave a distinctive
signature in LSS observables as well. We will investigate it in future work.

###### Acknowledgements.

We thank Giovanni Cabass, William Coulton, and Matias Zaldarriaga for
insightful discussions and Fabian Schmidt for valuable comments on the draft.
KA is supported by JSPS Overseas Research Fellowships. TO acknowledges support
from the Ministry of Science and Technology of Taiwan under Grants Nos. MOST
110-2112-M-001-045- and 111-2112-M-001-061- and the Career Development Award,
Academia Sinina (AS-CDA-108-M02) for the period of 2019 to 2023. Numerical
computation was carried out on the Helios and the Typhon cluster at the
Institute for Advanced Study and Popeye-Simons cluster at San Diego
Supercomputer Center. The Flatiron Institute is supported by the Simons
Foundation.

##  Appendix A From global FLRW coordinates to local CFC coordinates

Here we outline the mapping from the global perturbed FLRW metric to the local
CFC coordinates. More detailed discussion can be found in Refs. Pajer  _et
al._ (2013); Schmidt  _et al._ (2014); Dai  _et al._ (2015b). We start from
the perturbed FLRW metric,

$\displaystyle ds^{2}=a^{2}(\eta)\left[-(1+h_{00}(\eta,{\bf
x}))d\eta^{2}+\left(\delta^{\textsc{k}}_{ij}+h_{ij}(\eta,{\bf
x})\right)dx^{i}dx^{j}\right],$ (86)

where we neglect $h_{0i}$ component. The coordinate transformation is given by

$\displaystyle x^{0}$
$\displaystyle=x_{F}^{0}+\frac{1}{2}\int_{0}^{x^{0}_{F}}h_{00}(\tilde{\eta})\differential\tilde{\eta}+v_{i}x_{F}^{i}-\frac{1}{4}h^{\prime}_{ij}x_{F}^{i}x_{F}^{j},$
(87) $\displaystyle x^{i}$
$\displaystyle=x_{F}^{i}+v^{i}(x^{0}_{F}-\eta_{F})-\frac{1}{2}h^{i}_{\
j}x_{F}^{j}-\frac{1}{4}\left[h^{i}_{\ j,k}+h^{i}_{\ k,j}-h_{jk}^{\ \
,i}\right]x_{F}^{j}x_{F}^{k},$ (88)

where $v_{i}$ is the coordinate velocity of the central geodesic and
$h_{\mu\nu}$ is evaluated along the central geodesics. The CFC metric can be
obtained by using the transformation law of the metric,
$g_{\mu\nu}^{F}(\eta_{F},{\bf x}_{F})=g_{\alpha\beta}(\eta,{\bf x})(\partial
x^{\alpha}/\partial x^{\mu}_{F})(\partial x^{\beta}/\partial x^{\nu}_{F})$,
yielding

$\displaystyle g_{00}^{F}$
$\displaystyle=-a_{F}^{2}(\eta_{F})\left[1-\frac{1}{2}\left(h_{00,ij}+\mathcal{H}h^{\prime}_{ij}+h^{\prime\prime}_{ij}\right)x_{F}^{i}x_{F}^{j}\right],$
(89) $\displaystyle a_{F}(x^{0}_{F})$
$\displaystyle=a\left(\eta=\eta_{F}+\frac{1}{2}\int_{0}^{\eta_{F}}h_{00}(\tilde{\eta})\differential\tilde{\eta}\right).$
(90)

In Eq. (89), setting $h_{00}=0$ coincides to the main text (Eqs. (9)-(10)). On
the other hand, the case where $h_{00}=-2\Phi_{L}$ and $h_{ij}=0$ corresponds
to the usual tidal separate universe simulation picture (see App. B).

##  Appendix B The tidal responses to the scalar tides from tidal separete
universe simulations

Figure 13:  Growth response of matter auto-power spectrum to the scalar tides,
$R^{\rm scalar;mm}_{\rm growth}(k)$ for various redshifts, measured from the
simulations. The blue dashed line corresponds to the perturbation theory
prediction, $R^{\rm scalar}_{\rm growth;mm}=8/7$ Akitsu  _et al._ (2017).

In this appendix, we summarize the tidal responses to the scalar tidal field.
These include not only the tidal response of the matter auto-power spectrum
and the linear shape bias, which are already presented in Refs. Akitsu  _et
al._ (2021a); Stücker  _et al._ (2021), but also the tidal responses of halo-
matter and halo auto-power spectra and the halo tidal bias.

Figure 14:  Growth responses involving the halo density field as a function of
$k$. In the left panels the responses for various resfhits are plotted while
in the right panels the responses for various halo masses are plotted. The top
panels show the response of the halo-matter cross-power spectrum, $R^{\rm
scalar}_{\rm hm;growth}(k)$, The middle panels show the response of the halo
auto-power spectrum, $R^{\rm scalar}_{\rm hh;growth}(k)$, normalized by the
halo auto-power spectrum with the shot noise, and the bottom panels show the
response of the halo auto-power spectrum, $R^{\rm scalar}_{\rm hh;growth}(k)$,
normalized by the halo auto-power spectrum without the shot noise.

First, we sketch the construction of the tidal separate universe in the scalar
tides case, focusing on the difference to the GWs case, although we refer the
reader to Ref. Akitsu  _et al._ (2021a) for details. The difference between
the GWs case and scalar-tides case in the induced tidal field in the local
region, $\tau_{ij}$, results in the different initial conditions and
anisotropic scale factors. In the scalar-tides case, the induced tides become

$\displaystyle\tau_{ij}(\eta)=-\frac{1}{2}h_{00,ij}$
$\displaystyle=\Phi_{L,ij}$ (91) $\displaystyle=\frac{1}{2}\Omega_{\rm
m}(\eta)\mathcal{H}^{2}\delta_{L}(\eta)\delta^{\textsc{k}}_{ij}+\frac{3}{2}\Omega_{\rm
m}(\eta)\mathcal{H}^{2}s_{L,ij}(\eta),$ (92)

where we have used Poisson equation and decomposed into the trace (large-scale
overdensity, $\delta_{L}$) and the traceless (large-scale pure tidal field,
$s_{L,ij}$). The construction of the tidal separate universe with the scalar
tides corresponds to replacing Eq. (10) with Eq. (92) and repeating the
analysis in Sec. II.2, Sec. III.1, and Sec. III.2. Focusing on the pure tidal
mode ($s_{L,ij}$), the equation that governs the evolution of the anisotropic
scale factors, $\Delta_{i}$, is now

$\displaystyle\Delta^{\prime\prime}_{i}(\eta)+\mathcal{H}\Delta_{i}^{\prime}(\eta)=-\frac{3}{2}\Omega_{\rm
m}(\eta)\mathcal{H}^{2}s_{L,i}(\eta),$ (93)

whose solution is given by

$\displaystyle\Delta_{i}(\eta)=-s_{L,i}(\eta_{0})\frac{D(\eta)}{D(\eta_{0})}.$
(94)

The perturbative prediction of the growth response is

$\displaystyle\lim_{k\to 0}R_{\rm mm;growth}^{\rm scalar}(k)$
$\displaystyle=\frac{8}{7},$ (95)

where the response with respect to the scalar tides is defined through

$\displaystyle\left.\frac{\differential{\ln P_{\cal
G}}}{\differential{s_{L,ij}}}\right|_{{\bf k}_{\cal
G}}=\left.\frac{\differential{\ln P_{\cal
L}}}{\differential{s_{L,ij}}}\right|_{{\bf k}_{\cal G}}=$
$\displaystyle\left.\frac{\partial\ln P_{\cal L}}{\partial
s_{L,ij}}\right|_{{\bf k}_{\cal L}}+\left.\frac{\partial\ln P_{\cal
L}}{\partial\ln k_{{\cal
L},\ell}}\right|_{s_{L,ij}}\left.\frac{\differential{\ln k_{{\cal
L},\ell}}}{\partial s_{L,ij}}\right|_{{\bf k}_{\cal G}}$ $\displaystyle\equiv$
$\displaystyle\,\hat{k}_{i}\hat{k}_{j}\left[R^{\rm scalar}_{\rm
growth}(k)+R^{\rm scalar}_{\rm dilation}(k)\right],$ (96)

which is analogue to Eq. (49). One important difference of Eq. (96) from Eq.
(49) is that in Eq. (96) we define the response with respect to the scalar
tides at the same epoch rather than the initail epoch. As a result, the
perturbative prediction remains the same for all redshifts. Also note that
these results are independent of the wavenumber of long-modes.

Fig. 13 shows the tidal response of the matter auto-power spectrum to the
scalar tidal field as a function of $k$ for various redshifts, measured from
tidal separate universe simulation with the scalar tides. This should be
contrasted with Fig. 3 in the main text. In Fig. 14 we show the tidal
responses of the halo-matter cross-power spectrum (the top panels), the halo
auto-power spectrum with and without the short noise (the middle and the
bottom panels respectively) for various redshifts (the left panels) and for
various halo masses (the right panels). This is analogue to Figs. 5-6 in the
main text. The differences of the responses for various redshifts and halo
masses on large scales should be explained by the halo biases (the combination
of $b_{1}(z;M)$ and $b^{\rm scalar}_{s^{2}}(z;M)$). The peaky feature appeared
in the halo auto response normalized by the halo auto-power without the shot
noise is due to the non-Poissonian feature of the shot noise term and thus not
physical (see the discussion in the last paragraph in Sec. IV.3).

Figure 15:  Halo tidal bias measured from the halo-matter power spectrum
response. The left: The tidal bias as a function of halo mass at $z=0.5$. The
right: The tidal bias as a function of the linear bias $b_{1}$ from various
redshifts and halo masses. The black dashed line displays the Lagrangian
local-in-matter-density (LLIMD) prediction: $b_{s^{2}}^{\rm
scalar}=-\frac{2}{7}(b_{1}-1)$.  Figure 16:  Linear shape bias induced by the
scalar tides as a function of halo mass for various redshifts.

We measure the halo tidal bias $b^{\rm scalar}_{s^{2}}$ in the same way as Eq.
(66) and Eq. (67) in the main text. The result is presented in Fig. 15 where
the left panel shows $b_{s^{2}}^{\rm scalar}$ at $z=0.5$ for various halo
masses and the right panel shows $b_{s^{2}}^{\rm scalar}$ as a function of
$b_{1}$. We compare the result with the Lagrangian local-in-matter-density
(LLIMD) prediction Desjacques  _et al._ (2018): $b_{s^{2}}^{\rm
scalar}=-\frac{2}{7}(b_{1}-1)$, which is plotted as the black-dashed line. In
general the LLIMD prediction fails to capture the behaivour of $b_{s^{2}}^{\rm
scalar}$, in particular at high-mass end, which is consistent with Refs.
Lazeyras and Schmidt (2018); Abidi and Baldauf (2018).

In Fig. 16 we show the linear shape bias (or the linear alignment coefficient)
induced by the scalar tides, $b_{K}^{\rm scalar}$, (introduced in Eq. (69)) as
a function of halo mass for various redshifts. The estimator for $b_{K}^{\rm
scalar}$ used here is the same as Eq. (74). Fig. 16 is a counterpart of Fig.
9.

##  Appendix C The modification in the drift operator

Figure 17:  Ratio of the approximate drift integral to the exact one as a
function of the scale factor. Left: the ratio for different wavenumber of GWs
with the same starting redshift $z_{\rm ini}=99$. Right: the ratio for
different starting redshifts with the same wavenumber of GWs
$k_{L}=0.02~{}h/{\rm Mpc}$.

In the tidal separate universe simulation, the drift operator changes from the
usual $N$-body simulation as

$\displaystyle
x_{i;n+1}=x_{i;n}+\frac{P_{i}}{m}\int^{a_{n+1}}_{a_{n}}\frac{\differential
a}{a_{i}^{2}\mathcal{H}}\equiv x_{i;n}+\frac{P_{i}}{m}~{}{\cal D}^{\rm
exact}_{i;n},$ (97)

where $P_{i}$ is the conjugate momenta of $x_{i}$, $i=x,y,z$ and $n$
represents the time step. In the original Gadget-2 code Springel (2005),
instead of computing the drift integral in Eq. (97) at each time step, first
it prepares the following table

$\displaystyle{\cal D}_{i}[j]\equiv\int^{a[j]}_{a_{\rm
ini}}\frac{\differential a}{a_{i}^{2}\mathcal{H}},$ (98)

where $1\leq j\leq N_{\rm drift}$ with $N_{\rm drift}$ being the length of the
table. $a[j]$ is the $j$-th scale factor that is sampled equally spaced in
logarithm from $a_{\rm ini}$ to $a=1$, regardless of the actual time step. The
actual drift integral at each time step is then evaluated by linearly
interpolating this drift table as

$\displaystyle{\cal D}^{\rm approx.}_{i;n}=$
$\displaystyle\int^{a_{n+1}}_{a_{\rm ini}}\frac{\differential
a}{a_{i}^{2}\mathcal{H}}-\int^{a_{n}}_{a_{\rm ini}}\frac{\differential
a}{a_{i}^{2}\mathcal{H}}$ $\displaystyle\simeq$ $\displaystyle\ {\cal
D}_{i}[j+1]+\left(a_{n+1}-a[j+1]\right)\left({\cal D}_{i}[j+2]-{\cal
D}_{i}[j+1]\right)$ $\displaystyle-{\cal
D}_{i}[j]-\left(a_{n}-a[j]\right)\left({\cal D}_{i}[j+1]-{\cal
D}_{i}[j]\right)$ (99)

with $j=\lfloor n\rfloor$.

This prescription works well for a monotonic integrand, which is the case for
the usual cosmological simulation. We found, however, that this approximation
for the drift operator breaks down for the tidal separate universe simulation
with GWs where the integrand oscillates. Fig. 17 compares the drift factor
evaluated by the above procedure and the direct calculation of the integral in
Eq. (97). Although the default length of the drift table is 1000 we increased
it to $N_{\rm drift}=200000$. Even with this large table, the approximated
drift factor fails to capture the exact result. Therefore we modified the
drift operator so that at each step the drift integral is directly evaluated
without using the drift table or interpolation. This modification is
particularly important for $k_{L}\gtrsim 0.02~{}h/{\rm Mpc}$ when $z_{\rm
ini}=99$. In fact, without this modification the results suffers from the
artifact.

## References

  * Starobinsky (1979) A. A. Starobinsky, JETP Lett. 30, 682 (1979).
  * Sato (1981) K. Sato, Mon. Not. Roy. Astron. Soc. 195, 467 (1981).
  * Guth (1981) A. H. Guth, Phys. Rev. D 23, 347 (1981).
  * Linde (1982) A. D. Linde, Phys. Lett. B 108, 389 (1982).
  * Maggiore (2000) M. Maggiore, Phys. Rept. 331, 283 (2000), arXiv:gr-qc/9909001 .
  * Komatsu  _et al._ (2011) E. Komatsu _et al._ (WMAP), Astrophys. J. Suppl. 192, 18 (2011), arXiv:1001.4538 [astro-ph.CO] .
  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Ade  _et al._ (2018) P. A. R. Ade _et al._ (BICEP2, Keck Array), Phys. Rev. Lett. 121, 221301 (2018), arXiv:1810.05216 [astro-ph.CO] .
  * Adachi  _et al._ (2022) S. Adachi _et al._ (POLARBEAR), Astrophys. J. 931, 101 (2022), arXiv:2203.02495 [astro-ph.CO] .
  * Abazajian  _et al._ (2016) K. N. Abazajian _et al._ (CMB-S4), (2016), arXiv:1610.02743 [astro-ph.CO] .
  * Zaldarriaga and Seljak (1997) M. Zaldarriaga and U. Seljak, Phys. Rev. D 55, 1830 (1997), arXiv:astro-ph/9609170 .
  * Kamionkowski  _et al._ (1997) M. Kamionkowski, A. Kosowsky, and A. Stebbins, Phys. Rev. D 55, 7368 (1997), arXiv:astro-ph/9611125 .
  * Dodelson (2003) S. Dodelson,  _Modern Cosmology_ (Academic Press, Amsterdam, 2003).
  * Suto and Sasaki (1991) Y. Suto and M. Sasaki, Phys. Rev. Lett. 66, 264 (1991).
  * Makino  _et al._ (1992) N. Makino, M. Sasaki, and Y. Suto, Phys. Rev. D 46, 585 (1992).
  * Bernardeau  _et al._ (2002) F. Bernardeau, S. Colombi, E. Gaztanaga, and R. Scoccimarro, Phys. Rept. 367, 1 (2002), arXiv:astro-ph/0112551 .
  * Baumann  _et al._ (2012) D. Baumann, A. Nicolis, L. Senatore, and M. Zaldarriaga, JCAP 07, 051 (2012), arXiv:1004.2488 [astro-ph.CO] .
  * Springel (2005) V. Springel, Mon. Not. Roy. Astron. Soc. 364, 1105 (2005), arXiv:astro-ph/0505010 .
  * Potter  _et al._ (2016) D. Potter, J. Stadel, and R. Teyssier, (2016), arXiv:1609.08621 [astro-ph.IM] .
  * Springel  _et al._ (2021) V. Springel, R. Pakmor, O. Zier, and M. Reinecke, Mon. Not. Roy. Astron. Soc. 506, 2871 (2021), arXiv:2010.03567 [astro-ph.IM] .
  * Angulo and Hahn (2021) R. E. Angulo and O. Hahn,  (2021), 10.1007/s41115-021-00013-z, arXiv:2112.05165 [astro-ph.CO] .
  * Kaiser and Jaffe (1997) N. Kaiser and A. H. Jaffe, Astrophys. J. 484, 545 (1997), arXiv:astro-ph/9609043 .
  * Dodelson  _et al._ (2003) S. Dodelson, E. Rozo, and A. Stebbins, Phys. Rev. Lett. 91, 021301 (2003), arXiv:astro-ph/0301177 .
  * Yoo  _et al._ (2009) J. Yoo, A. L. Fitzpatrick, and M. Zaldarriaga, Phys. Rev. D 80, 083514 (2009), arXiv:0907.0707 [astro-ph.CO] .
  * Jeong  _et al._ (2012) D. Jeong, F. Schmidt, and C. M. Hirata, Phys. Rev. D 85, 023504 (2012), arXiv:1107.5427 [astro-ph.CO] .
  * Schmidt and Jeong (2012a) F. Schmidt and D. Jeong, Phys. Rev. D 86, 083513 (2012a), arXiv:1205.1514 [astro-ph.CO] .
  * Jeong and Schmidt (2012) D. Jeong and F. Schmidt, Phys. Rev. D 86, 083512 (2012), arXiv:1205.1512 [astro-ph.CO] .
  * Masui and Pen (2010) K. W. Masui and U.-L. Pen, Phys. Rev. Lett. 105, 161302 (2010), arXiv:1006.4181 [astro-ph.CO] .
  * Dai  _et al._ (2013) L. Dai, D. Jeong, and M. Kamionkowski, Phys. Rev. D 88, 043507 (2013), arXiv:1306.3985 [astro-ph.CO] .
  * Schmidt  _et al._ (2014) F. Schmidt, E. Pajer, and M. Zaldarriaga, Phys. Rev. D 89, 083507 (2014), arXiv:1312.5616 [astro-ph.CO] .
  * Adamek  _et al._ (2013) J. Adamek, D. Daverio, R. Durrer, and M. Kunz, Phys. Rev. D 88, 103527 (2013), arXiv:1308.6524 [astro-ph.CO] .
  * Adamek  _et al._ (2016a) J. Adamek, D. Daverio, R. Durrer, and M. Kunz, Nature Phys. 12, 346 (2016a), arXiv:1509.01699 [astro-ph.CO] .
  * Adamek  _et al._ (2016b) J. Adamek, D. Daverio, R. Durrer, and M. Kunz, JCAP 07, 053 (2016b), arXiv:1604.06065 [astro-ph.CO] .
  * Domènech and Sasaki (2021) G. Domènech and M. Sasaki, Phys. Rev. D 103, 063531 (2021), arXiv:2012.14016 [gr-qc] .
  * Sirko (2005) E. Sirko, Astrophys. J. 634, 728 (2005), arXiv:astro-ph/0503106 .
  * Gnedin  _et al._ (2011) N. Y. Gnedin, A. V. Kravtsov, and D. H. Rudd, The Astrophysical Journal Supplement Series 194, 46 (2011).
  * Baldauf  _et al._ (2011) T. Baldauf, U. Seljak, L. Senatore, and M. Zaldarriaga, JCAP 10, 031 (2011), arXiv:1106.5507 [astro-ph.CO] .
  * Li  _et al._ (2014) Y. Li, W. Hu, and M. Takada, Phys. Rev. D 89, 083519 (2014), arXiv:1401.0385 [astro-ph.CO] .
  * Wagner  _et al._ (2015) C. Wagner, F. Schmidt, C.-T. Chiang, and E. Komatsu, Mon. Not. Roy. Astron. Soc. 448, L11 (2015), arXiv:1409.6294 [astro-ph.CO] .
  * Lazeyras  _et al._ (2016) T. Lazeyras, C. Wagner, T. Baldauf, and F. Schmidt, JCAP 02, 018 (2016), arXiv:1511.01096 [astro-ph.CO] .
  * Li  _et al._ (2016) Y. Li, W. Hu, and M. Takada, Phys. Rev. D 93, 063507 (2016), arXiv:1511.01454 [astro-ph.CO] .
  * Baldauf  _et al._ (2016) T. Baldauf, U. Seljak, L. Senatore, and M. Zaldarriaga, JCAP 09, 007 (2016), arXiv:1511.01465 [astro-ph.CO] .
  * Stücker  _et al._ (2021) J. Stücker, A. S. Schmidt, S. D. M. White, F. Schmidt, and O. Hahn, Mon. Not. Roy. Astron. Soc. 503, 1473 (2021), arXiv:2003.06427 [astro-ph.CO] .
  * Masaki  _et al._ (2020) S. Masaki, T. Nishimichi, and M. Takada, Mon. Not. Roy. Astron. Soc. 496, 483 (2020), arXiv:2003.10052 [astro-ph.CO] .
  * Akitsu  _et al._ (2021a) K. Akitsu, Y. Li, and T. Okumura, JCAP 04, 041 (2021a), arXiv:2011.06584 [astro-ph.CO] .
  * Dai  _et al._ (2015a) L. Dai, E. Pajer, and F. Schmidt, JCAP 10, 059 (2015a), arXiv:1504.00351 [astro-ph.CO] .
  * Akitsu  _et al._ (2017) K. Akitsu, M. Takada, and Y. Li, Phys. Rev. D 95, 083522 (2017), arXiv:1611.04723 [astro-ph.CO] .
  * Akitsu and Takada (2018) K. Akitsu and M. Takada, Phys. Rev. D 97, 063527 (2018), arXiv:1711.00012 [astro-ph.CO] .
  * Akitsu  _et al._ (2019) K. Akitsu, N. S. Sugiyama, and M. Shiraishi, Phys. Rev. D 100, 103515 (2019), arXiv:1907.10591 [astro-ph.CO] .
  * Li  _et al._ (2018) Y. Li, M. Schmittfull, and U. Seljak, Journal of Cosmology and Astroparticle Physics 2018, 022 (2018).
  * Taruya and Akitsu (2021) A. Taruya and K. Akitsu, JCAP 11, 061 (2021), arXiv:2106.04789 [astro-ph.CO] .
  * Hu  _et al._ (2016) W. Hu, C.-T. Chiang, Y. Li, and M. LoVerde, Phys. Rev. D 94, 023002 (2016), arXiv:1605.01412 [astro-ph.CO] .
  * Chiang  _et al._ (2016) C.-T. Chiang, Y. Li, W. Hu, and M. LoVerde, Phys. Rev. D 94, 123502 (2016), arXiv:1609.01701 [astro-ph.CO] .
  * Chiang  _et al._ (2018) C.-T. Chiang, W. Hu, Y. Li, and M. Loverde, Phys. Rev. D 97, 123526 (2018), arXiv:1710.01310 [astro-ph.CO] .
  * Ade  _et al._ (2016) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO] .
  * Manasse and Misner (1963) F. K. Manasse and C. W. Misner, J. Math. Phys. 4, 735 (1963).
  * Dai  _et al._ (2015b) L. Dai, E. Pajer, and F. Schmidt, JCAP 11, 043 (2015b), arXiv:1502.02011 [gr-qc] .
  * Crocce  _et al._ (2006) M. Crocce, S. Pueblas, and R. Scoccimarro, Mon. Not. Roy. Astron. Soc. 373, 369 (2006), arXiv:astro-ph/0606505 .
  * Knollmann and Knebe (2009) S. R. Knollmann and A. Knebe, Astrophys. J. Suppl. 182, 608 (2009), arXiv:0904.3662 [astro-ph.CO] .
  * Pajer  _et al._ (2013) E. Pajer, F. Schmidt, and M. Zaldarriaga, Phys. Rev. D 88, 083502 (2013), arXiv:1305.0824 [astro-ph.CO] .
  * Smith  _et al._ (2007) R. E. Smith, R. Scoccimarro, and R. K. Sheth, Phys. Rev. D 75, 063512 (2007), arXiv:astro-ph/0609547 .
  * Hamaus  _et al._ (2010) N. Hamaus, U. Seljak, V. Desjacques, R. E. Smith, and T. Baldauf, Phys. Rev. D 82, 043515 (2010), arXiv:1004.5377 [astro-ph.CO] .
  * Baldauf  _et al._ (2013) T. Baldauf, U. Seljak, R. E. Smith, N. Hamaus, and V. Desjacques, Phys. Rev. D 88, 083507 (2013), arXiv:1305.2917 [astro-ph.CO] .
  * Catelan  _et al._ (2001) P. Catelan, M. Kamionkowski, and R. D. Blandford, Mon. Not. Roy. Astron. Soc. 320, L7 (2001), arXiv:astro-ph/0005470 .
  * Hirata and Seljak (2004) C. M. Hirata and U. Seljak, Phys. Rev. D 70, 063526 (2004), [Erratum: Phys.Rev.D 82, 049901 (2010)], arXiv:astro-ph/0406275 .
  * Desjacques  _et al._ (2018) V. Desjacques, D. Jeong, and F. Schmidt, Phys. Rept. 733, 1 (2018), arXiv:1611.09787 [astro-ph.CO] .
  * Biagetti and Orlando (2020) M. Biagetti and G. Orlando, JCAP 07, 005 (2020), arXiv:2001.05930 [astro-ph.CO] .
  * Okumura and Taruya (2020) T. Okumura and A. Taruya, Mon. Not. R. Astron. Soc. 493, L124 (2020), arXiv:1912.04118 [astro-ph.CO] .
  * Kurita  _et al._ (2021) T. Kurita, M. Takada, T. Nishimichi, R. Takahashi, K. Osato, and Y. Kobayashi, Mon. Not. Roy. Astron. Soc. 501, 833 (2021), arXiv:2004.12579 [astro-ph.CO] .
  * Akitsu  _et al._ (2021b) K. Akitsu, T. Kurita, T. Nishimichi, M. Takada, and S. Tanaka, Phys. Rev. D 103, 083508 (2021b), arXiv:2007.03670 [astro-ph.CO] .
  * Kurita and Takada (2022) T. Kurita and M. Takada, Phys. Rev. D 105, 123501 (2022), arXiv:2202.11839 [astro-ph.CO] .
  * Schmidt and Jeong (2012b) F. Schmidt and D. Jeong, Phys. Rev. D 86, 083527 (2012b), arXiv:1204.3625 [astro-ph.CO] .
  * Kaiser (1987) N. Kaiser, Mon. Not. Roy. Astron. Soc. 227, 1 (1987).
  * Endlich  _et al._ (2013) S. Endlich, A. Nicolis, and J. Wang, JCAP 10, 011 (2013), arXiv:1210.0569 [hep-th] .
  * Domènech  _et al._ (2017) G. Domènech, T. Hiramatsu, C. Lin, M. Sasaki, M. Shiraishi, and Y. Wang, JCAP 05, 034 (2017), arXiv:1701.05554 [astro-ph.CO] .
  * Jeong and Kamionkowski (2012) D. Jeong and M. Kamionkowski, Phys. Rev. Lett. 108, 251301 (2012), arXiv:1203.0302 [astro-ph.CO] .
  * Dimastrogiovanni  _et al._ (2014) E. Dimastrogiovanni, M. Fasiello, D. Jeong, and M. Kamionkowski, JCAP 12, 050 (2014), arXiv:1407.8204 [astro-ph.CO] .
  * Masui  _et al._ (2017) K. W. Masui, U.-L. Pen, and N. Turok, Phys. Rev. Lett. 118, 221301 (2017), arXiv:1702.06552 [astro-ph.CO] .
  * Dodelson (2010) S. Dodelson, Phys. Rev. D 82, 023522 (2010), arXiv:1001.5012 [astro-ph.CO] .
  * Alizadeh and Hirata (2012) E. Alizadeh and C. M. Hirata, Phys. Rev. D 85, 123540 (2012), arXiv:1201.5374 [astro-ph.CO] .
  * Chisari  _et al._ (2014) N. E. Chisari, C. Dvorkin, and F. Schmidt, Phys. Rev. D 90, 043527 (2014), arXiv:1406.4871 [astro-ph.CO] .
  * Philcox and Johnson (2022) O. H. E. Philcox and M. C. Johnson, (2022), arXiv:2206.07054 [astro-ph.CO] .
  * Lazeyras and Schmidt (2018) T. Lazeyras and F. Schmidt, JCAP 09, 008 (2018), arXiv:1712.07531 [astro-ph.CO] .
  * Abidi and Baldauf (2018) M. M. Abidi and T. Baldauf, JCAP 07, 029 (2018), arXiv:1802.07622 [astro-ph.CO] .