#  Cosmological simulation in tides: power spectra, halo shape responses, and
shape assembly bias

Kazuyuki Akitsu  Yin Li  and Teppei Okumura

###### Abstract

The well-developed separate universe technique enables accurate calibration of
the response of any observable to an isotropic long-wavelength density
fluctuation. The large-scale environment also hosts tidal modes that perturb
all observables anisotropically. As in the separate universe, both the long
tidal and density modes can be absorbed by an effective anisotropic
background, on which the interaction and evolution of the short modes change
accordingly. We further develop the tidal simulation method, including proper
corrections to the second order Lagrangian perturbation theory (2LPT) to
generate initial conditions of the simulations. We measure the linear tidal
responses of the matter power spectrum, at high redshift from our modified
2LPT, and at low redshift from the tidal simulations. Our results agree
qualitatively with previous works, but exhibit quantitative differences in
both cases. We also measure the linear tidal response of the halo shapes, or
the shape bias, and find its universal relation with the linear halo bias, for
which we provide a fitting formula. Furthermore, analogous to the assembly
bias, we study the secondary dependence of the shape bias, and discover for
the first time the dependence on the halo concentration and axis ratio. Our
results provide useful insights for studies of the intrinsic alignment as a
source of either contamination or information. These effects need to be
correctly taken into account when one uses intrinsic alignments of galaxy
shapes as a precision cosmological tool.

##  1 Introduction

In large-scale structure (LSS) surveys, what we expect to observe is the long-
range correlation of biased tracers (e.g. galaxy number density field)
mediated by long-wavelength perturbations. Long-wavelength perturbations have
an impact on the formation and evolution of the small-scale structure via
nonlinear mode-couplings induced by gravity. Because of the equivalence
principle, the leading-order effects of the long-wavelength gravitational
potential on the local physics arise from the second derivative of the
gravitational potential, which can be decomposed into the large-scale
overdensity and tidal fields. Therefore, it has been a fundamental task of LSS
cosmology to investigate how the large-scale overdensity and tidal field
generate the long-range correlation of biased tracers.

The separate universe simulations provide us with a powerful means to
accurately measure or calibrate the local response of various statistics, such
as the power spectrum and halo mass function, to the large-scale overdensity
[1, 2, 3, 4, 5]. The homogeneity and isotropy of the large-scale overdensity
leads to the simple prescription for the modification, i.e., just changing the
cosmological parameters according to the amplitude of the local overdensity.
On the other hand, large-scale tidal field breaks isotropy and thus the local
background embedded into the large-scale tides is no longer
Friedmann–Lemaître–Robertson–Walker (FLRW) universe but becomes rather
anisotropic [6, 7]. Recently there have appeared some works that incorporate
the large-scale tidal field in $N$-body simulations by introducing anisotropic
scale factors [8, 9, 10, 11]. These simulations enable us to isolate effects
of the large-scale tidal field since we can impose homogeneous tides in the
entire simulation box. In other words, utilizing such simulations we can
robustly measure tidal responses separately from other effects.

In this paper, we implement the large-scale tidal field in $N$-body
simulations together with the appropriate initial condition generator where we
solve the second order Lagrangian perturbation theory (2LPT) in an anisotropic
background. Using our simulations, we measure two kinds of tidal responses:
the matter power spectrum response and the halo shape response. The large-
scale tidal field makes the local clustering pattern anisotropic, which
potentially mimics other anisotropic signatures such as the redshift-space
distortion and Alcock-Paczyński effect [12, 13, 14]. The tidal response of the
matter power spectrum measured from our simulations can give a rough estimate
of such contaminants on small scales. The response also allows us to compute
the covariance of weak lensing power spectra [15]. Since no conclusive results
have yet been reached about the amplitude or scale-dependence of the tidal
response of the matter power spectrum at high redshifts due to possible
numerical artifacts [9, 11], we show the tidal response of the matter power
spectrum from our modified 2LPT.

The density field is not the only one affected by the large-scale tidal field;
shapes of galaxies and halos are naturally affected as well. Indeed, these
shapes are aligned with each other even at large separation due to the large-
scale tidal field, known as “intrinsic alignment” [16, 17, 18]. The intrinsic
alignment is usually considered as a source of the error in measuring weak
lensing signals from galaxy imaging surveys. However, given that its origin is
the underlying large-scale tidal field, it should have their own cosmological
information. For example, the intrinsic alignment potentially probes the
energy budget of the universe [19, 20], the stochastic gravitational waves
background [21, 22], and the angular-dependent primordial non-Gaussianity [23,
24, 25]. In order to extract the cosmological information from measurements of
intrinsic alignments, it is of importance to test theoretical models of the
intrinsic alignment. Our simulation allows us to directly examine the simplest
model of the intrinsic alignment, the so-called linear alignment or tidal
alignment model [16, 17], which predicts shapes of galaxies or halos linearly
aligned with large-scale tides. Using the well-controlled tidal simulations,
we quantify the strength of alignments over a wide redshift and mass range and
explore the secondary dependence of the alignment strength on halo properties
other than halo mass.

This paper is laid out as follows. In Section 2, we present how to absorb the
large-scale tidal and density fields into the simulation background and
modifications of the 2LPT initial condition in anisotropic background. Section
3 summarizes our simulation specifications. In Section 4, we show results on
both the tidal response of the matter power spectrum and the halo shapes. We
give conclusion and discussion in Section 5. Appendices A-C provide the
details of computation and our modifications. We show the convergence test of
our simulations by comparing ours with the usual separate universe simulations
in Appendix D. Results of the different shape definitions used in the main
text are summarized in Appendix E.

##  2 Methodology

A uniform (DC) tidal or density field preserves the translational symmetry,
and can be modeled effectively as a time-dependent coordinate transformation.
By doing this, we separate out the DC modes, and absorb these including their
evolution into an effective background, which we can include in $N$-body
simulations [26]. In this section we present analytical derivations and
numerical implementations of this method. Many of the results have already
been obtained in the recent literature [9]. Here we simplify the derivations,
and present for the first time the modulation of the tidal modes on the second
order Lagrangian perturbation theory (2LPT).

###  2.1 Model uniform tidal field by coordinate transformation

In Newtonian cosmology, the large-scale effect of infinitely long-wavelength
(i.e. DC) density and tidal modes on a dark matter particle can be absorbed by
a coordinate transformation

$r_{i}=a_{ij}x_{j},$ (2.1)

where ${\bm{r}}$ is the physical coordinate of the dark matter particle, and
$a_{ij}$ is a symmetric matrix that absorbs the DC modes so that the large-
scale displacement is isotropic in the ${\bm{x}}$ coordinate. We normalize
$a_{ij}$ to the scale factor of global expansion $a\,\delta_{ij}$ in the
absence of any DC-mode, in which case ${\bm{x}}$ is reduced to the usual
comoving coordinates.

From Eq. (2.1) we can immediately separate the physical velocity $u_{i}$ into
the expansion of a local background and a peculiar component.

$u_{i}\equiv\dot{r}_{i}=H_{ij}r_{j}+v_{i},$ (2.2)

where the overdot denotes a time derivative,
$H_{ij}\equiv\dot{a}_{ik}[a^{-1}]_{kj}$ describes a local anisotropic Hubble
expansion, and $v_{i}\equiv a_{ij}\dot{x}_{j}$ is the peculiar velocity.

The dark matter particles follow the Newtonian equation of motion

$\dot{u}_{i}=-\frac{\partial}{\partial r_{i}}(\Phi+\phi),$ (2.3)

where we split the gravitational potential into an effective background
potential $\Phi$ and a peculiar potential $\phi$. By plugging Eq. (2.2) into
Eq. (2.3), we see that the acceleration also splits into a local background
expansion and a peculiar piece, which are respectively driven by $\Phi$ and
$\phi$,

$\displaystyle\ddot{a}_{ik}[a^{-1}]_{kj}r_{j}$
$\displaystyle=-\partial\Phi/\partial r_{i},$ (2.4)
$\displaystyle\dot{v}_{i}+H_{ij}v_{j}$ $\displaystyle=-\partial\phi/\partial
r_{i}.$ (2.5)

One can regard Eq. (2.4) as a modified Friedmann equation. The large-scale
stress due to the DC modes is absorbed into $\Phi$, leaving $\phi$ sourced
only by local structures,

$\nabla_{\bm{r}}^{2}\phi=4\pi G{\bar{\rho}_{\mathrm{m}}}(1+\Delta_{0})\delta,$
(2.6)

where ${\bar{\rho}_{\mathrm{m}}}$ is the mean density of matter, $\Delta_{0}$
is the large-scale overdensity relative to ${\bar{\rho}_{\mathrm{m}}}$, and
$\delta$ denotes the overdensity with respect to the local background density
${\bar{\rho}_{\mathrm{m}}}(1+\Delta_{0})$. The mass conservation between $a\to
0$ and a later time requires
${\bar{\rho}_{\mathrm{m}}}(1+\Delta_{0})\mathrm{d}^{3}{\bm{r}}={\bar{\rho}_{\mathrm{m}}}a^{3}\mathrm{d}^{3}{\bm{x}}$,
or

$1+\Delta_{0}\equiv\frac{a^{3}}{\det a_{ij}}.$ (2.7)

Without loss of generality, we can simplify the equations and the numerical
implementation by rotating the simulation box to align with the principal axes
of the DC tides, so that $a_{ij}=a_{i}\delta_{ij}$ and
$H_{ij}=H_{i}\delta_{ij}=\dot{a}_{i}\delta_{ij}/a_{i}$ with their off-diagonal
degrees of freedom eliminated. Let us define $\Delta_{i}$ as the relative
difference of $a_{i}$ to $a$:

$1+\Delta_{i}\equiv\frac{a_{i}}{a},\quad i=1,2,3.$ (2.8)

Combining the above two equations implies the mass conservation,

$(1+\Delta_{0})(1+\Delta_{1})(1+\Delta_{2})(1+\Delta_{3})=1.$ (2.9)

We also have

$H_{i}=H+\frac{\dot{\Delta}_{i}}{1+\Delta_{i}}\simeq H+\dot{\Delta}_{i},$
(2.10)

where $H$ is the global expansion rate. The approximation holds at high
redshifts when $\Delta_{i}\ll 1$.

This rotation also diagonalizes the DC tidal field, the traceless part of the
Hessian of the large-scale potential, to $\tau_{ij}=\tau_{i}\delta_{ij}$. Now
we can write down the effective background potential that absorbs the DC
density and tidal modes

$\Phi=\frac{2}{3}\pi
G{\bar{\rho}_{\mathrm{m}}}(1+\Delta_{0})r^{2}-\frac{\Lambda}{6}r^{2}+2\pi
G{\bar{\rho}_{\mathrm{m}}}\tau_{i}r_{i}^{2}.$ (2.11)

The DC density modulation $1+\Delta_{0}$ is determined through Eq. (2.7), and
the last term is sourced by the DC tidal mode. One can easily verify that
substituting the first two terms without $\Delta_{0}$ into Eq. (2.4) gives
rise to the usual Friedmann equation.

To determine the evolution of $\Delta_{i}$ in the presence of $\tau_{i}$, we
plug Eq. (2.11) into Eq. (2.4) and subtract the usual Friedmann equation. Then
it reads

$\ddot{\Delta}_{i}+2H\dot{\Delta}_{i}=-4\pi
G\bar{\rho}_{m}\left(\tau_{i}+\frac{1}{3}\Delta_{0}\right)(1+\Delta_{i}).$
(2.12)

Note that the addition of $\Delta_{0}$ in Eq. (2.11) is valid even at
nonlinear level. This can be verified by setting $\tau_{i}=0$ to reproduce the
evolution of the spherical collapse model (e.g. [2]). For anisotropic
simulations we compute anisotropic scale factors in the following subsections
by solving Eq. (2.12) and Eq. (2.9) numerically, using the matter dominated
initial conditions.

Before proceeding to the next subsection, we can derive some analytic
solutions for better understanding. Linearizing Eq. (2.12) drops the
$1+\Delta_{i}$ factor on the right hand side, and yields an equation with the
same form as that of the linear growth function $D$, if we replace
$-(\tau_{i}+\Delta_{0}/3)$ with $\Delta_{i}$. Because both $\tau_{i}$ and
$\Delta_{0}$ are proportional to $D$, the linear-order solution is simply
given by

$\Delta_{i}^{(1)}=-\tau_{i}-\frac{\Delta_{0}^{(1)}}{3}.$ (2.13)

Unsurprisingly, the linear anisotropic correction to the scale factor is the
sum of the isotropic one $-\Delta_{0}^{(1)}/3$ and the tidal mode. Summing
over $i$ and using the traceless constraint $\sum_{i}\tau_{i}=0$, one can
verify

$\Delta_{0}^{(1)}+\sum_{i}\Delta_{i}^{(1)}=0,$ (2.14)

which is equivalent to Eq. (2.9) in the linear theory limit.

The linear tide $\tau_{i}$ also induces the second order overdensity that is
consistent with the perturbation theory prediction. We derive the second order
solution for $\Delta_{0}$ in App. A, and show here the solution for the matter
dominated era,

$\Delta_{0}^{(2)}=\frac{5}{7}{\Delta_{0}^{(1)}}^{2}+\frac{2}{7}\sum_{i}{\Delta_{i}^{(1)}}^{2}=\frac{17}{21}{\Delta_{0}^{(1)}}^{2}+\frac{2}{7}\sum_{i}\tau_{i}^{2}.$
(2.15)

This is consistent with the result from the second-order standard perturbation
theory (e.g. [27]).

###  2.2 Lagrangian perturbation theory and initial conditions

Due to gravity, the long-wavelength modes are coupled to the short ones and
affect their growth. Here we solve the leading-order anisotropic perturbations
to the first and second order Lagrangian displacement, and use the results to
generate initial conditions for $N$-body simulations. The long modes
$\Delta_{0}$ and $\Delta_{i}$ are assumed to be first order in the following
derivations.

The Lagrangian perturbation theory follows the evolution of the displacement
field ${\bm{\Psi}}({\bm{q}})$, a mapping from a particle’s Lagrangian position
${\bm{q}}$ to its Eulerian position ${\bm{x}}$:

$x_{i}=q_{i}+\Psi_{i}({{\bm{q}}}).$ (2.16)

Before shell crossing, the overdensity is simply related to its Jacobian
determinant by

$\delta=\Bigl{|}\frac{\partial{{\bm{x}}}}{\partial{{\bm{q}}}}\Bigr{|}^{-1}-1=\bigl{|}\delta_{ij}+\Psi_{i,j}\bigr{|}^{-1}-1,$
(2.17)

where $\Psi_{i,j}\equiv\partial\Psi_{i}/\partial q_{j}$.

To derive the master equation of the Lagrangian perturbation theory, we
substitute Eq. (2.16) into Eq. (2.5),

$\ddot{\Psi}_{i}^{(1)}+2\sum_{i}H_{i}\dot{\Psi}_{i}^{(1)}=-\sum_{i}\frac{1}{a_{i}}\frac{\partial\phi}{\partial
r_{i}}.$ (2.18)

Taking the derivative with respect to $x_{i}$ and then summing over $i$, we
obtain the master equation

$\sum_{ij}~{}\Bigl{|}\frac{\partial{{\bm{x}}}}{\partial{{\bm{q}}}}\Bigr{|}\bigl{[}\delta_{ij}+\Psi_{i,j}\bigr{]}^{-1}\Bigl{[}\ddot{\Psi}_{i,j}+2H_{i}\dot{\Psi}_{i,j}\Bigr{]}=4\pi
G{\bar{\rho}_{\mathrm{m}}}(1+\Delta_{0})\Bigl{(}\Bigl{|}\frac{\partial{{\bm{x}}}}{\partial{{\bm{q}}}}\Bigr{|}-1\Bigr{)},$
(2.19)

where we have used Eq. (2.6) and Eq. (2.17), and the chain rule
$\partial/\partial x_{i}=\bigl{[}\partial q_{j}/\partial
x_{i}\bigr{]}\partial/\partial q_{j}$.

Now let us start with the Zel’dovich approximation (ZA; the linear Lagrangian
perturbation theory [28]). Keeping only the leading order displacement terms
in Eq. (2.19) leads to

$\sum_{i}\ddot{\Psi}_{i,i}^{(1)}+2\sum_{i}H_{i}\dot{\Psi}_{i,i}^{(1)}=\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)(1+\Delta_{0})\sum_{i}\Psi_{i,i}^{(1)}.$
(2.20)

In deriving the above equation we have used
$\bigl{|}\partial{{\bm{x}}}/\partial{{\bm{q}}}\bigr{|}^{-1}\simeq
1+\sum_{i}\Psi_{i,i}$ and
${\Omega_{\mathrm{m}}}(a)\equiv{\bar{\rho}_{\mathrm{m}}}/{\rho_{\mathrm{c}}}=8\pi
G{\bar{\rho}_{\mathrm{m}}}/3H^{2}$, with ${\rho_{\mathrm{c}}}$ being the
critical density.

At linear order the vorticity in $\Psi_{i,j}^{(1)}$ decays, so the growing
displacement solution is a potential flow
$\Psi_{i}^{(1)}=-\partial\psi_{W}^{(1)}/\partial
q_{i}\equiv-\psi_{W,i}^{(1)}$, with the potential $\psi_{W}^{(1)}$ sourced by
the overdensity in Lagrangian space:

$\nabla_{{\bm{q}}}^{2}\psi_{W}^{(1)}\equiv\sum_{i}\psi_{W,ii}^{(1)}=\delta^{(1)}.$
(2.21)

In terms of $\psi_{W}^{(1)}$, Eq. (2.20) is

$\sum_{i}\ddot{\psi}^{(1)}_{W,ii}+2\sum_{i}H_{i}\dot{\psi}^{(1)}_{W,ii}-\frac{3}{2}{\Omega_{\mathrm{m}}}(a)(1+\Delta_{0})\sum_{i}\psi^{(1)}_{W,ii}=0.$
(2.22)

The subscript $W$ here denotes local quantities inside a _window_ , within
which the DC density and tidal modes can be nonzero.

This is in contrast with the usual linear growth equation, which describes the
evolution of the short modes in the global background where the long modes
vanish. It can be obtained by setting $H_{i}\to H$ and $\Delta_{0}\to 0$ in
the above equation,

$\sum_{i}\ddot{\psi}^{(1)}_{ii}+2H\sum_{i}\dot{\psi}^{(1)}_{ii}-\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)\sum_{i}\psi^{(1)}_{ii}=0.$
(2.23)

The solution to this equation gives the usual time dependence by the linear
growth function, $\psi^{(1)}\propto D(t)$.

In the presence of the long modes, the displacement potential $\psi^{(1)}_{W}$
receives corrections of order ${\cal O}(\psi^{(1)}\Delta)$. We denote this
correction as $\epsilon^{(1)}$ so that the solution of Eq. (2.22) can be
written as

$\psi_{W}^{(1)}=\psi^{(1)}+\epsilon^{(1)}.$ (2.24)

Then using Eq. (2.10), Eq. (2.22), and Eq. (2.23) one can show that
$\epsilon^{(1)}$ satisfies

$\sum_{i}\ddot{\epsilon}^{(1)}_{,ii}+2H\sum_{i}\dot{\epsilon}^{(1)}_{,ii}-\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)\sum_{i}\epsilon^{(1)}_{,ii}=-2\sum_{i}\dot{\psi}_{,ii}^{(1)}\dot{\Delta}_{i}+\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)\Delta_{0}\sum_{i}\psi_{,ii}^{(1)}.$
(2.25)

We can solve this equation by rewriting it in Fourier space

$\ddot{\epsilon}^{(1)}+2H\dot{\epsilon}^{(1)}-\frac{3}{2}{\Omega_{\mathrm{m}}}(a)\epsilon^{(1)}=-2\dot{\psi}^{(1)}\sum_{i}\hat{p}_{i}^{2}\dot{\Delta}_{i}+\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)\psi^{(1)}\Delta_{0}.$
(2.26)

Here in this subsection (and App. B) we use ${\bm{p}}$ to denote the local
Lagrangian space wavevector, with $\hat{\bm{p}}$ being its direction. The
above equation clearly shows that the effect of the long modes manifests in
the quadrupolarly direction-dependent Hubble drag, whose coefficients depend
on the growth history of the long modes $\Delta_{i}(t)$.

To solve the above equation we can first decompose it as

$\epsilon^{(1)}(a,\hat{{\bm{p}}})=\sum_{i}\hat{p}_{i}^{2}\epsilon^{(1)}_{i}(a),$
(2.27)

with each component $\epsilon^{(1)}_{i}(a)$ satisfying

$\ddot{\epsilon}_{i}^{(1)}+2H\dot{\epsilon}_{i}^{(1)}-\frac{3}{2}{\Omega_{\mathrm{m}}}(a)\epsilon_{i}^{(1)}=-2\dot{\psi}^{(1)}\dot{\Delta}_{i}+\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)\psi^{(1)}\Delta_{0}.$
(2.28)

Note that $\epsilon_{i}^{(1)}$ here is different from $\epsilon^{(1)}_{,ii}$
in Eq. (2.25). For the matter dominated era, assuming that $H=2/3t$,
${\Omega_{\mathrm{m}}}(a)=1$, and the long modes are well sub-horizon
($\Delta_{0},\Delta_{i}\propto D^{(1)}$), we obtain

$\epsilon^{(1)}_{i}=-\frac{4}{7}\Delta_{i}\psi^{(1)}+\frac{3}{7}\Delta_{0}\psi^{(1)}.$
(2.29)

Now we can write down the Fourier-space correction to ZA due to the long modes

$\epsilon^{(1)}=\psi^{(1)}\Bigl{(}\frac{3}{7}\Delta_{0}-\frac{4}{7}\sum_{i}\Delta_{i}\hat{p}_{i}^{2}\Bigr{)}.$
(2.30)

This can be understood as a direction-dependent modulation on the linear
growth function

$D_{W}(t,\hat{\bm{p}})=D(t)\Bigl{(}1+\frac{3}{7}\Delta_{0}-\frac{4}{7}\sum_{i}\Delta_{i}\hat{p}_{i}^{2}\Bigr{)},$
(2.31)

where $D_{W}(t,\hat{\bm{p}})$ is the modified linear growth function for a
Fourier mode along $\hat{\bm{p}}$. Eq. (2.31) can be directly compared to the
results derived with the standard perturbation theory in the Einstein de-
Sitter universe. For an isotropic perturbation, $\Delta_{i}=-\Delta_{0}/3$, so
$D^{(1)}_{W}=D^{(1)}(1+13\Delta_{0}/21)$ and for pure tides, $\Delta_{0}=0$
and $\Delta_{i}=-\tau_{i}$, so
$D^{(1)}_{W}=D^{(1)}(1+4\sum_{i}\tau_{i}\hat{p}_{i}^{2}/7)$.

The above derivations have shown that the tidal effect on ZA is simpler in
Fourier space, and can be captured by a direction-dependent modulation on the
linear growth function. However, this is not the case for 2LPT, for which we
find the correction more straightforward in configuration space. Again we
define the second order displacement potential and its correction by
$\psi_{W}^{(2)}=\psi^{(2)}+\epsilon^{(2)}$. For the matter dominated era, we
find the following solution

$\displaystyle\sum_{i}\epsilon^{(2)}_{,ii}(t,{\bm{q}})=$
$\displaystyle\frac{1}{4}\left[-\frac{16}{9}\sum_{i}\psi^{(2)}_{,ii}(t,{\bm{q}})\Delta_{i}+\frac{8}{9}\sum_{ij}\psi^{(1)}_{,ij}(t,{\bm{q}})\psi^{(1)}_{,ji}(t,{\bm{q}})\Delta_{i}\right]$
$\displaystyle+\frac{1}{6}\left[\sum_{i}\psi_{,ii}^{(2)}-\frac{1}{2}\sum_{i}\left(\psi^{(1)}_{,ii}\right)^{2}-\frac{1}{2}\sum_{ij}\psi^{(1)}_{,ij}\psi^{(1)}_{,ji}\right]\Delta_{0}$
$\displaystyle+\frac{1}{4}\left[-\frac{2}{3}\sum_{i}\psi^{(1)}_{,ii}(t,{\bm{q}})\sum_{j}\epsilon_{,jj}^{(1)}(t,{\bm{q}})+\frac{20}{9}\sum_{ij}\psi^{(1)}_{,ij}(t,{\bm{q}})\epsilon_{,ij}^{(1)}(t,{\bm{q}})\right].$
(2.32)

We present the general formula and its derivation in App. B.

Having derived the leading order perturbations to ZA and 2LPT by the DC modes,
we can implement them numerically to generate initial conditions for our tidal
simulations. We modify the initial condition code 2LPTIC [29] to include the
$\epsilon^{(1)}$ and $\epsilon^{(2)}$ corrections to the displacements, as
well as the corresponding corrections to the initial velocities. Since
$\epsilon^{(1)}\propto D^{2}$ and $\epsilon^{(2)}\propto D^{3}$ during matter
domination, the corrections to ZA and 2LPT velocities are respectively given
by

$\displaystyle a\dot{\epsilon}^{(1)}_{,i}$
$\displaystyle=2Hf_{1}\epsilon^{(1)}_{,i},$ $\displaystyle
a\dot{\epsilon}^{(2)}_{,i}$ $\displaystyle=3Hf_{1}\epsilon^{(2)}_{,i}.$ (2.33)

with $f_{1}\equiv\mathrm{d}\ln D/\mathrm{d}\ln
a\simeq[{\Omega_{\mathrm{m}}}(a)]^{5/9}$. Since we generate initial conditions
at $a_{\mathrm{i}}$ deep in the matter dominated era, the above approximations
should be accurate.

There is one other thing to note about the velocity modifications. Here we
derive the corrections to the peculiar velocity $a\dot{\bm{x}}$. However, the
series of Gadget codes uses the canonical momentum ${\bm{P}}\equiv
a^{2}m\dot{\bm{x}}$ as the internal velocity variable, and another velocity
variable ${\bm{u}}\equiv\sqrt{a}\dot{{\bm{x}}}$ in their data format when
saving and loading data. Therefore, given that we are dealing with anisotropic
scale factors, we need to be careful when converting among peculiar velocity,
momentum, and velocity variable on disks. We multiplied
$\sqrt{a}\dot{\Psi}_{i}$ by $\sqrt{1+\Delta_{i}}$ when generating and saving
initial conditions, and then modified the conversion factor $a\sqrt{a}$ from
${\bm{u}}$ to ${\bm{P}}$ as $a_{i}\sqrt{a_{i}}$, which eventually results in
the modified canonical momentum $P_{i}=a_{i}^{2}\dot{x}_{i}$.

###  2.3 Particle-mesh and tree forces

Because the new effective background evolves anisotropically, the
gravitational force and the equation of motion, isotropic in physical
coordinates, need to be modified and expressed in the local comoving
coordinates. We describe the modification of the force law here, and explain
the time integration in the next subsection.

We focus on the TreePM method which computes gravitational force efficiently
by splitting it into the long-range and short-range contributions, computed by
the particle mesh (PM) method and the tree algorithm [30, 31], respectively.
The PM forces can be solved efficiently in Fourier space, and the tree forces
of nearby particles are summed with the help of a tree data structure. In the
absence of the DC modes,

$\displaystyle\phi$ $\displaystyle=\phi^{\mathrm{PM}}+\phi^{\mathrm{T}},$
$\displaystyle\phi^{\mathrm{PM}}({\bm{k}})$ $\displaystyle=-4\pi
G{\bar{\rho}_{\mathrm{m}}}a^{2}\frac{\delta({\bm{k}})}{k^{2}}e^{-k^{2}x_{\mathrm{s}}^{2}},$
$\displaystyle\phi^{\mathrm{T}}({\bm{x}})$
$\displaystyle=-\frac{Gm}{a}\sum_{n}\frac{1}{|{\bm{x}}-{\bm{x}}_{n}|}\operatorname{erfc}\Bigl{(}\frac{|{\bm{x}}-{\bm{x}}_{n}|}{2x_{\mathrm{s}}}\Bigr{)},$
(2.34)

where ${\bm{k}}$ is the wavevector, ${\bm{x}}_{n}$ denotes the position of the
$n$-th particle, and the overdensity field is determined by the spatial
distribution of the particles

${\bar{\rho}_{\mathrm{m}}}\bigl{(}1+\delta({\bm{x}})\bigr{)}=\frac{m}{a^{3}}\sum_{n}\operatorname{\delta^{\textsc{d}}}({\bm{x}}-{\bm{x}}_{n}).$
(2.35)

The long- and short-range forces in (2.34) are split with a Gaussian kernel of
comoving width $x_{\mathrm{s}}$. One can verify the above force splitting,
using the fact that $e^{-k^{2}x_{\mathrm{s}}^{2}}/k^{2}$ and ${\rm
erf}(x/2x_{\mathrm{s}})/4\pi x$ are a 3D Fourier transform pair. The
acceleration due to the tree force is

$-\nabla_{\bm{x}}\phi^{\mathrm{T}}=-\frac{Gm}{a}\sum_{n}\frac{{\bm{x}}-{\bm{x}}_{n}}{|{\bm{x}}-{\bm{x}}_{n}|^{3}}\biggl{[}\operatorname{erfc}\Bigl{(}\frac{|{\bm{x}}-{\bm{x}}_{n}|}{2x_{\mathrm{s}}}\Bigr{)}+\frac{|{\bm{x}}-{\bm{x}}_{n}|}{x_{\mathrm{s}}\sqrt{\pi}}\exp\Bigl(-\frac{|{\bm{x}}-{\bm{x}}_{n}|^{2}}{4x_{\mathrm{s}}^{2}}\Bigr{missing})\biggr{]}.$
(2.36)

Now let us modify the above conventional TreePM algorithm for an
anisotropically expanding universe. The modified Poisson equation (2.6),

$(1+\Delta_{i})^{-2}\frac{\partial^{2}}{\partial x_{i}^{2}}\phi=4\pi
G{\bar{\rho}_{\mathrm{m}}}a^{2}(1+\Delta_{0})\delta,$ (2.37)

implies that the TreePM potentials in the presence of the DC modes should be

$\displaystyle\phi^{\mathrm{PM}}({\bm{p}})$ $\displaystyle=-4\pi
G{\bar{\rho}_{\mathrm{m}}}a^{2}(1+\Delta_{0})\frac{\delta({\bm{p}})}{k^{2}}e^{-p^{2}x_{\mathrm{s}}^{2}},$
$\displaystyle\phi^{\mathrm{T}}({\bm{x}})$
$\displaystyle=-\frac{Gm}{a}\sum_{n}\frac{a}{|{\bm{r}}-{\bm{r}}_{n}|}-\phi^{\rm
PM}({\bm{x}}),$ (2.38)

where we have introduced ${\bm{p}}$ as the local comoving wavevector111Though
this Eulerian ${\bm{p}}$ is technically different from the Lagrangian one in
Sec. 2.2, they are both Fourier conjugates to the local comoving coordinates,
therefore denoted by the same symbol., related to the global comoving
wavevector by

$k_{i}=\frac{p_{i}}{1+\Delta_{i}}.$ (2.39)

And recall $r_{i}=a_{i}x_{i}$ is the physical coordinates. As expected, both
the PM and tree forces above are manifestly isotropic in the physical or
global comoving coordinates. In App. C we provide the functional form of
$\phi^{\rm PM}({\bm{x}})$, which involves integral and no longer has the
simple form as in Eq. (2.34).

Note that there are two choices on the force-splitting scale: isotropic in
physical scales, i.e. anisotropic in local comoving scales, or isotropic in
local comoving scales, i.e. anisotropic in physical scales. In the former case
the equations become quite simple. However, we have found this choice
introduces numerical artifacts, e.g. on the second order responses to tides of
the halo abundance. This is probably because in this case the force-splitting
boundary in local comoving scales (i.e., the simulation coordinates) is
neither isotropic nor constant in time, and can interact with the anisotropic
PM force artifacts at grid scales. To avoid this problem we choose to split
force isotropically in local comoving scales, following Ref. [9]. The tree
acceleration now is

$-\frac{\partial\phi^{\mathrm{T}}}{\partial
x_{i}}=-\frac{Gm}{a}\sum_{n}\frac{aa_{i}[{\bm{r}}-{\bm{r}}_{n}]_{i}}{|{\bm{r}}-{\bm{r}}_{n}|^{3}}+\partialderivative{\phi^{\rm
PM}}{x_{i}}.$ (2.40)

Since $\partialderivative*{\phi^{\rm PM}}{x_{i}}$ is computationally expensive
to exactly evaluate in simulations, we expanded this in terms of $\Delta_{i}$
and included up to the second order terms in $\Delta_{i}$ as in Ref. [9].
Details are given in App. C.

###  2.4 Time integration

From Eq. (2.5), the equation of motion for the peculiar part takes a simple
form of

$\frac{\dot{P}_{i}}{m}=-\frac{\partial\phi}{\partial x^{i}},$ (2.41)

where ${\bm{P}}=a^{2}m\dot{\bm{x}}$ is the canonical momentum of the $N$-body
Hamiltonian

$H=\sum_{n}\frac{{\bm{P}}_{n}^{2}}{2ma^{2}}+\sum_{n\neq
n^{\prime}}\frac{m^{2}\varphi({\bm{x}}_{n}-{\bm{x}}_{n^{\prime}})}{2a},$
(2.42)

in which $\varphi$ is the potential of a unit-mass particle in a box of
comoving size $L$ at $a=1$:

$\nabla_{\bm{x}}^{2}\varphi=4\pi
G\Bigl{(}\operatorname{\delta^{\textsc{d}}}({\bm{x}}-{\bm{x}}^{\prime})-\frac{1}{L^{3}}\Bigr{)}.$
(2.43)

Note that we have made the time dependence in the potential explicit by
introducing $\varphi$, which otherwise is simply related to $\phi$ by Eq.
(2.6) and Eq. (2.35).

$N$-body simulations use the computed gravitational forces to update the
particle velocities and then them to evolve the particle positions in time.
This time integration is performed using the kick and drift leapfrog operators
[32, 33]

$\displaystyle\mathrm{Kick}$
$\displaystyle:\quad\frac{{\bm{P}}}{m}\to\frac{{\bm{P}}}{m}-\nabla_{\bm{x}}\phi\int_{t}^{t+\Delta
t}\frac{\mathrm{d}t}{a},$ (2.44) $\displaystyle\mathrm{Drift}$
$\displaystyle:\quad{\bm{x}}\to{\bm{x}}+\frac{{\bm{P}}}{m}\int_{t}^{t+\Delta
t}\frac{\mathrm{d}t}{a^{2}}.$ (2.45)

On the anisotropic background, the Hamiltonian leading to the modified EoM Eq.
(2.41) is

$H=\sum_{n}\sum_{i=1}^{3}\frac{P_{n,i}^{2}}{2ma_{i}^{2}}+\sum_{n\neq
n^{\prime}}\frac{m^{2}\varphi({\bm{x}}_{n}-{\bm{x}}_{n^{\prime}})}{2a},$
(2.46)

where $P_{n,i}\equiv a_{i}^{2}m\dot{x}_{n,i}$ are the conjugate momenta of
$x_{n,i}$. Thus we have

$\displaystyle\mathrm{Kick}$
$\displaystyle:\quad\frac{{\bm{P}}}{m}\to\frac{{\bm{P}}}{m}-\nabla_{{\bm{x}}}\phi\int_{t}^{t+\Delta
t}\frac{dt}{a},$ (2.47) $\displaystyle\mathrm{Drift}$ $\displaystyle:\quad
x_{i}\to x_{i}+\frac{P_{i}}{m}\int_{t}^{t+\Delta t}\frac{dt}{a_{i}^{2}}.$
(2.48)

Note that the drift operator changes due to the modified canonical momentum,
whereas the kick operator is unchanged because the anisotropic effect on force
calculations is already accounted in Sec. 2.3. We implemented these
modifications of the TreePM forces and time integration based on L-Gadget2
[33].

##  3 Simulations

type of simulations | $\Delta_{0}^{(1)}$ | $\Delta_{\rm p}$ | $\Delta_{\rm e}$ | $L~{}[{\rm Mpc}/h]$ | $N_{\rm p}$ | realizations
---|---|---|---|---|---|---
|  |  |  | 250 | $1024^{3}$ | 6
fiducial | 0 | 0 | 0 | 1000
|  |  |  | 3000
|  |  |  | 250
$\Delta_{0}$-type | $\pm 0.09$ | 0 | 0 | 1000
|  |  |  | 3000
|  |  |  | 250
$\Delta_{\rm p}$-type | 0 | $\pm 0.15$ | 0 | 1000
|  |  |  | 3000
|  |  |  | 250
$\Delta_{\rm e}$-type | 0 | 0 | $\pm 0.1$ | 1000
|  |  |  | 3000
Table 1: Summary of our 126 $N$-body simulations, of different types of mean
strains (related to the DC density and tidal modes) in the box. The
$\Delta_{0}$-type simulations have isotropic strain (the DC density
perturbation); the $\Delta_{\mathrm{p}}$-type and $\Delta_{\mathrm{e}}$-type
perturb the two DC tidal modes separately as decomposed in Eq. (3.1); and the
fiducial type are conventional simulations without any DC modes on global
background cosmology. $\Delta_{0}^{(1)}$, $\Delta_{\mathrm{p}}$, and
$\Delta_{\mathrm{e}}$ specify the strain amplitudes of their respective types;
$L$ are the simulation box sizes; and $N_{\rm p}$ is the total number of
particle.

We perform $N$-body simulations in the tidal backgrounds as described in the
previous section, using the Planck 2015 cosmology [34]:
${\Omega_{\mathrm{m}}}=0.3089$, $\Omega_{\Lambda}=0.6911$, $h=0.6774$,
$n_{s}=0.9667$, and $\sigma_{8}=0.8159$. We generate initial conditions at
redshift $z_{\mathrm{i}}=49$ with $1024^{3}$ particles using CLASS [35] and
our modified 2LPTIC 222https://cosmo.nyu.edu/roman/2LPT/, and then run the
modified L-Gadget2 [33] with $N_{\rm mesh}=2048^{3}$ TreePM grid. To cover a
large range of halo masses, we use 3 different box sizes:
$L=250~{}\mathrm{Mpc}/h$, $L=1~{}\mathrm{Gpc}/h$, and $L=3~{}\mathrm{Gpc}/h$.
For each type of simulations, we set the force softening scales to $4\%$ of
the mean particle distances, which means $9.77~{}{\rm kpc}/h$, $39.1~{}{\rm
kpc}/h$, and $117~{}{\rm kpc}/h$ for $L=250~{}\mathrm{Mpc}/h$,
$L=1~{}\mathrm{Gpc}/h$, and $L=3~{}\mathrm{Gpc}/h$ simulations, respectively.
About the accuracy of the time integration, ErrTolIntAccuracy is set to 0.05.
These parameters on the force accuracy potentially have an impact on the halo
shape as recently presented in Ref. [36].

So far we have used $\Delta_{1}$, $\Delta_{2}$, and $\Delta_{3}$ to
parametrize the remaining three eigenvalues of the background strain.
Alternatively one can parametrize it by one isotropic dilation and two
anisotropic scalings [6], with the former equal to the negative DC overdensity
at linear order,

$\displaystyle-\Delta^{(1)}_{0}$
$\displaystyle=\Delta^{(1)}_{1}+\Delta^{(1)}_{2}+\Delta^{(1)}_{3}=-(\tau_{1}+\tau_{2}+\tau_{3}),$
$\displaystyle\Delta_{\mathrm{e}}$
$\displaystyle\equiv\frac{\Delta^{(1)}_{1}-\Delta^{(1)}_{2}}{2}=-\frac{\tau_{1}-\tau_{2}}{2},$
$\displaystyle\Delta_{\mathrm{p}}$
$\displaystyle\equiv\Delta^{(1)}_{3}-\frac{\Delta^{(1)}_{1}+\Delta^{(1)}_{2}}{2}=-\tau_{3}+\frac{\tau_{1}+\tau_{2}}{2}.$
(3.1)

The subscripts of the two anisotropic scaling modes stand for ellipticity and
prolaticity333We do not impose ordering on $\Delta_{i}$’s to compress the
parameter space (cf. [6])..

To focus on tidal effects, we can remove the linear order overdensity mode by
setting $\Delta_{0}^{(1)}=0$, and run “pure” tidal simulations by perturbing
the remaining two degrees of freedom, $\Delta_{\rm p}$ and $\Delta_{\rm e}$.
We also vary the two modes separately, further dividing those pure tidal
simulations into $\Delta_{\mathrm{p}}$\- and $\Delta_{\mathrm{e}}$-type
simulations. For $\Delta_{\mathrm{p}}$-type simulations we chose the
background strain to satisfy $\Delta^{(1)}_{0}=0$ and $\Delta_{\rm e}=0$ so
that the remaining degree of freedom is only $\Delta_{\rm
p}=3\Delta^{(1)}_{3}/2$. For $\Delta_{\mathrm{e}}$-type simulations we chose
the background strain to satisfy $\Delta^{(1)}_{0}=0$ and $\Delta_{\rm p}=0$
so that the remaining degree of freedom is only $\Delta_{\rm
e}=\Delta^{(1)}_{1}$. This means the configuration of the tides and linear
strain takes the form $(\tau_{1},\tau_{2},\tau_{3})=(1,1,-2)\Delta_{\rm p}/3$
for $\Delta_{\mathrm{p}}$-type simulations, and
$(\tau_{1},\tau_{2},\tau_{3})=(-1,1,0)\Delta_{\rm e}$ for
$\Delta_{\mathrm{e}}$-type simulations.

Our modification for the tidal background also works for the isotropic
configuration, where $\Delta_{\rm p}=\Delta_{\rm e}=0$ and only
$\Delta_{0}^{(1)}$ is varied. In fact this configuration is equivalent to the
usual isotropic separate universe simulation at the nonlinear level. We
validate our simulation pipeline by running the $\Delta_{0}$-type simulations
and comparing them to the usual separate universe simulations. We find their
power spectrum responses and halo biases are fully consistent with each other,
as summarized in App. D.

As reference and normalization for the response estimations, we also run some
conventional fiducial simulations without any variation in the DC density or
tidal modes. For each of type of simulations, each box size, and each sign
(positive, negative, or fiducial) of the DC modes, we have performed six
simulations, and in total used 126 simulations to produce our results.
Simulations with different DC modes shares the same random phases in the
initial conditions so that the sample variances in the response functions are
suppressed. We summarize all our simulations in Table 1.

##  4 Results

In this section, we present analytic expressions and numerical calibrations of
the power spectrum tidal response using 2LPT and $N$-body simulations for high
and low redshift, respectively; we also present measurements of the halo shape
response to the tidal field. For the latter, we study its dependence on the
halo mass and other (secondary) halo properties, including the concentration
and axis ratio. This is analogous to the halo assembly bias, the dependence of
halo abundance on variables beyond the halo mass.

###  4.1 Power spectrum responses

The long-wavelength modes modulate the evolution of the short-wavelength ones,
in their amplitudes and scales, known as the growth and dilation effects [1],
respectively. The growth effect focuses on the changes in the short mode
amplitude in the local comoving coordinate system, as considered in Sec. 2.2.
On the other hand, the dilation effect arises from the (anisotropic) expansion
of the local comoving coordinate system with respect to the global one, as
described in Sec. 2.1. In addition to the growth and dilation effects, the
change in the reference mean density by $(1+\Delta_{0})$ in power spectrum
estimation also contributes to its responses, as shown below.

Let us derive analytically this separation of the growth and dilation effects
in the power spectrum response functions. As in Sec. 2.3, we use ${\bm{p}}$ to
denote a wavevector in the local comoving coordinates, which is related to the
global comoving wavevector ${\bm{k}}$ by Eq. (2.39). At linear order, the
density and tidal perturbations induce anisotropic responses in the monopole
and quadrupole. Therefore we define the total response functions
$R^{\mathrm{E}}$, including both growth and dilation effects, and the growth
only response functions $R^{\mathrm{L}}$, by

$\displaystyle-\frac{\partial\ln
P}{\partial\Delta^{(1)}_{i}}\bigg{|}_{{\bm{k}},{\bm{\Delta}}=0}$
$\displaystyle\equiv
R^{\mathrm{E}}_{0}(k)+\operatorname{\mathcal{L}}_{2}(\hat{k}_{i})R^{\mathrm{E}}_{2}(k),$
$\displaystyle-\frac{\partial\ln
P_{W}}{\partial\Delta^{(1)}_{i}}\bigg{|}_{{\bm{p}},{\bm{\Delta}}=0}$
$\displaystyle\equiv
R^{\mathrm{L}}_{0}(p)+\operatorname{\mathcal{L}}_{2}(\hat{p}_{i})R^{\mathrm{L}}_{2}(p),$
(4.1)

where $P_{W}$ is the power spectrum measured in local comoving space. The
derivatives with respect to $\Delta_{i}$ is taken holding the other two
$\Delta_{j}$’s ($j\neq i$) fixed, at ${\bm{\Delta}}=0$ with ${\bm{\Delta}}$
being $\Delta_{0}$ and $\Delta_{i}$’s. The function
$\operatorname{\mathcal{L}}_{2}$ is the second order Legendre polynomial.

Our notation for the power spectrum responses bear a resemblance to that for
the halo biases. The linear halo bias can also be measured as responses of the
halo abundance with respect to the long-wavelength density mode. One measures
the Eulerian bias $b^{\mathrm{E}}$ or Lagrangian bias $b^{\mathrm{L}}$,
depending on whether the measurement is carried out in either global or local
comoving space (see e.g. [4]), respectively. Therefore this “Lagrangian”
superscript in $R^{\mathrm{L}}$ is not to be confused with _Lagrangian_
perturbation theory in Sec. 2.2.

Because the variances should be conserved when transforming between Fourier-
space volume elements,

$P({\bm{k}};{\bm{\Delta}})\,\mathrm{d}^{3}{\bm{k}}=(1+\Delta_{0})^{2}P_{W}({\bm{p}};{\bm{\Delta}})\,\mathrm{d}^{3}{\bm{p}},$
(4.2)

where the $(1+\Delta_{0})^{2}$ factor on the right hand side is due to a
change of the reference density in $P_{W}$. Therefore, for the dimensionless
power spectra we have

$k^{3}P({\bm{k}};{\bm{\Delta}})=(1+\Delta_{0})^{2}p^{3}P_{W}({\bm{p}};{\bm{\Delta}}).$
(4.3)

Using Eq. (2.39) and the chain rule we can relate the total and growth
responses in Eq. (4.1) by

$\displaystyle-\frac{\partial\ln
k^{3}P}{\partial\Delta^{(1)}_{i}}\bigg{|}_{\bm{k}}=2-\frac{\partial\ln
p^{3}P_{W}}{\partial\Delta^{(1)}_{i}}\bigg{|}_{\bm{k}}$
$\displaystyle=2-\frac{\partial\ln
P_{W}}{\partial\Delta^{(1)}_{i}}\bigg{|}_{\bm{p}}-\frac{\partial\ln
p^{3}P_{W}}{\partial p_{j}}\bigg{|}_{{\bm{\Delta}}}\frac{\partial
p_{j}}{\partial\Delta^{(1)}_{i}}\bigg{|}_{{\bm{k}}}$
$\displaystyle=2-\frac{\partial\ln
P_{W}}{\partial\Delta^{(1)}_{i}}\bigg{|}_{\bm{p}}-\hat{k}_{i}^{2}\frac{\mathrm{d}\ln
k^{3}P}{\mathrm{d}\ln k},$ $\displaystyle=2-\frac{\partial\ln
P_{W}}{\partial\Delta^{(1)}_{i}}\bigg{|}_{\bm{p}}-\frac{1}{3}\frac{\mathrm{d}\ln
k^{3}P}{\mathrm{d}\ln
k}-\operatorname{\mathcal{L}}_{2}(\hat{k}_{i})\frac{2}{3}\frac{\mathrm{d}\ln
k^{3}P}{\mathrm{d}\ln k},$ (4.4)

where $\hat{k}_{i}=k_{i}/k$ and ${\bm{\Delta}}=0$ is assumed.

Hence, the total (Eulerian) and growth (Lagrangian) responses are related with
each other through

$\displaystyle R^{\mathrm{E}}_{0}$
$\displaystyle=2+R^{\mathrm{L}}_{0}-\frac{1}{3}\frac{\mathrm{d}\ln
k^{3}P}{\mathrm{d}\ln k},$ $\displaystyle R^{\mathrm{E}}_{2}$
$\displaystyle=0+R^{\mathrm{L}}_{2}-\frac{2}{3}\frac{\mathrm{d}\ln
k^{3}P}{\mathrm{d}\ln k}.$ (4.5)

As explained in the beginning of this subsection, we have decomposed the total
responses $R^{\mathrm{E}}$ into three contributions: the constant two due to
the change of reference density, the growth responses $R^{\mathrm{L}}$, and
the dilation term proportional to the slope of the dimensionless power
spectrum. Note the quadrupole response is not affected by the reference
density, which is emphasized above with 0.

In the linear regime $P_{W}\propto D_{W}^{2}$ and we have calculated the
modified growth factor $D_{W}$ in Sec. 2.2. Rewriting Eq. (2.31) in Legendre
polynomial,

$D_{W}=D\Bigl{(}1+\frac{13}{21}\Delta_{0}-\frac{8}{21}\operatorname{\mathcal{L}}_{2}(\hat{p}_{i})\Delta^{(1)}_{i}\Bigr{)},$
(4.6)

so that the tree-level response functions are

$\displaystyle R^{\mathrm{L}}_{0}=\frac{26}{21},\qquad$ $\displaystyle
R^{\mathrm{E}}_{0}=\frac{68}{21}-\frac{1}{3}\frac{\mathrm{d}\ln
k^{3}P}{\mathrm{d}\ln k};$ $\displaystyle
R^{\mathrm{L}}_{2}=\frac{16}{21},\qquad$ $\displaystyle
R^{\mathrm{E}}_{2}=\frac{58}{21}-\frac{2}{3}\frac{\mathrm{d}\ln
k^{3}P}{\mathrm{d}\ln k}.$ (4.7)

These results are consistent with the standard perturbation theory
calculations, e.g. [12, 13].

In the nonlinear regime, we need to numerically calibrate the response
functions with simulations. We focus on the growth response because the other
contributions in Eq. (4.5) are well understood, with the dilation term readily
computable from the nonlinear power spectrum. A pure long tidal mode has only
two degrees of freedom, as parametrized in Eq. (3.1) by $\Delta_{\mathrm{p}}$
and $\Delta_{\mathrm{e}}$. One can show that they modulate the power spectrum
respectively by

$\displaystyle\delta P_{W}({\bm{p}};\Delta_{\mathrm{p}})$
$\displaystyle=-P(p)R^{\mathrm{L}}_{2}(p)\Delta_{\mathrm{p}}\operatorname{\mathcal{L}}_{2}(\hat{p}_{3}),$
$\displaystyle\delta P_{W}({\bm{p}};\Delta_{\mathrm{e}})$
$\displaystyle=-P(p)R^{\mathrm{L}}_{2}(p)\Delta_{\mathrm{e}}\bigl{[}\operatorname{\mathcal{L}}_{2}(\hat{p}_{1})-\operatorname{\mathcal{L}}_{2}(\hat{p}_{2})\bigr{]}.$
(4.8)

In practice, we measure $R_{2}^{\mathrm{L}}$ from those
$\Delta_{\mathrm{p}}$-type and $\Delta_{\mathrm{e}}$-type of simulations.

It is more straightforward to extract $R_{2}^{\mathrm{L}}$ from the
$\Delta_{\mathrm{p}}$-type simulations. Eq. (4.8) implies that the relevant
changes in $P_{W}$ lie in its quadrupole along the $z$-axis (denoted by
$\ell_{\mathrm{p}}$ below),

$P^{\ell_{\mathrm{p}}=2}_{W}(p;\Delta_{\rm p})\equiv
5\int\frac{\mathrm{d}^{2}\hat{\bm{p}}}{4\pi}P_{W}({\bm{p}};\Delta_{\rm
p})\operatorname{\mathcal{L}}_{2}(\hat{p}_{3})=-P(p)R_{2}^{\rm
L}(p)\Delta_{\rm p},$ (4.9)

due to the orthonormality of the Legendre polynomials and the fact that $P(p)$
is istropic. Therefore, the estimator for $R_{2}^{\rm L}(p)$ can be
constructed as

$R_{2}^{\rm
L}(p,z)=-\frac{P_{W}^{\ell_{\mathrm{p}}=2}(p,z;\Delta_{\mathrm{p}}=+\epsilon)-P_{W}^{\ell_{\mathrm{p}}=2}(p,z;\Delta_{\mathrm{p}}=-\epsilon)}{2\epsilon
D(z)P(p,z)}.$ (4.10)

Similarly, we can also estimate the growth response from the power spectrum
quadrupoles along both $x$ and $y$ axes, from the $\Delta_{\mathrm{e}}$-type
simulations,

$P^{\ell_{\mathrm{e}}=2}_{W}(p;\Delta_{\rm e})\equiv
5\int\frac{\mathrm{d}^{2}\hat{\bm{p}}}{4\pi}P_{W}({\bm{p}};\Delta_{\rm
e})\bigl{[}\operatorname{\mathcal{L}}_{2}(\hat{p}_{1})-\operatorname{\mathcal{L}}_{2}(\hat{p}_{2})\bigr{]}=-P(p)R_{2}^{\rm
L}(p)\Delta_{\rm e},$ (4.11)

where we have used the orthonormality of the Legendre polynomials and

$\int\frac{\mathrm{d}^{2}\hat{\bm{p}}}{4\pi}\operatorname{\mathcal{L}}_{2}(\hat{p}_{1})\operatorname{\mathcal{L}}_{2}(\hat{p}_{2})=-\frac{1}{10}.$
(4.12)

Thus, we can estimate $R_{2}^{\rm L}(p)$ by

$R_{2}^{\rm
L}(p,z)=-\frac{P_{W}^{\ell_{\mathrm{e}}=2}(p,z;\Delta_{\mathrm{e}}=+\epsilon)-P_{W}^{\ell_{\mathrm{e}}=2}(p,z;\Delta_{\mathrm{e}}=-\epsilon)}{2\epsilon
D(z)P(p)}.$ (4.13)

####  4.1.1 Responses at high redshifts from 2LPT

At high redshifts 2LPT works accurately with more modes in the linear regime.
Using results derived in Sec. 2.2, we have modified the initial condition
generator to incorporate the leading order impact of the long modes to 2LPT.
This allows us to measure the power spectrum responses reliably at high
redshifts. For this purpose, we generated eight pairs of 2LPT realizations
with $\Delta_{\rm p}$-type which contain $1024^{3}$ particles in $100~{}{\rm
Mpc}/h$ boxes at $z=49,15,10,\textrm{and }7$.

We first compare the isotropic response $R^{\mathrm{L}}_{0}$ measured from our
2LPT to that from the usual separate universe simulations, in order to find
out the scale above which our 2LPT responses converge. Fig. 9 in App. D shows
the results of this convergence test. The isotropic 2LPT responses are
accurate to within $5\%$ of the separate universe results, for $k\leq
9~{}h/{\rm Mpc}$ at $z=15$, $k\leq 4~{}h/{\rm Mpc}$ at $z=10$, and $k\leq
2~{}h/{\rm Mpc}$ at $z=7$.

We measure the 2LPT response $R_{2}^{\mathrm{L}}$ at each redshift, and show
the results in Fig. 1. The valid scale obtained from the comparison of
$R_{0}^{\rm L}$ are shown as vertical dotted lines. Note that we have scaled
the vertical axis for easy comparison with previous studies, $3R^{\rm
L}_{2}(p)/2=G(k)$ used in Refs. [9, 10, 11]. It is clear that at these high
redshifts $R_{2}^{\rm L}(p)$ grows more than the prediction of the tree-level
perturbation theory on small scales for $k\geq 2~{}h/{\rm Mpc}$. While the
overall trend is in agreement with the results of Ref. [11], our responses
have quantitatively less enhancement than that in Ref. [11], even after taking
into account the possible error of 2LPT. This may be attributed to the
difference on details of the implementations, e.g. generating initial
conditions at different orders, and needs the further examination.

Figure 1:  Power spectrum growth response function to the large-scale tidal
field, $R_{2}^{\mathrm{L}}$, measured at high redshifts from 2LPT
realizations. The shaded region shows the $1\sigma$ error. At $z=49$, the
response agrees well with the linear perturbation theory prediction on all
scales, which is shown by the horizontal dashed line, while at later times the
responses are enhanced on small scales. Vertical dotted lines mark for each
redshift the convergence scale, above which isotropic 2LPT responses
$R_{0}^{\mathrm{L}}$ converge to within $5\%$ of the validation separate
universe simulations. Though of similar trend, our results show quantitatively
less enhancement than that in Ref. [11], even after taking into account the
possible error of 2LPT.

####  4.1.2 Responses at low redshifts from $N$-body

Although Refs. [9, 10] already investigated the linear tidal response of a
matter power spectrum from their simulations, we measure the growth response
$R_{2}^{\rm L}(p)$ at low redshift from our $N$-body simulations to check the
validity of our numerical implementation. In this subsection we show results
from $1~{}{\rm Gpc}/h$ simulations to cover both linear and nonlinear regimes.

First, since we ran two different types of the tidal simulations ($\Delta_{\rm
p}$\- and $\Delta_{\rm e}$-type), we test whether the simulations with these
different types give converging results. In the left panel of Fig. 2 we
present $R_{2}^{\rm L}(p)$ from both $\Delta_{\rm p}$-type simulations (Eq.
(4.10)) and $\Delta_{\rm e}$-type simulations (Eq. (4.13)). Both results agree
well up to $k=10~{}h/{\rm Mpc}$. Therefore we combine $\Delta_{\rm e}$\- and
$\Delta_{\rm p}$-type simulations to estimate $R_{2}^{\rm L}(p)$ in the
results below.

The right panel of Fig. 2 shows $R_{2}^{\rm L}(p)$ for several redshifts:
$z=3,~{}2,~{}1,~{}0.5,$ and $0$. At these low redshifts, $R_{2}^{\rm L}(p)$
decreases over all the scales as the redshift decreases. This is likely
because strong nonlinearity tends to erase the memory of large-scale tidal
field gradually. Although these features are in general consistent with
previous studies, quantitatively there is a small difference. For instance, at
$z=3$ tidal response from our simulations takes a maximum $3R_{2}^{\rm
L}/2\sim 1.5$ while in Ref. [11] the maximum at $z=3$ is less than $1.5$. In
addition, Ref. [9] reports the different behaviour of $R_{2}^{\rm L}$ at $z=2$
from ours. These disagreements may arise from the halo sample variance in the
nonlinear regime, which is also seen in the separate universe simulations
(e.g. Ref. [1]) or the difference in details of the numerical implementations
and require further studies.

Figure 2: Power spectrum responses as in the previous figure, measured at low
redshifts from $N$-body simulations with $L=1~{}{\rm Gpc}/h$. The left plot
shows measurements from $\Delta_{\rm p}$\- and $\Delta_{\rm e}$-type
simulations at $z=0$, and the right one combines $\Delta_{\mathrm{p}}$ and
$\Delta_{\mathrm{e}}$ types at multiple low redshifts. The shaded region shows
the $1\sigma$ error. The horizontal dashed line is the tree-level perturbation
theory prediction. The vertical dotted line on the right panel depicts the
particle Nyquist wavenumber for $L=1~{}{\rm Gpc}/h$ simulations.

###  4.2 Halo shape response

The tidal field is known to cause intrinsic alignment of galaxy shapes.
Likewise, halo shapes respond to the large-scale tidal modes, with the
sensitivity captured by the shape bias. This is analogous to the halo bias
which is a response of the halo abundance to the large-scale density mode. In
this section, we present measurements of the shape bias using our tidal
simulations. We show its universal behavior as a function of the linear bias,
and its dependence on properties beyond mass, which we call the shape assembly
bias.

We identify dark matter halos from our simulations with AHF [37]. It uses
adaptively refined meshes to find local density peaks as centers of
prospective halos. It then defines the halos as spherical overdensity (SO)
regions $\Delta_{\mathrm{h}}$ times denser than the mean matter density
${\bar{\rho}_{\mathrm{m}}}$. We choose $\Delta_{\mathrm{h}}=200$, and disable
the gravitational unbinding procedure to find the host halos with more than
400 particles. While we need to identify SO halos in the global comoving
coordinates, AHF by default uses the simulation coordinates that are local
comoving, so that halos are not identified as the spherical overdensity in the
global ones. Therefore, we modify the AHF code to use the Euclidean metric for
the halo identification in the global coordinates.

Given an SO halo, we define its quadrupole shape in two ways, the inertia
tensor and the reduced inertia tensor. The former is defined as

$I_{ij}\equiv\sum_{i=1}^{N_{\mathrm{p}}}m_{\mathrm{p}}x_{i}x_{j},$ (4.14)

summing over all $N_{\mathrm{p}}$ particles of the halo. $m_{\mathrm{p}}$ is
the particle mass and $x_{i}$ is $i$-th components of the particle location
with respect to the halo center. In the literature $I_{ij}$ is sometimes
normalized by the halo mass, which however does not affect our response
measurement presented below.

The reduced inertia tensor is defined similarly but with additional radius
weighting

$J_{ij}\equiv\sum_{i=1}^{N_{\mathrm{p}}}m_{\mathrm{p}}\frac{x_{i}x_{j}}{x^{2}},$
(4.15)

where $x$ is the distance of a particle to the halo center. This definition
uses a dimensionless ratio and therefore weight each mass equally only by
angular position regardless of radial distance $x$. Compared to Eq. (4.14),
$J_{ij}$ upweights the inner masses and thus should be more strongly
correlated with properties of galaxies which reside in the halo.444See, e.g.,
Refs. [38, 39, 40] for the detection of large misalignments ($\sim 30~{}{\rm
deg}$) between the major axes of central galaxies and their host halos when
$I_{ij}$ is used to define halo shapes. Therefore, in the following main text,
we use $J_{ij}$ to estimate halo shapes. The results from $I_{ij}$ are
summarized in App. E. Ref. [41] gives detailed discussion on the dependence of
the halo shape on its definition. Ref. [42] also presented comparisons of
$I_{ij}$ and $J_{ij}$ in the context of intrinsic alignments of galaxies (see
their Appendix B).

According to the linear alignment model, at leading order halo shapes responds
to the external tidal field as

$J_{ij}=J_{0}\left[\frac{1}{3}\delta_{ij}^{\rm K}+b_{K}K_{ij}\right],$ (4.16)

where $J_{0}$ is the trace component of the shape tensor: $J_{0}={\rm
Tr}[J_{ij}]$, $K_{ij}$ is the DC tidal field:
$K_{ij}=(\partial_{i}\partial_{j}/\partial^{2}-\delta^{\rm K}_{ij}/3)\delta$
and $b_{K}$ is the dimensionless linear shape bias parameter, which is related
to the conventionally used linear alignment coefficient $C_{1}$ through
$b_{K}=-a^{3}C_{1}\bar{\rho}(a)/D(a)$. The shape bias $b_{K}$ represents the
strength of the response or alignment and thus an analogous parameter to the
linear bias $b_{1}$, which describes the response of the number density of
halos to the spherically symmetric long-wavelength perturbation: $\delta_{\rm
h}=b_{1}\delta$.

Following the decomposition of the traceless components of the background
strain into $\Delta_{\rm p}$ and $\Delta_{\rm e}$, it is convenient to define
the following two quantities

$\displaystyle J_{\rm p}$ $\displaystyle\equiv
J_{33}-\frac{J_{11}+J_{22}}{2},$ (4.17) $\displaystyle J_{\rm e}$
$\displaystyle\equiv\frac{J_{11}-J_{22}}{2}.$ (4.18)

Then, Eq. (4.16) implies

$\displaystyle b_{K}(M,z)$ $\displaystyle=-\frac{J_{\rm p}(M,z;\Delta_{\rm
p}=+\epsilon)-J_{\rm p}(M,z;\Delta_{\rm p}=-\epsilon)}{2\epsilon
D(z)J_{0}(M,z)},$ (4.19) $\displaystyle b_{K}(M,z)$
$\displaystyle=-\frac{J_{\rm e}(M,z;\Delta_{\rm e}=+\epsilon)-J_{\rm
e}(M,z;\Delta_{\rm e}=-\epsilon)}{2\epsilon D(z)J_{0}(M,z)},$ (4.20)

for $\Delta_{\rm p}$-type and $\Delta_{\rm e}$-type simulations, respectively.
Having these relations, we can estimate $b_{K}$ from our tidal simulations by
measuring the averaged trace of the shape tensor, $J_{0}$, and the averaged
traceless components $J_{\rm p}$ and $J_{\rm e}$ from fiducial, $\Delta_{\rm
p}$-type, and $\Delta_{\rm e}$-type simulations respectively.

Note that the linear shape bias $b_{K}$ measured in this way should be
regarded as the Lagrangian shape bias since we do not take into account the
volume distortion due to the background strain. For linear shape bias,
however, there is no difference between the Lagrangian shape bias $b_{K}^{\rm
L}$ and the Eulerian one $b_{K}^{\rm E}$ unlike the number density bias where
the Lagrangian linear bias $b^{\rm L}_{1}$ is related to the Eulerian one
$b^{\rm E}_{1}$ through $b_{1}^{\rm E}=b_{1}^{\rm L}+1$. This is because the
pure tidal field does not induce the volume distortion at linear order of the
tides and can be explicitly shown by considering the conservation laws:
$J^{\rm E}_{ij}({\bm{x}})d^{3}{\bm{x}}=J^{\rm L}_{ij}({\bm{q}})d^{3}{\bm{q}}$
and $(1+\delta^{\rm E}({\bm{x}}))d^{3}{\bm{x}}=d^{3}{\bm{q}}$555 For the
second order shape bias, the Lagrangian shape bias is no longer identical to
the Eulerian one. See the discussion in [43]. Thus, in this paper we do not
distinguish $b_{K}^{\rm L}$ from $b_{K}^{\rm E}$ and the linear shape bias is
just written as $b_{K}$.

####  4.2.1 Convergence on the resolution and external tides

Figure 3: Linear alignment coefficient, $b_{K}$, for the reduced inertial
tensor, $J_{ij}$, at $z=0$. The left plot shows measurements from simulations
of different box sizes, and the right one shows measurements from different
tidal types.

Before showing the redshift- and environment-dependence of $b_{K}$, here we
discuss the convergence of measured $b_{K}$ for different resolutions and
different kinds of tides.

The left panel of Figure 3 shows $b_{K}$ from different boxsize simulations,
meaning the different resolutions since we fix the number of particles. The
results are in agreement with each other over all the mass range except for
the 6th-8th mass bins where the $250~{}{\rm Mpc}/h$ and $1~{}{\rm Gpc}/h$
simulations give slightly different results. For the 6th and 7th mass bins,
these differences can be attributed to the insufficient number of particles in
the inner regime of halos in $1~{}{\rm Gpc}/h$ simulation to determine halo
shapes, given that the results from $I_{ij}$ are converged at these mass bins
(see Fig. 11 in App. E). On the other hand, at the 8th mass bin $b_{K}$ from
$250~{}{\rm Mpc}/h$ simulations are not in agreement with that from $1~{}{\rm
Gpc}/h$ for both $J_{ij}$ and $I_{ij}$ results. This could happen due to the
small number of halos at this mass bin in $250~{}{\rm Mpc}/h$ simulations.
Considering these results, in the following we use $250~{}{\rm Mpc}/h$,
$1~{}{\rm Gpc}/h$, and $3~{}{\rm Gpc}/h$ simulations for the 1st-7th, 8th-9th,
and 10th-11th mass bins, respectively.

In the right panel of Fig. 3 we show $b_{K}$ from different kinds of tides,
namely $\Delta_{\rm p}$-type and $\Delta_{\rm e}$-type tides. They are in good
agreement with each other, which implies that the validity of the linear
alignment model is irrelevant to the substructure of the cosmic web such as
knots, filaments, or pancakes. Because the results from the two different
tides are converged over all the mass range, in the following we combine two
kinds of simulations to estimate $b_{K}$.

####  4.2.2 Redshift-dependence: the relation between $b_{K}$ and $b_{1}^{\rm
E}$

Figure 4: Linear alignment coefficient, $b_{K}$, for the reduced inertial
tensor, $J_{ij}$, at various redshifts. Here we have combined results from
difference box sizes and different tides.

Here we discuss the redshift-dependence of the linear alignment coefficient
and show that there seems an universal relation between $b_{K}$ and $b_{1}$.

Fig. 4 shows the linear alignment coefficient, $b_{K}$, for different
redshifts. It is clear that the absolute value of $b_{K}$ is greater at more
massive halos and at higher redshift. This means that more massive halos align
stronger than less massive ones and the strength of the alignment becomes
larger as redshift increases for all mass range. These trends are supposed to
originate from the fact that the alignment of halo shape is also affected by
the surrounding matter distribution of each halos; less massive halos are
susceptible to their surroundings and as time evolves the impact of their
surroundings becomes greater.

These trends are similar to the linear bias $b_{1}$ so it is interesting to
explore the relation between $b_{K}$ and $b_{1}$. In Fig. 5 we plot $b_{K}$ as
a function of the Eulerian linear bias $b_{1}^{\rm E}$ for various redshifts.
$b_{1}^{\rm E}$ is estimated from the Lagrangian linear bias $b_{1}^{\rm L}$,
which is directly measured as the response of the halo number in our
$\Delta_{0}$-type simulations (see also App. D), using $b_{1}^{\rm
E}=b_{1}^{\rm L}+1$. We find the relation between $b_{K}$ and $b_{1}^{\rm E}$
shows an universal behaviour over the range $z=0\sim 2$. This universal
relation is also found when using $I_{ij}$ (see App. E) and thus is not
relevant to how to measure the halo shapes. This strongly suggests that
$b_{K}$ is also uniquely determined by some quantity depending mass as
$b_{1}^{\rm E}$ does by the variance of the dark matter density field. Since
our simulations enable both $b_{K}$ and $b_{1}$ to be measured very
accurately, here we provide the fitting formula in the form of
$b_{K}(b_{1}^{\rm E})$ for convenience. We combine results from all redshifts
and obtain a fitting formula of the $b_{K}$-$b_{1}$ relation using a very
simple rational function:

$b_{K}=\frac{0.06461-0.09322b_{1}^{\mathrm{E}}}{1+0.3073b_{1}^{\mathrm{E}}}$
(4.21)

This fitting function is shown as dashed curve in Fig. 5.

Figure 5: Linear alignment coefficient, $b_{K}$, as a function of the Eulerian
linear bias, $b_{1}^{E}$, combining all redshift and mass information. The
dashed curve is the fitting given in Eq. (4.21).

####  4.2.3 Secondary halo shape responses

Studies of the halo assembly bias show that the halo bias depends on
properties other than the mass. Likewise, it is natural that the halo shape
response also possesses rich dependences beyond the halo mass. Here we study
dependences of $b_{K}$ on the halo concentration and the eccentricity of
inertial tensor of halos.

We use the AHF halo finder to measure the halo concentration parameter $c$.
Instead of fitting a NFW halo profile to each halo, AHF measures the ratio
$v_{\mathrm{max}}/v_{200}$. $v_{\mathrm{max}}$ is the maximum circular
velocity, $v_{\mathrm{max}}\equiv\max\sqrt{GM(<R)/R}$, and $v_{200}$ is the
circular velocity at virial radius
$v_{\mathrm{200}}\equiv\sqrt{GM_{200}/R_{200}}$. Assuming a NFW halo profile,
$v_{\mathrm{max}}/v_{200}$ is related to $c$ as given by [44], and thus is
used in AHF to determine the concentration.

At each redshift and mass-bin we divided halo samples into those with greater
than the median of the concentration and those with smaller concentration.
Then we measured $b_{K}$ from each group. The results are shown in Fig. 6. For
all redshifts halos with lower concentration tend to have a large amplitude of
$|b_{K}|$ at high mass, while the difference is likely to become small at low
mass. Since halos with the high concentration are expected to be formed from
highly curved peaks [45], the process of collapse into halo is not expected to
be much affected by large-scale tidal field. Further, since they were formed
earlier they tend to lose their memory on large-scale tidal field through
interacting with the local surroundings for a long time.

Next we discuss a dependence of the axis-ratio of halo shapes on $b_{K}$. We
introduce the axis-ratio as the ratio of the major axis to the minor axis of
the shape tensor: $q\equiv J_{3}/J_{1}$ where $J_{1}$, $J_{2}$, and $J_{3}$
are the eigenvalues of $J_{ij}$ satisfying $J_{1}\geq J_{2}\geq J_{3}$. As is
done in the concentration, we divided halo samples into two groups: those
above the median of the axis-ratio and those below the median, and then
measured $b_{K}$ from each group. Fig. 7 presents the results. For all
redshifts and mass bins |$b_{K}$| from lower $q$ samples is greater than
higher $q$ samples, which means that halos with rounder shapes do not respond
to the large-scale tidal field as strongly as halos with greater ellipticity.
This implies that distortion of halo shapes is indeed accelerated by the
large-scale tidal field. The same trend is found for a study of intrinsic
alignments in Ref. [40]: more elongated halos are more tightly aligned with
the surrounding matter distribution. Thus, the existence of the DC mode would
not only bias measurements of the cosmic shear power spectrum in weak lensing
surveys but also affect the cosmological application of intrinsic alignment
itself.

Figure 6: Upper panels: The halo concentration dependence of the linear
alignment coefficient, $b_{K}$, at various redshifts and masses: $b_{K}$ from
high concentration (red) and low concentration (blue). Lower panels: The
median concentration by which we divided halo samples.

Figure 7:  Upper panels: The axis-ratio ($q=J_{3}/J_{1}$ with $J_{1}\geq
J_{3}$) dependence of the linear alignment coefficient, $b_{K}$, at various
redshifts and masses: $b_{K}$ from high $q$ (low ellipticity, blue) and low
$q$ (high ellipticity, red). Lower panels: The median axis-ratio on which we
divided halo samples.

##  5 Discussion

In this paper, we have implemented mean (DC) tidal and density fluctuations
into cosmological $N$-body simulations by absorbing them in effective
anisotropic background expansion. We have improved upon previous works [8, 9,
10, 11] on generating initial conditions for the tidal simulations, with full
second order Lagrangian dynamics properly solved in general anisotropic
background.

The 2LPT in anisotropic background can be used to investigate the linear tidal
response of the matter power spectrum, $R_{2}^{\rm L}$, on quasi-nonlinear
scales. Since in simulations $R_{2}^{\rm L}$ at high redshifts suffer from
various numerical artifacts [9], it is of particular interest to measure
$R_{2}^{\rm L}$ on small scales at high redshifts from the 2LPT. We have found
$R_{2}^{\mathrm{L}}$ is enhanced upon the tree-level perturbation theory
prediction at $z=7\sim 15$. Our 2LPT results should be more robust since they
are not affected by high-redshift numerical artifacts that affect $N$-body
simulations. Though our findings are qualitatively consistent with the result
of Ref. [11], our enhancement is quantitatively weaker than that in Ref. [11].

We have also investigated the effect of large-scale tidal field on three-
dimensional halo shapes using our simulations. The linear alignment model
predicts that the halo shape responds to the large-scale tidal field and thus
linearly related with each other: $J_{ij}=b_{K}K_{ij}$. Our tidal simulations
allow us to directly test this relation and measure the proportional
coefficient $b_{K}$ accurately. We have found that the dependence of $b_{K}$
on redshifts and halo mass is similar to that of the linear halo bias $b^{\rm
E}_{1}$; i.e., at the same halo mass $|b_{K}|$ is getting smaller as redshift
decrease and more massive halos have greater $|b_{K}|$. Furthermore, we have
noticed that the relation between $b_{K}$ and $b_{1}^{\rm E}$ shows the
universal behaviour over a wide range of redshifts and masses; $z=0\sim 2$ and
$M_{\rm h}=10^{12}\sim 10^{15}~{}M_{\odot}/h$. This implies that we can
construct an analytical, physical model that can properly describe the mass-
and redshift-dependence of $b_{K}$ as done for the linear bias using the peak
theory and excursion set approach. This kind of the theoretical prediction on
$b_{K}$ is quite useful especially when treating the intrinsic alignments as
the signal. In particular, such a model is of crucial importance for exploring
the angular-dependent (quadrupolar) primordial non-Gaussianity (PNG) from
observations of the intrinsic alignments. It is pointed out that the intrinsic
alignments can uniquely probe the quadrupolar PNG in Refs. [23, 24]. However,
since the bias induced by the quarupolar PNG on halo or galaxy shapes is
completely degenerated with the quadrupolar PNG signal, the lack of such a
model makes it impossible to extract the information on the quadrupolar PNG
from measurements of the intrinsic alignments. Therefore to develop a theory
for the linear alignment coefficient $b_{K}$, analogous to the linear bias
$b_{1}$ case, is an urgent issue and worth exploring in future works.

In addition, we have measured for the first time the secondary dependence of
$b_{K}$ on halo properties other than the mass; it also depends on the halo
concentration and axis-ratio. This can be seen as the shape “assembly bias” as
in the case of the number density bias 666Refs. [46, 47] discuss the impact of
the intrinsic alignment or galaxy shape on the density tracer as the assembly
bias. This should be distinguished with ours that is the assembly bias of the
intrinsic alignment.. These findings will help to understand how halo shapes
are determined in the hierarchical structure formation.

## Acknowledgments

KA acknowledge support from JSPS Research Fellowship for Young Scientists and
JSPS KAKENHI Grant Numbers JP19J12254 and JP19H00677. YL acknowledge support
from Fellowships at Simons Foundation, the Kavli IPMU established by World
Premier International Research Center Initiative (WPI) of the MEXT Japan, and
at the Berkeley Center for Cosmological Physics. TO acknowledges support from
the Ministry of Science and Technology of Taiwan under Grants No. MOST 109-
2112-M-001-027- and the Career Development Award, Academia Sinica (AS-
CDA-108-M02) for the period of 2019 to 2023.

##  Appendix A Solving DC density mode at second order

At second order, Eq. (2.12) gives

$\ddot{\Delta}_{i}^{(2)}+2H\dot{\Delta}_{i}^{(2)}=4\pi
G{\bar{\rho}_{\mathrm{m}}}\Bigl{(}{\Delta_{i}^{(1)}}^{2}-\frac{1}{3}\Delta_{0}^{(2)}\Bigr{)}.$
(A.1)

Also expand Eq. (2.9) at second order and use Eq. (2.14) to derive

$\Delta_{0}^{(2)}=-\sum_{i}\Delta_{i}^{(2)}+\sum_{i<j}\Delta_{i}^{(1)}\Delta_{j}^{(1)}+\sum_{i}{\Delta_{i}^{(1)}}^{2}=-\sum_{i}\Delta_{i}^{(2)}+\frac{1}{2}{\Delta_{0}^{(1)}}^{2}+\frac{1}{2}\sum_{i}{\Delta_{i}^{(1)}}^{2}.$
(A.2)

Sum over $i$ of the first equation and plug in the second one, one can show
that at matter dominated era

$\sum_{i}\Delta_{i}^{(2)}=\frac{3}{14}\sum_{i}{\Delta_{i}^{(1)}}^{2}-\frac{3}{14}{\Delta_{0}^{(1)}}^{2}.$
(A.3)

Thus

$\Delta_{0}^{(2)}=\frac{5}{7}{\Delta_{0}^{(1)}}^{2}+\frac{2}{7}\sum_{i}{\Delta_{i}^{(1)}}^{2}=\frac{17}{21}{\Delta_{0}^{(1)}}^{2}+\frac{2}{7}\sum_{i}\tau_{i}^{2},$
(A.4)

where the last equality follows from Eq. (2.13).

##  Appendix B Second order Lagrangian perturbation theory in an anisotropic
background

For 2LPT, the Jacobian determinant and matrix inverse in the master equation
Eq. (2.19) can be expanded as

$\displaystyle\Bigl{|}\frac{\partial{\bm{x}}}{\partial{\bm{q}}}\Bigr{|}\simeq
1+\sum_{i}\Psi_{i,i}^{(1)}+\sum_{i}\Psi_{i,i}^{(2)}+\frac{1}{2}\left[\sum_{i}\left(\Psi_{i,i}^{(1)}\right)^{2}-\sum_{ij}\Psi_{i,j}^{(1)}\Psi_{j,i}^{(1)}\right]$
(B.1)
$\displaystyle\left[\delta_{ij}+\Psi_{i,j}\right]^{-1}\simeq\delta_{ij}-\Psi_{i,j}^{(1)},$
(B.2)

leading to the second-order equation

$\sum_{i}\ddot{\Psi}^{(2)}_{i,i}+2\sum_{i}H_{i}\dot{\Psi}^{(2)}_{i,i}-\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)(1+\Delta_{0})\sum_{i}\Psi_{i,i}^{(2)}\\\
=-\sum_{i}\Psi_{i,i}^{(1)}\sum_{j}\left[\ddot{\Psi}^{(1)}_{j,j}+2H_{j}\dot{\Psi}_{j,j}^{(1)}\right]+\sum_{ij}\Psi_{i,j}^{(1)}\left[\ddot{\Psi}^{(1)}_{i,j}+2H_{i}\dot{\Psi}_{i,j}^{(1)}\right]\\\
+\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)(1+\Delta_{0})\left[\frac{1}{2}\sum_{i}\left(\Psi_{i,i}^{(1)}\right)^{2}-\frac{1}{2}\sum_{ij}\Psi_{i,j}^{(1)}\Psi_{j,i}^{(1)}\right].$
(B.3)

Similar to the linear order, we introduce the second order displacement
potential through $\Psi_{i}^{(2)}=\partial\psi_{W}^{(2)}/\partial
q_{i}\equiv\psi^{(2)}_{W,i}$. In the absence of the long modes, the equation
for $\psi^{(2)}$ reduces to

$\sum_{i}\ddot{\psi}^{(2)}_{,ii}+2H\sum_{i}\dot{\psi}^{(2)}_{,ii}-\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)\sum_{i}\psi^{(2)}_{,ii}=-\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)\left[\frac{1}{2}\left(\sum_{i}\psi^{(1)}_{,ii}\right)^{2}-\frac{1}{2}\sum_{ij}\psi^{(1)}_{,ij}\psi^{(1)}_{,ji}\right],$
(B.4)

where we used the linear equation Eq. (2.23). In this usual case, we denote
the time-dependent part of $\psi^{(2)}$ as $D^{(2)}(t)$, which obeys

$\displaystyle\ddot{D}^{(2)}+2H\dot{D}^{(2)}-\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)D^{(2)}=-\frac{3}{2}H^{2}\Omega_{\rm
m}{D^{(1)}}^{2}.$ (B.5)

In the matter-domination, we have $D^{(2)}=3{D^{(1)}}^{2}/7$. The correction
induced by the long modes, which is expressed by
$\epsilon^{(2)}(t,{\bm{q}})\equiv\psi_{W}^{(2)}(t,{\bm{q}})-\psi(t,{\bm{q}})$,
follows

$\displaystyle\sum_{i}\ddot{\epsilon}_{,ii}^{(2)}+$ $\displaystyle
2H\sum_{i}\dot{\epsilon}^{(2)}_{,ii}-\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)\sum_{i}\epsilon_{,ii}^{(2)}=$
$\displaystyle-2\sum_{i}\dot{\psi}^{(2)}_{,ii}\dot{\Delta}_{i}+2\sum_{ij}\dot{\Delta}_{i}\dot{\psi}_{,ij}^{(1)}\psi_{,ji}^{(1)}-\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)\sum_{j}\epsilon^{(1)}_{,jj}\sum_{i}\psi^{(1)}_{,ii}$
$\displaystyle+\sum_{ij}\psi^{(1)}_{,ij}\left[\ddot{\epsilon}_{,ij}^{(1)}+2H\dot{\epsilon}_{,ij}^{(1)}\right]$
$\displaystyle+\frac{3}{2}H^{2}{\Omega_{\mathrm{m}}}(a)\Delta_{0}\left[\sum_{i}\psi_{,ii}^{(2)}-\frac{1}{2}\sum_{i}\left(\psi^{(1)}_{,ii}\right)^{2}-\frac{1}{2}\sum_{ij}\psi^{(1)}_{,ij}\psi^{(1)}_{,ji}\right],$
(B.6)

where we have neglected $\mathcal{O}(\Delta_{i}^{2})$ terms. Notice that
$\epsilon^{(1)}$ is $\mathcal{O}(\Delta_{i})$. For the matter dominated era,
the solution for Eq. (B.6) is given by

$\displaystyle\sum_{i}\epsilon^{(2)}_{,ii}(t,{\bm{q}})=$
$\displaystyle\frac{1}{4}\sum_{i}\left[-\frac{16}{9}\psi^{(2)}_{,ii}(t,{\bm{q}})+\frac{8}{9}\sum_{j}\psi^{(1)}_{,ij}(t,{\bm{q}})\psi^{(1)}_{,ji}(t,{\bm{q}})\right]\Delta_{i}$
$\displaystyle+\frac{1}{6}\left[\sum_{i}\psi_{,ii}^{(2)}-\frac{1}{2}\sum_{i}\left(\psi^{(1)}_{,ii}\right)^{2}-\frac{1}{2}\sum_{ij}\psi^{(1)}_{,ij}\psi^{(1)}_{,ji}\right]\Delta_{0}$
$\displaystyle+\frac{1}{4}\left[-\frac{2}{3}\sum_{i}\psi^{(1)}_{,ii}(t,{\bm{q}})\sum_{j}\epsilon_{,jj}^{(1)}(t,{\bm{q}})+\frac{20}{9}\sum_{ij}\psi^{(1)}_{,ij}(t,{\bm{q}})\epsilon_{,ij}^{(1)}(t,{\bm{q}})\right].$
(B.7)

Although the modified second order growth factor, $D^{(2)}_{W}$, due to the
long modes can be identified as
$D^{(2)}_{W}(t,{\bm{p}})=D^{(2)}(1+\Delta_{0}/6-4\hat{p}_{i}^{2}\Delta_{i}/9)$,
the local gravitational tides cannot be neglected at second order.

##  Appendix C Force computation

In this appendix, we review how to evaluate the tree force, especially the
real-space counterparts of the PM force. $\phi^{\rm PM}({\bm{x}})$ satisfies

$\displaystyle\sum_{i}(1+\Delta_{i})^{-2}\partialderivative{x_{i}}^{2}\phi^{\rm
PM}({\bm{x}})=-4\pi
G\bar{\rho}_{m}a^{2}(1+\Delta_{0})\int\differential^{3}{\bm{x}}^{\prime}\delta({\bm{x}}^{\prime})\cdot\left[\frac{1}{8\pi\sqrt{\pi}x^{3}_{s}}\exp\left(-\frac{|{\bm{x}}-{\bm{x}}^{\prime}|^{2}}{4x_{s}^{2}}\right)\right],$
(C.1)

where the function inside the last bracket corresponds to the Fourier
transform of $\exp(-p^{2}x_{s}^{2})$, which is the Gaussian smoothing kernel
used in Eq. (2.38) to split force. The PM potential is related to the tree
potential as $\phi^{\rm PM}=\phi-\phi^{\rm T}$ and the solution for the PM
potential is found to be

$\displaystyle\phi^{\rm
PM}({\bm{x}})=-\frac{Gm}{2a\sqrt{\pi}x_{s}}\int_{0}^{\infty}\frac{\exp[-\frac{1}{4x_{s}^{2}}\left(\frac{(1+\Delta_{1})^{2}x^{2}}{(1+\Delta_{1})^{2}+\lambda}+\frac{(1+\Delta_{2})^{2}y^{2}}{(1+\Delta_{2})^{2}+\lambda}+\frac{(1+\Delta_{3})^{2}z^{2}}{(1+\Delta_{3})^{2}+\lambda}\right)]}{\sqrt{((1+\Delta_{1})^{2}+\lambda)((1+\Delta_{2})^{2}+\lambda)((1+\Delta_{3})^{2}+\lambda)}}~{}\differential\lambda,$
(C.2)

where we used $\bar{\rho}_{\rm
m}(1+\Delta_{0})\left[1+\delta({\bm{x}})\right]=m(1+\Delta_{0})/a^{3}\cdot\sum_{n}\delta^{\rm
D}({\bm{x}}-{\bm{x}}_{n})$. Although (C.2) has no closed analytic form, we can
approximate this potential by Taylor expansion in $\Delta_{i}$ 777Our
expansion is different from that in Ref. [9], where $\phi^{\rm PM}({\bm{x}})$
is expanded in $1+\Delta_{i}-\sqrt[3]{1+\Delta_{0}}$.. Using the following
identity

$I_{k}(r)\equiv\int_{0}^{\infty}\frac{\exp[-\frac{1}{4x_{s}^{2}}\frac{r^{2}}{a^{2}(1+\lambda)}]}{\sqrt{(1+\lambda)^{k}}}~{}\differential\lambda=\biggl{(}\frac{r}{2ax_{s}}\biggr{)}^{2-k}\gamma\biggl{(}\frac{k-2}{2},\frac{r^{2}}{4a^{2}x_{s}^{2}}\biggr{)}$
(C.3)

where $\gamma$ is the lower incomplete gamma function. We can express the
approximated $\phi^{\rm PM}({\bm{x}})$ up to the second order of $\Delta_{i}$
as

$\displaystyle\phi^{\rm PM}({\bm{r}})=$
$\displaystyle-\frac{Gm}{2a\sqrt{\pi}x_{s}}\left[I_{3}(r)+\sum_{i}\left(-I_{5}(r)+\frac{r_{i}^{2}}{2a^{2}x_{s}^{2}}I_{7}(r)\right)\Delta_{i}\right.$
$\displaystyle+\frac{1}{2}\sum_{i}\left(-I_{5}(r)+\left(3+\frac{r_{i}^{2}}{2a^{2}x_{s}^{2}}\right)I_{7}(r)-\frac{3r_{i}^{2}}{a^{2}x_{s}^{2}}I_{9}(r)+\frac{r_{i}^{4}}{4a^{4}x_{s}^{4}}I_{11}(r)\right)\Delta_{i}^{2}$
$\displaystyle\left.+\frac{1}{2}\sum_{i\neq
j}\left(I_{7}(r)-\frac{r_{i}^{2}+r_{j}^{2}}{2a^{2}x_{s}^{2}}I_{9}(r)+\frac{r_{i}^{2}r_{j}^{2}}{4a^{4}x_{s}^{4}}I_{11}(r)\right)\Delta_{i}\Delta_{j}\right]+{\cal
O}(\Delta_{i}^{3}).$ (C.4)

In order to derive the force from this potential, we must be careful that the
derivative should be taken with respect to the local comoving coordinate
${\bm{x}}$, not to ${\bm{r}}$. Thus, the force from the PM potential is
computed as

$\displaystyle\frac{\partial\phi^{\rm PM}}{\partial
x^{k}}=a(1+\Delta_{k})\frac{\partial\phi^{\rm PM}}{\partial r^{k}},$ (C.5)

and

$\displaystyle\partialderivative{\phi^{\rm PM}}{r^{k}}=$
$\displaystyle-\frac{Gm}{2a\sqrt{\pi}x_{s}}\left[\frac{r_{k}}{r}I^{\prime}_{3}(r)+\frac{r_{k}}{r}\sum_{i}\left(-I^{\prime}_{5}(r)+\frac{r_{i}^{2}}{2a^{2}x_{s}^{2}}I_{7}^{\prime}(r)\right)\Delta_{i}+\frac{r_{k}}{a^{2}x_{s}^{2}}I_{7}(r)\Delta_{k}\right.$
$\displaystyle+\frac{1}{2}\frac{r_{k}}{r}\sum_{i}\left(-I_{5}^{\prime}(r)+\left(3+\frac{r_{i}^{2}}{2a^{2}x_{s}^{2}}\right)I_{7}^{\prime}(r)-\frac{3r_{i}^{2}}{a^{2}x_{s}^{2}}I^{\prime}_{9}(r)+\frac{r_{i}^{4}}{4a^{4}x_{s}^{4}}I^{\prime}_{11}(r)\right)\Delta_{i}^{2}$
$\displaystyle\hskip
22.76228pt+\left(\frac{r_{k}}{2a^{2}x_{s}^{2}}I_{7}(r)-\frac{3r_{k}}{a^{2}x_{s}^{2}}I_{9}(r)+\frac{r_{k}^{3}}{2a^{4}x_{s}^{4}}I_{11}(r)\right)\Delta_{k}^{2}$
$\displaystyle+\frac{1}{2}\frac{r_{k}}{r}\sum_{i\neq
j}\left(I^{\prime}_{7}(r)-\frac{r_{i}^{2}+r_{j}^{2}}{2a^{2}x_{s}^{2}}I^{\prime}_{9}(r)+\frac{r_{i}^{2}r_{j}^{2}}{4a^{4}x_{s}^{4}}I^{\prime}_{11}(r)\right)\Delta_{i}\Delta_{j}$
$\displaystyle\left.\hskip 22.76228pt+\Delta_{k}\sum_{i\neq
k}\left(-\frac{r_{k}}{2a^{2}x_{s}^{2}}I_{9}(r)+\frac{r^{2}_{i}r_{k}}{4a^{4}x_{s}^{4}}I_{11}(r)\right)\Delta_{i}\right],$
(C.6)

where $I^{\prime}_{k}(r)$ is the derivative of $I_{k}(r)$ given by

$I_{k}^{\prime}(r)=\partialderivative{I_{k}(r)}{r}=\frac{2}{r}\biggl{[}\exp\biggl(-\frac{r^{2}}{4a^{2}x_{s}^{2}}\biggr{missing})-\frac{k-2}{2}\biggl{(}\frac{r}{2ax_{s}}\biggr{)}^{2-k}\gamma\biggl{(}\frac{k-2}{2},\frac{r^{2}}{4a^{2}x_{s}^{2}}\biggr{)}\biggr{]}.$
(C.7)

##  Appendix D Comparison between our simulation and the conventional separate
universe simulation

In this appendix, we show the result of the convergence test for the isotropic
background by comparing our simulations with the conventional separate
universe simulations.

###  D.1 Recap of the usual separate universe simulations

One way to incorporate the isotropic super-box mode is to change the
background parameters according to its value. This technique is based on the
fact that the flat FLRW universe with the spherically homogenious density
perturbation $\Delta_{0}$ is equivalent to the curved FLRW universe without
$\Delta_{0}$. This means $\Delta_{0}$ can be absorbed into the background
parameters in cosmological simulations. The relation of cosmological
parameters between the global and local universe can be characterized by the
ratio of Hubble parameters,

$\delta
h\equiv\frac{h_{W}}{h}=\sqrt{1-\frac{3}{5}\frac{{\Omega_{\mathrm{m}}}}{D(t_{0})}\Delta_{0}^{(1)}(t_{0})},$
(D.1)

where $h_{W}$ is the local Hubble parameter and hereafter the subscript W
denotes local quantities. In terms of $\delta h$, other cosmological paramters
are given by

$\displaystyle\Omega_{{\rm m}W}$ $\displaystyle={\Omega_{\mathrm{m}}}\delta
h^{-2},$ (D.2) $\displaystyle\Omega_{\Lambda W}$
$\displaystyle=\Omega_{\Lambda}\delta h^{-2},$ (D.3)
$\displaystyle\Omega_{KW}$ $\displaystyle=1-\delta h^{-2}.$ (D.4)

In general, the time and comoving coordinates are also different among the
global and local universes. Each cosmology has its own expansion history and
thus $a_{W}\neq 1$ when $a=1$. Therefore we need to find the relation between
the global and local scale factors at the same physical time $t$. We can
compute the difference between $a_{W}$ and $a$ by numerically solving Eq.
(2.12) with $\tau_{i}=0$, since the difference of the Friedmann equations in
the two cosmology solves the spherical collapse. We have to be careful to this
mapping of time in generating the initial conditions and determining the
output time.

As for the comoving length, it is common to use the unit of $\mathrm{Mpc}/h$
so the simulation box are given by $L~{}{\rm Mpc}/h$ and $L_{W}~{}{\rm
Mpc}/h_{W}$ in each cosmology. The choice of $L_{W}$ depends on what one wants
to measure directly from separate universe simulations. If one needs to obtain
the Eulerian response directly, $L_{W}$ is set to follow $a_{W}(t_{\rm
out})L_{W}/h_{W}=a(t_{\rm out})L/h$ with $t_{\rm out}$ being the output
physical tim. In order to get the Lagrangian response directly, one have to
set $L_{W}/h_{W}=L/h$. The former and the latter is called as the total
derivative method and the growth-dilation method respectively in Ref. [1]. We
employ the growth-dilation method where we set $L_{W}/h_{W}=L/h$ at all times
to share the ramdom fluctuations in the comoving scale in Mpc.

In this comparison study, we ran 6 pairs of separate universe simulations with
$\Delta_{0}^{(1)}=\pm 0.09$. The boxsize is $L=250~{}{\rm Mpc}/h$ and the
number of particles is $N_{\rm p}=1024^{3}$, which are the same as our high-
resolution simulations.

Figure 8:  Upper panels: Power spectrum responses to the large-scale
overdensity ($\Delta_{\rm 0}$) from our tidal separate universe simulations
(TSU, blue) and usual separate universe simulations (SU, orange)

at $z=0$ and 3. The horizontal dashed line presents the tree-level prediction
of the perturbation theory (Eq. (4.7)). Lower panels: the ratio of the two:
TSU/SU.

###  D.2 Power spectrum response

This subsection presents the comparison of the power spectrum response to the
density perturbation $\Delta_{0}$ from both our $N$-body simulation and 2LPT
with that from the usual separate universe simulations. We estimate
$R_{0}^{\rm L}(p)$ as

$\displaystyle R_{0}^{\rm
L}(p,z)=\frac{P_{W}(p,z;\Delta_{0}^{(1)}=+\epsilon)-P_{W}(p,z;\Delta_{0}^{(1)}=-\epsilon)}{2\epsilon
D(z)P(p)}.$ (D.5)

####  D.2.1 Convergence of $N-$body results

Fig. 8 shows the $R_{0}^{\rm L}(p)$ responses from our simulations and the
usual separate universe simulations. For both $z=3$ and $z=0$, our results are
in good agreement with the usual separate universe one, down to $k\simeq
10~{}h{\rm/Mpc}$.

####  D.2.2 On the valid scale of our 2LPT at high redshifts

Here we discuss the valid scales of our 2LPT by comparing it with the separate
univserse $N$-body results. Fig. 9 presents the $R_{0}^{\rm L}(p)$ responses
at $z=15,~{}10,$ and $7$ from our 2LPT and the usual separate universe
simulations. Note that the boxsize of our 2LPT is $L=100~{}{\rm Mpc}/h$ with
$N_{\rm p}=1024^{3}$ while the usual separate universe simulations have
$L=250~{}{\rm Mpc}/h$ with $N_{\rm p}=1024^{3}$. Setting the criteria to be
$5\%$ difference between our 2LPT and the separate universe simulations, we
conclude the responses from our 2LPT are reliable up to $k=~{}9,~{}4,~{}{\rm
and}~{}2h/{\rm Mpc}$ at $z=15,~{}10,$ and $7$, respectively.

Figure 9:  Upper panels: Power spectrum responses to the large-scale
overdensity ($\Delta_{\rm 0}$) from our modified 2LPT (blue) and usual
separate universe (SU) simulations (orange) at $z=15$, 10, and 7. Lower
panels: the ratio of the two: 2LPT/SU.

###  D.3 Linear bias

Finally we also compare results on the linear bias measured from our
simulations and the separate universe simulations. Fist, we directly measure
the Lagrangian linear bias $b_{1}^{\rm L}$ as

$\displaystyle b_{1}^{\rm L}(M,z)=\frac{N_{\rm
h}(M,z;\Delta_{0}^{\rm(1)}=+\epsilon)-N_{\rm
h}(M,z;\Delta_{0}^{\rm(1)}=-\epsilon)}{2\epsilon D(z)N_{\rm
h}(M,z;\Delta_{0}^{\rm(1)}=0)},$ (D.6)

where $N_{\rm h}(M,z)$ is the total number of halos at mass $M$ in the local
comoving volume. Then the Eulerian linear bias is computed as $b_{1}^{\rm
E}=b_{1}^{\rm L}+1$. Fig. 10 shows $b_{1}^{\rm E}$ at $z=0$ and $2$. For all
mass range our results agree with those from the usual separate universe
simulations. Together with the results about $R_{0}^{\rm L}(p)$, this suggests
that our implementation is correctly working.

Figure 10:  Upper panels: $b_{1}^{\rm E}$ from our tidal separate universe
simulations (TSU, blue) and usual separate universe simulations (SU, red).
Lower panels: the ratio of the two: TSU/SU.

##  Appendix E Results from the inertial tensor $I_{ij}$

In this appendix, we summarize the shape response results when using the
inertial tensor $I_{ij}$. Fig. 11 shows the boxsize or equivalently resolution
dependence of $b_{K}$. Unlike using $J_{ij}$, there is no difference between
$b_{K}$ at 6th and 7th mass-bin from $L=250~{}{\rm Mpc}/h$ and $L=1~{}{\rm
Gpc}/h$. This can be explained by the enough number of particle to determine
the halo shape when using $J_{ij}$ in $1~{}{\rm Gpc}/h$ simulations since the
number of particles used to define the halo shape is effectively higher in
$I_{ij}$ than in $J_{ij}$. Given these results, in the following in this
appendix we use $250~{}{\rm Mpc}/h$, $1~{}{\rm Gpc}/h$, and $3~{}{\rm Gpc}/h$
simulations for 1st-6th, 7th-9th, and 10th-11th mass-bin, respectively.

Fig. 11 also compares $b_{K}$ from different kinds of tides and Fig. 12 shows
the time evolution of $b_{K}$ when using $I_{ij}$. While the amplitude of
$b_{K}$ from $I_{ij}$ differs from $J_{ij}$, $I_{ij}$ results show the same
trend as $J_{ij}$. The universal behavior between $b_{1}^{\rm E}$ and $b_{K}$
is also found from $I_{ij}$ results as shown in Fig. 13. This implies that
this relation is indeed “universal”, regardless of the definition of shapes.
In $I_{ij}$ case we can fit this relation as

$\displaystyle
b_{K}=\frac{0.09302-0.1289b_{1}^{\mathrm{E}}}{1+0.3541b_{1}^{\mathrm{E}}}$
(E.1)

In Fig. 14 we provide the comparison of $b_{K}$ measured from $I_{ij}$ and
$J_{ij}$ for various redshifts and mass-bins. For all points, the amplitude of
$|b_{K}|$ is greater when using $I_{ij}$ than when using $J_{ij}$ and its
ratio does not change significantly over all the redshift and mass range. Thus
the choice of the definition of shapes does not change the dependence of
$b_{K}$ on mass or redshift as already seen in Fig. 12.

Finally, in Fig. 15 and Fig. 16, we present that the secondary dependence of
$b_{K}$ on halo concentration and the axis-ratio is also found when using
$I_{ij}$ with the same trend. This suggests that the secondary dependence of
$b_{K}$ is genuine.

Figure 11: Linear alignment coefficient, $b_{K}$, for the reduced inertial
tensor, $I_{ij}$, at $z=0$. The left plot shows measurements from simulations
of different box sizes, and the right one shows measurements from different
tidal types.  Figure 12: Linear alignment coefficient, $b_{K}$, for the
reduced inertial tensor, $I_{ij}$, at various redshifts. Here we combine the
results from difference box sizes and different tides.  Figure 13: Linear
alignment coefficient, $b_{K}$, from $I_{ij}$, as a function of the Eulerian
linear bias, $b_{1}^{\mathrm{E}}$, combining all redshift and mass
information. The dashed curve is the fitting given in Eq. (E.1).

Figure 14:  Upper panels: The shape definition dependence of the linear
alignment coefficient, $b_{K}$, at various redshifts and masses: $b_{K}$ from
$I_{ij}$ (blue) and $J_{ij}$ (red). Lower panels: The ratio of the two:
$I_{ij}/J_{ij}$.

Figure 15:  Upper panels: The halo concentration dependence of the linear
alignment coefficient, $b_{K}$, at various redshifts and masses: $b_{K}$ from
high concentration (red) and low concentration (blue). Lower panels: The
median concentration on which we divided halo samples.

Figure 16:  Upper panels: The axis-ratio ($q=I_{3}/I_{1}$ with $I_{1}\geq
I_{3}$) dependence of the linear alignment coefficient, $b_{K}$, at various
redshifts and masses: $b_{K}$ from high $q$ (low ellipticity, blue) and low
$q$ (high ellipticity, red). Lower panels: The median axis-ratio by which we
divided halo samples.

## References

  * [1] Y. Li, W. Hu and M. Takada, _Super-Sample Covariance in Simulations_ , _Phys. Rev. D_ 89 (2014) 083519 [1401.0385].
  * [2] C. Wagner, F. Schmidt, C.-T. Chiang and E. Komatsu, _Separate universe simulations_ , _Monthly Notices of the Royal Astronomical Society: Letters_ 448 (2015) L11.
  * [3] T. Baldauf, U. Seljak, L. Senatore and M. Zaldarriaga, _Linear response to long wavelength fluctuations using curvature simulations_ , _JCAP_ 09 (2016) 007 [1511.01465].
  * [4] Y. Li, W. Hu and M. Takada, _Separate Universe Consistency Relation and Calibration of Halo Bias_ , _Phys. Rev. D_ 93 (2016) 063507 [1511.01454].
  * [5] T. Lazeyras, C. Wagner, T. Baldauf and F. Schmidt, _Precision measurement of the local bias of dark matter halos_ , _JCAP_ 02 (2016) 018 [1511.01096].
  * [6] J. Bond and S. Myers, _The peak-patch picture of cosmic catalogs. i. algorithms_ , _The Astrophysical Journal Supplement Series_ 103 (1996) 1.
  * [7] K. Akitsu, M. Takada and Y. Li, _Large-scale tidal effect on redshift-space power spectrum in a finite-volume survey_ , _Phys. Rev. D_ 95 (2017) 083522 [1611.04723].
  * [8] A.S. Schmidt, S.D. White, F. Schmidt and J. Stücker, _Cosmological N-Body Simulations with a Large-Scale Tidal Field_ , _Mon. Not. Roy. Astron. Soc._ 479 (2018) 162 [1803.03274].
  * [9] J. Stücker, A. Schmidt, S.D. White, F. Schmidt and O. Hahn, _Measuring the Tidal Response of Structure Formation: Anisotropic Separate Universe Simulations using TreePM_ , 2003.06427.
  * [10] S. Masaki, T. Nishimichi and M. Takada, _Anisotropic separate universe simulations_ , _Mon. Not. Roy. Astron. Soc._ 496 (2020) 483 [2003.10052].
  * [11] S. Masaki, T. Nishimichi and M. Takada, _Impacts of pre-initial conditions on anisotropic separate universe simulations: a boosted tidal response in the epoch of reionization_ , 2007.08727.
  * [12] K. Akitsu and M. Takada, _Impact of large-scale tides on cosmological distortions via redshift-space power spectrum_ , _Phys. Rev. D_ 97 (2018) 063527 [1711.00012].
  * [13] Y. Li, M. Schmittfull and U. Seljak, _Galaxy power-spectrum responses and redshift-space super-sample effect_ , _Journal of Cosmology and Astroparticle Physics_ 2018 (2018) 022.
  * [14] K. Akitsu, N.S. Sugiyama and M. Shiraishi, _Super-sample tidal modes on the celestial sphere_ , _Phys. Rev. D_ 100 (2019) 103515 [1907.10591].
  * [15] A. Barreira, E. Krause and F. Schmidt, _Complete super-sample lensing covariance in the response approach_ , _JCAP_ 06 (2018) 015 [1711.07467].
  * [16] P. Catelan, M. Kamionkowski and R.D. Blandford, _Intrinsic and extrinsic galaxy alignment_ , _Mon. Not. Roy. Astron. Soc._ 320 (2001) L7 [astro-ph/0005470].
  * [17] C.M. Hirata and U. Seljak, _Intrinsic alignment-lensing interference as a contaminant of cosmic shear_ , _Phys. Rev. D_ 70 (2004) 063526 [astro-ph/0406275].
  * [18] R. Mandelbaum, C.M. Hirata, M. Ishak, U. Seljak and J. Brinkmann, _Detection of large scale intrinsic ellipticity-density correlation from the sloan digital sky survey and implications for weak lensing surveys_ , _Mon. Not. Roy. Astron. Soc._ 367 (2006) 611 [astro-ph/0509026].
  * [19] A. Taruya and T. Okumura, _Improving geometric and dynamical constraints on cosmology with intrinsic alignments of galaxies_ , 2001.05962.
  * [20] T. Kurita, M. Takada, T. Nishimichi, R. Takahashi, K. Osato and Y. Kobayashi, _Power spectrum of halo intrinsic alignments in simulations_ , 2004.12579.
  * [21] F. Schmidt and D. Jeong, _Large-Scale Structure with Gravitational Waves II: Shear_ , _Phys. Rev. D_ 86 (2012) 083513 [1205.1514].
  * [22] F. Schmidt, E. Pajer and M. Zaldarriaga, _Large-Scale Structure and Gravitational Waves III: Tidal Effects_ , _Phys. Rev. D_ 89 (2014) 083507 [1312.5616].
  * [23] F. Schmidt, N.E. Chisari and C. Dvorkin, _Imprint of inflation on galaxy shape correlations_ , _JCAP_ 10 (2015) 032 [1506.02671].
  * [24] K. Akitsu, T. Kurita, T. Nishimichi, M. Takada and S. Tanaka, _Imprint of anisotropic primordial non-Gaussianity on halo intrinsic alignments in simulations_ , 2007.03670.
  * [25] K. Kogai, K. Akitsu, F. Schmidt and Y. Urakawa, _Galaxy imaging surveys as spin-sensitive detector for cosmological colliders_ , 2009.05517.
  * [26] N.Y. Gnedin, A.V. Kravtsov and D.H. Rudd, _Implementing the DC Mode in Cosmological Simulations with Supercomoving Variables_ , _Astrophys. J. Suppl._ 194 (2011) 46 [1104.1428].
  * [27] B.D. Sherwin and M. Zaldarriaga, _Shift of the baryon acoustic oscillation scale: A simple physical picture_ , _Physical Review D_ 85 (2012) 103523.
  * [28] Y.B. Zel’Dovich, _Gravitational instability: An approximate theory for large density perturbations._ , _Astronomy and astrophysics_ 5 (1970) 84.
  * [29] M. Crocce, S. Pueblas and R. Scoccimarro, _Transients from initial conditions in cosmological simulations_ , _MNRAS_ 373 (2006) 369 [astro-ph/0606505].
  * [30] J.S. Bagla, _TreePM: A Code for Cosmological N-Body Simulations_ , _Journal of Astrophysics and Astronomy_ 23 (2002) 185 [astro-ph/9911025].
  * [31] J.S. Bagla and S. Ray, _Performance characteristics of TreePM codes_ , _New A_ 8 (2003) 665 [astro-ph/0212129].
  * [32] T. Quinn, N. Katz, J. Stadel and G. Lake, _Time stepping N-body simulations_ , _ArXiv Astrophysics e-prints_ (1997) [astro-ph/9710043].
  * [33] V. Springel, S.D. White, A. Jenkins, C.S. Frenk, N. Yoshida, L. Gao et al., _Simulations of the formation, evolution and clustering of galaxies and quasars_ , _nature_ 435 (2005) 629.
  * [34] Planck Collaboration, P.A.R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont et al., _Planck 2015 results. XIII. Cosmological parameters_ , _A &A_ 594 (2016) A13 [1502.01589].
  * [35] D. Blas, J. Lesgourgues and T. Tram, _The cosmic linear anisotropy solving system (class). part ii: approximation schemes_ , _Journal of Cosmology and Astroparticle Physics_ 2011 (2011) 034.
  * [36] P. Mansfield and C. Avestruz, _How Biased Are Halo Properties in Cosmological Simulations?_ , _Mon. Not. Roy. Astron. Soc._ 500 (2020) 3309 [2008.08591].
  * [37] S.R. Knollmann and A. Knebe, _Ahf: Amiga’s halo finder_ , _The Astrophysical Journal Supplement Series_ 182 (2009) 608.
  * [38] T. Okumura, Y.P. Jing and C. Li, _Intrinsic Ellipticity Correlation of SDSS Luminous Red Galaxies and Misalignment with Their Host Dark Matter Halos_ , _ApJ_ 694 (2009) 214 [0809.3790].
  * [39] A. Faltenbacher, C. Li, S.D.M. White, Y.-P. Jing, Shu-DeMao and J. Wang, _Alignment between galaxies and large-scale structure_ , _Research in Astronomy and Astrophysics_ 9 (2009) 41 [0811.1995].
  * [40] T. Okumura and Y.P. Jing, _The Gravitational Shear-Intrinsic Ellipticity Correlation Functions of Luminous Red Galaxies in Observation and in the $\Lambda$CDM Model_, _ApJ_ 694 (2009) L83 [0812.2935].
  * [41] M. Zemp, O.Y. Gnedin, N.Y. Gnedin and A.V. Kravtsov, _On determining the shape of matter distributions_ , _Astrophys. J. Suppl._ 197 (2011) 30 [1107.5582].
  * [42] J. Shi, T. Kurita, M. Takada, K. Osato, Y. Kobayashi and T. Nishimichi, _Power Spectrum of Intrinsic Alignments of Galaxies in IllustrisTNG_ , _arXiv e-prints_ (2020) arXiv:2009.00276 [2009.00276].
  * [43] D.M. Schmitz, C.M. Hirata, J. Blazek and E. Krause, _Time evolution of intrinsic alignments of galaxies_ , _JCAP_ 07 (2018) 030 [1805.02649].
  * [44] F. Prada, A.A. Klypin, A.J. Cuesta, J.E. Betancort-Rijo and J. Primack, _Halo concentrations in the standard $\lambda$ cold dark matter cosmology_, _Monthly Notices of the Royal Astronomical Society_ 423 (2012) 3018.
  * [45] N. Dalal, M. White, J. Bond and A. Shirokov, _Halo Assembly Bias in Hierarchical Structure Formation_ , _Astrophys. J._ 687 (2008) 12 [0803.3453].
  * [46] A. Obuljen, N. Dalal and W.J. Percival, _Anisotropic halo assembly bias and redshift-space distortions_ , _JCAP_ 10 (2019) 020 [1906.11823].
  * [47] A. Obuljen, W.J. Percival and N. Dalal, _Detection of anisotropic galaxy assembly bias in BOSS DR12_ , _JCAP_ 10 (2020) 058 [2004.07240].