Also at ]Key Laboratory of Artificial Structures and Quantum Control (Ministry
of Education), School of Physics and Astronomy, Shanghai Jiao Tong University,
800 Dong Chuan Road, Shanghai 200240, China Also at ]Key Laboratory of
Artificial Structures and Quantum Control (Ministry of Education), School of
Physics and Astronomy, Shanghai Jiao Tong University, 800 Dong Chuan Road,
Shanghai 200240, China

# Proposal for measuring Newtonian constant of gravitation at an exceptional
point in an optomechanical system

Lei Chen  [<EMAIL_ADDRESS>Jian Liu  [  Ka-di Zhu  Key
Laboratory of Artificial Structures and Quantum Control (Ministry of
Education), School of Physics and Astronomy, Shanghai Jiao Tong University,
800 Dong Chuan Road, Shanghai 200240, China

###### Abstract

We develop a quantum mechanical method of measuring the Newtonian constant of
gravitation, $G$. In this method, an optomechanical system consisting of two
cavities and two membrane resonators is used. The added source mass would
induce the shifts of the eigenfrequencies of the supermodes. Via detecting the
shifts, we can perform our measurement of $G$. Furthermore, our system can
features exceptional point (EP) which are branch point singularities of the
spectrum and eigenfunctions. In the paper, we demonstrate that operating the
system at EP can enhance our measurement of $G$. In addition, we derive the
relationship between EP enlarged eigenfrequency shift and the Newtonian
constant. This work provides a way to engineer EP-assisted optomechanical
devices for applications in the field of precision measurement of $G$

††preprint: APS/123-QED

##  I Introduction

Newton’s law of gravitation is often written as

$F=G\frac{m_{1}m_{2}}{r^{2}},$ (1)

where $m_{1}$ and $m_{2}$ are the masses of two particles, $r$ is the distance
between them and $G$ is the gravitational constant. Though the absolute value
of $G$ has been measured by many experiments [1, 2, 3], two big puzzles still
exist nowadays. One is the value of G ramains the least preciously known of
the fundamental constants [4]. The Committee on Data for Science and
Technology (CODATA) recommended value of $G$ based on the 2014 adjustment is
$6.67408(31)\times{10}^{-11}m^{3}{kg}^{-1}s^{-2}$, and the relative
measurement uncertainty is as high as $4.7\times{10}^{-5}$ [5]. The other
puzzle is the experimental values of $G$ are not consistent with each
other[6]. To resolve these puzzles strongly motivates us to develop different
measurement methods [7].

In this paper, we develop a new method for determining $G$. In our method, a
quantum system consisting of two cavities and two membranes is considered.
Gravitational force gradient originating from the source mass causes shifts of
the resonant frequencies of two membranes, resulting the shifts of the
eigenfrequencies of the supermodes emerging in our system. Based on this, we
establish our measurement principles. Exceptional points are branch point
singularities of the spectrum and eigenfunctions [8] which occur generically
in eigenvalue problems that depend on a parameter[9]. Because EPs in open
quantum and wave systems can exhibit a strong spectral response to
perturbations[10], their potential applications in precision measurement have
been investigated both theoretically[11, 12] and experimentally[13, 14]. In
this paper, EP effect is investigated numerically. The relating results
indicate that this effect can enhance the detection of the eigenfrequency
shift. Furthermore, the relationship between eigenfrequency shift and
Newtonian constant is established provided that our system is operated at EP.
In summary, we propose an EP-based quantum mechanical method for measuring
$G$. Finally we expect our work can enrich the experimental approaches of
determining $G$.

The remainder of the paper is organized as follows: In Sec. II we demonstrate
the theoretical framework, in Sec. III we present the measurement principles,
in Sec. IV we summarize the paper and provide an outlook.

Figure 1:  Schematic setup. Two Cavities, whose end-mirrors are two membrane
resonators, are driven with two detuned lasers.

##  II Theoretical framework

We consider a system composed of two optomechanical cavities the two end-
mirrors of which are two nanosized membrane resonators [15]. Cavity 1(2) is
driven with a red (blue) detuned laser. By symmetrically driving the cavities
we can engineer either mechanical gain or mechanical loss [16]. Note that at
the gain and loss balance the system features EP. In the rotating frame of the
driving fields, the Hamiltonian of the system could be written as ($\hbar=1$)

$H=H_{1}+H_{2}+H_{3},$ (2)

where

$\displaystyle H_{1}$
$\displaystyle=\sum_{j=1,2}\omega_{j}b_{j}^{+}b_{j}-\Delta_{j}a_{j}^{+}a_{j}-ga_{j}^{+}a_{j}(b_{j}^{+}+b_{j})$
$\displaystyle H_{2}$
$\displaystyle=-\varepsilon(b_{1}b_{2}^{+}+b_{1}^{+}b_{2})$ $\displaystyle
H_{3}$ $\displaystyle=\sum_{j=1,2}E(a_{j}^{+}+a_{j}).$ (3)

Here $a_{j}$ and $b_{j}$ are the annihilation bosonic field operators
describing the optical and mechanical resonators. $\omega_{j}$ is the
mechanical frequency of the $j^{th}$ resonator and
$\Delta_{j}=\omega_{P}^{j}-\omega_{c}^{j}$ is the optical detuning between the
optical ($\omega_{P}^{j}$) and the cavity ($\omega_{c}^{j}$) frequencies.
$\varepsilon$ is the coupling strength between the two mechanical resonators,
$g$ is the vacuum optomechanical coupling strength. The amplitude of the
driving pump is $E$.

The quantum Langevin equations (QLEs) for the operations of the optical and
the mechanical modes are derived from Eq. (3) as

$\displaystyle{\dot{a}}_{j}$
$\displaystyle=[i(\Delta_{j}+g(b_{j}^{+}+b_{j}))-\frac{\kappa}{2}]a_{j}-i\sqrt{\kappa}(\alpha^{in}+\xi_{a_{j}})$
$\displaystyle{\dot{b}}_{j}$
$\displaystyle=-(i\omega_{j}+\frac{\gamma_{m}}{2})b_{j}+i\varepsilon
b_{k}+iga_{j}^{+}a_{j}+\sqrt{\gamma_{m}}\xi_{b_{j}},$ (4)

where optical ($\kappa$) and mechanical ($\gamma_{m}$) dissipations have been
added, $k$ is defined as $k=3-j$, and the amplitude of the driving pump has
been substituted as $E=\sqrt{\gamma_{m}}\alpha^{in}$ in order to account for
losses. The term $\xi_{a_{j}}$ ($\xi_{b_{j}}$) denotes the optical (thermal)
Langevin noise at room temperature. We seek to investigate in the classical
limit, where photon and phonon numbers are assumed large in the system, and
noise terms can be neglected in our analysis. Then we can derive the following
set of nonlinear equations:

$\displaystyle{\dot{\alpha}}_{j}$
$\displaystyle=[i(\Delta_{j}+g(\beta_{j}^{*}+\beta_{j}))-\frac{\kappa}{2}]\alpha_{j}-i\sqrt{\kappa}\alpha^{in}$
$\displaystyle{\dot{\beta}}_{j}$
$\displaystyle=-(i\omega_{j}+\frac{\gamma_{m}}{2})\beta_{j}+i\varepsilon\beta_{k}+ig\alpha_{j}^{*}\alpha_{j}.$
(5)

Here $\alpha_{j}=\langle a_{j}\rangle$ and $\beta_{j}=\langle b_{j}\rangle$.

To identify the EP feature, we approach the limit cycle oscillations by the
ansatz, $\beta_{j}(t)={\bar{\beta}}_{j}+A_{j}exp(-i\omega_{lock}t)$.
${\bar{\beta}}_{j}$ is a constant shift in the origin of the movement, $A_{j}$
is taken to be a slowly varying function of time, and $\omega_{lock}$ is the
mechanical locked frequency. According to [17], we can derive

$\displaystyle i\frac{\partial\Psi}{\partial t}=H_{eff}\Psi,$ (6)

where $\Psi=(\beta_{1},\beta_{2})^{T}$ is the state vector and the effective
Hamiltonian is

$H_{eff}=\begin{bmatrix}\omega^{1}_{eff}-i{\frac{\gamma^{1}_{eff}}{2}}&-\varepsilon\\\
-\varepsilon&\omega^{2}_{eff}-i{\frac{\gamma^{2}_{eff}}{2}}\end{bmatrix}$

. Here $\omega^{j}_{eff}=\omega_{j}+\Delta\omega_{j}$ and
$\gamma^{j}_{eff}=\gamma_{m}+\gamma_{opt}^{j}$ define the effective
frequencies and the effective damping respectively. The optical spring effect
$\Delta\omega_{j}$ and the optical damping $\gamma_{opt}^{j}$ are given by

$\displaystyle\Delta\omega_{j}=-\frac{2\kappa(g\alpha^{in})^{2}}{\omega_{lock}\zeta_{j}}Re(\sum_{n}\frac{J_{n+1}(-\zeta_{j})J_{n}(-\zeta_{j})}{h_{n+1}^{j*}{h_{n}^{j}}})$
(7)

and

$\displaystyle\gamma_{opt}^{j}=\frac{2(\kappa
g\alpha^{in})^{2}}{\zeta_{j}}\sum_{n}\frac{J_{n+1}(-\zeta_{j})J_{n}(-\zeta_{j})}{{|h_{n+1}^{j*}{h_{n}^{j}}|}^{2}}.$
(8)

$\zeta_{j}=\frac{2gRe(A_{j})}{\omega_{lock}}$,${\tilde{\Delta}}_{j}=\Delta_{j}+2gRe({\bar{\beta}}_{j})$,
$h_{n}^{j}=i(n\omega_{lock}-{\tilde{\Delta}}_{j})+\frac{\kappa}{2}$, and
$J_{n}$ is the Bessel function of the first kind of order $n$. The eigenvalues
of the effective Hamiltonian can be described as

$\displaystyle\tau_{\pm}=\frac{\omega^{1}_{eff}+\omega^{2}_{eff}}{2}-i\frac{\gamma^{1}_{eff}+\gamma^{2}_{eff}}{4}\pm\frac{\sqrt{\varDelta}}{2},$
(9)

with

$\displaystyle\varDelta=[(\omega^{1}_{eff}-\omega^{2}_{eff})+\frac{i}{2}(\gamma^{2}_{eff}-\gamma^{1}_{eff})]^{2}+4\varepsilon^{2}.$
(10)

The EP of our system appears if $\Delta=0$. The eigenfrequencies and the
dissipations of the supermodes are defined as the real and imaginary parts of
the eigenvalues, i.e.,$\nu_{\pm}=Re(\tau_{\pm})$, and
$\Upsilon_{\pm}=Im(\tau_{\pm})$. At the EP, both pairs of eigenfrequencies and
dissipations coalesce, i.e., $\nu_{+}=\nu_{-}$ and
$\Upsilon_{+}=\Upsilon_{-}$.

For simplicity, we assume the two resonators are degenerated, causing
$\omega_{1}=\omega_{2}=\omega_{r}$. We adopt a weak driving regime, resulting
$\Delta\omega_{j}\ll\omega_{r}$. We further assume
$\gamma_{opt}^{j}=(\alpha^{in})^{2}\eta_{j}$ with $\eta_{1}\neq\eta_{2}$. Then
the eigenvalues of the Hamiltonian can be written as

$\displaystyle\tau_{\pm}\approx$
$\displaystyle\omega_{r}-i\frac{2\gamma_{m}+(\eta_{1}+\eta_{2})(\alpha^{in})^{2}}{4}$
$\displaystyle\pm\sqrt{\varepsilon^{2}-\frac{(\alpha^{in})^{4}(\eta_{2}-\eta_{1})^{2}}{16}}.$
(11)

As a result, the EP condition is tranformed into
$\varepsilon^{2}-\frac{(\alpha^{in})^{4}(\eta_{2}-\eta_{1})^{2}}{16}=0$, from
which we derive that
$\alpha^{in}=\alpha_{EP}^{in}=\sqrt{\frac{4\varepsilon}{|\eta_{2}-\eta_{1}|}}$.
Now we use an example to illustrate the EP feature (see Fig. 2). In this
example, the parameters used are
$\eta_{1}={10}^{-6},\eta_{1}=2\times{10}^{-6},\varepsilon={10}^{-2}\omega_{r},\gamma_{m}={10}^{-4}\omega_{r}.$
The corresponding EP is $\alpha_{EP}^{in}=200\omega_{r}^{1/2}$. From the
figure, it is seen that two eigenfrequencies coalesce at the EP as well as two
dissipations. Based on the above analysis about our system, we present the
principles of measurement in the next section.

Figure 2:  Illustration of EP feature. Eigenfrequencies $\nu_{\pm}$ and
dissipations $\Upsilon_{\pm}$ versus the driving strength $\alpha^{in}$. the
parameters used are
$\eta_{1}={10}^{-6},\eta_{1}=2\times{10}^{-6},\varepsilon={10}^{-2}\omega_{r},\gamma_{m}={10}^{-4}\omega_{r}.$

##  III Principles of measurement

Figure 3:  (a)The locations of two membranes and the sphere . (b) A cross
section of the diagram in (a).

In our proposed method, a sphere with uniformly distributed mass $m$, radius
$R$ and density $\rho$ is utilized. This sphere and two membranes of our
system are placed as shown in Fig. 3(a) Note that the sizes of two membrane
resonators are much smaller than $R$. The distances between two membranes and
the center of the sphere are $a_{1}$ and $a_{2}$ respectively, which is shown
in Fig. 3(b). We assume that exotic forces such as non-Newtonian gravitylike
forces applied to the two membranes can be neglected. Because of the
gravitational forces generated by the source mass, the resonant frequencies of
two resonators are modified by $\delta\omega_{1}$ and $\delta\omega_{2}$
respectively. Here $\delta\omega_{1}$ and $\delta\omega_{2}$ are defined as
$\delta\omega_{j}=\omega_{j}^{\prime}-\omega_{j}$, where $\omega_{j}^{\prime}$
is the perturbed resonant frequency. From [18] we obtain

$\frac{\delta\omega_{j}}{\omega_{j}}\approx\frac{1}{2m_{j}\omega_{j}^{2}}\frac{\partial
F(a_{j})}{\partial a_{j}},$ (12)

where $m_{j}$ is the mass of $j$th membrane, and the gravitational force
acting on $j$th membrane is

$F(a_{j})=\frac{GMm_{j}}{a_{j}^{2}}.$ (13)

From the above two equations, we derive

$\delta\omega_{j}\approx-\frac{GM}{\omega_{j}a_{j}^{3}}.$ (14)

From section 2 we know $\omega_{1}=\omega_{2}=\omega_{r}$. Then we assume that
$a_{1}\gg a_{2}$, resulting $|\delta\omega_{1}|\ll|\delta\omega_{2}|$.

Figure 4:  Eigenfrequencies $\nu_{\pm}$ and eigenfrequencies shifts
$\Delta\nu_{\pm}$ versus the driving field $\alpha_{in}$ around EP. The black
curves are for $\nu_{+}$ and $\Delta\nu_{+}$, while red ones for $\nu_{-}$ and
$\Delta\nu_{-}$. (a)-(b)The parameters used are
$\eta_{1X}={10}^{-6},\eta_{2X}=2\times{10}^{-6},\varepsilon_{X}={10}^{-2}\omega_{r},\gamma_{mX}={10}^{-4}\omega_{r}.$
(c)-(d)The parameters used are
$\eta_{1Y}=3\times{10}^{-7},\eta_{2Y}=7\times{10}^{-7},\varepsilon_{Y}={10}^{-3}\omega_{r},\gamma_{mY}={10}^{-3}\omega_{r}$.
(e)-(f) The parameters used are
$\eta_{1Z}={10}^{-5},\eta_{2Z}=4\times{10}^{-5},\varepsilon_{Z}=3\times{10}^{-3}\omega_{r},\gamma_{mZ}=2\times{10}^{-3}\omega_{r}.$

As a result of resonant frequency modification of two membranes, the
eigenvalues of the Hamiltonian can be rewritten as

$\displaystyle\tau_{\pm}^{\prime}\approx$
$\displaystyle\omega_{r}+\frac{\delta\omega_{2}}{2}-i\frac{2\gamma_{m}+(\alpha^{in})^{2}(\eta_{1}+\eta_{2})}{4}$
$\displaystyle\pm\frac{\sqrt{[-\delta\omega_{2}+\frac{i}{2}(\alpha^{in})^{2}(\eta_{2}-\eta_{1})]^{2}+4\varepsilon^{2}}}{2}$
(15)

For convenience, we substitute $\delta\omega$ for $\delta\omega_{2}$ in the
remaining section. The shifts of the eigenfrequencies of the supermodes caused
by the source mass can be described as

$\Delta\nu_{\pm}=Re(\tau_{\pm}^{\prime})-Re(\tau_{\pm}).$ (16)

Then from Eqs. (11), (15)-(16) we derive that

$\displaystyle\Delta\nu_{\pm}=Re(\frac{\delta\omega}{2}$
$\displaystyle\pm\sqrt{\varepsilon^{2}-\frac{(\alpha^{in})^{4}(\eta_{2}-\eta_{1})^{2}}{16}+\frac{(\delta\omega)^{2}-i\delta\omega(\alpha^{in})^{2}(\eta_{2}-\eta_{1})}{4}}$
$\displaystyle\mp\sqrt{\varepsilon^{2}-\frac{(\alpha^{in})^{4}(\eta_{2}-\eta_{1})^{2}}{16}}$
(17)

In the following, we consider three cases denoted by X, Y, and Z where we
operate our system near the EP.

In the case of X, the parameters used are
$\eta_{1X}={10}^{-6},\eta_{2X}=2\times{10}^{-6},\varepsilon_{X}={10}^{-2}\omega_{r},\gamma_{mX}={10}^{-4}\omega_{r},$
and the corresponding EP is $\alpha_{EPX}^{in}=200\omega_{r}^{1/2}$. In Y,
$\eta_{1Y}=3\times{10}^{-7},\eta_{2Y}=7\times{10}^{-7},\varepsilon_{Y}={10}^{-3}\omega_{r},\gamma_{mY}={10}^{-3}\omega_{r},\alpha_{EPY}^{in}=100\omega_{r}^{1/2}.$
In Z,
$\eta_{1Z}={10}^{-5},\eta_{2Z}=4\times{10}^{-5},\varepsilon_{Z}=3\times{10}^{-3}\omega_{r},\gamma_{mZ}=2\times{10}^{-3}\omega_{r},\alpha_{EPZ}^{in}=20\omega_{r}^{1/2}.$
Figure 4 illustrates how eigenfrequencies ($\nu_{\pm}$) and eigenfrequencies
shifts ($\Delta\nu_{\pm}$) undergo the EP corresponding to different
perturbation. (a)-(b), (c)-(d) and (e)-(f) are for X, Y, Z respectively. The
black curves are for $\nu_{+}$ and $\Delta\nu_{+}$, while red ones for
$\nu_{-}$ and $\Delta\nu_{-}$. From this figure, we find that in all nine
situations (different perturbation and different three cases) we would obtain
a relative bigger positive frequency shift ($|\Delta\nu_{+}|$) and a maximum
negative eigenfrequency shift ($|\Delta\nu_{-}|$) at EPs.

Figure 5:  Nine ratios ($\Gamma$)correspond to three cases(X, Y, and Z) and
different perturbation ($|\delta\omega|$).

To further demonstrate the EP effect, we define ratio $\Gamma$ as
$\Gamma=\frac{|\Delta\nu_{-}|_{EP}}{|\Delta\nu_{-}|_{min}}$, where
$|\Delta\nu_{-}|_{EP}$ is the value of $|\Delta\nu_{-}|$ related to
$\alpha_{EP}^{in}$ and $|\Delta\nu_{-}|_{min}$ is the minimum value of
$|\Delta\nu_{-}|$. Nine values of $\Gamma$ correspond to three cases(X, Y, and
Z) and different perturbation ($|\delta\omega|$) are shown in Fig. 5. where
square, circles and triangle correspond to X, Y and Z respectively. We can
find that in any of three cases a smaller $|\delta\omega|$ corresponds to a
bigger $\Gamma$. Since $\delta\omega$ in practical experiments would be very
tiny, we conclude that operating the system at EP may make a large
contribution to our measurement.

We assume our system is operated at EP and investigate the relationship
between EP enlarged eigenfrequency shift and the Newtonian constant.
Accordingly there is
$\alpha^{in}=\alpha^{in}_{EP}=\sqrt{\frac{4\varepsilon}{|\eta_{2}-\eta_{1}|}}$.
With the assumption $\eta_{2}>\eta_{1}$, according to Eq. (17) we derive that

$\Delta\nu_{\pm}=Re(\frac{\delta\omega}{2}\pm\sqrt{\frac{(\delta\omega)^{2}-i\delta\omega
4\varepsilon}{4}}).$ (18)

Then, from Eq. (18) we derive

$\Delta\nu_{\pm}=\frac{\delta\omega}{2}(1\mp\sqrt{\frac{1+\sqrt{1+\frac{16\varepsilon^{2}}{(\delta\omega)^{2}}}}{2}}).$
(19)

Next we focus on $\Delta\nu_{-}$. From the above, we obtain

$\delta\omega\approx-\frac{GM}{\omega_{r}a_{2}^{3}}.$ (20)

From Eqs. (19)-(20), we derive

$\displaystyle\Delta\nu_{-}=-\frac{GM}{2\omega_{r}a_{2}^{3}}-\frac{GM}{2\omega_{r}a_{2}^{3}}\sqrt{\frac{1+\sqrt{1+\frac{16\varepsilon^{2}\omega_{r}^{2}a_{2}^{6}}{G^{2}M^{2}}}}{2}}.$
(21)

We assume $a_{2}-R\ll R$. Consequently there is $a_{2}\approx R$. Since
$M=\frac{4}{3}\pi R^{3}\rho$, equation (21) can be rewritten as

$\Delta\nu_{-}\approx-\frac{2\pi G\rho}{3\omega_{r}}-\frac{2\pi
G\rho}{3\omega_{r}}\sqrt{\frac{1+\sqrt{1+\frac{9\varepsilon^{2}\omega_{r}^{2}}{G^{2}\pi^{2}\rho^{2}}}}{2}}.$
(22)

Then we obtain

$|\Delta\nu_{-}|=\frac{2\pi G\rho}{3\omega_{r}}+\frac{2\pi
G\rho}{3\omega_{r}}\sqrt{\frac{1+\sqrt{1+\frac{9\varepsilon^{2}\omega_{r}^{2}}{G^{2}\pi^{2}\rho^{2}}}}{2}}.$
(23)

Till now, we have derived the relationship between eigenfrequency shift and
the Newtonian constant. Then we can determine the value of $G$ and the
relative uncertainty via the detection of $|\Delta\nu_{-}|$.

Figure 6:  (a) The adopted values of G in CODATA-2014 adjustment. (b)-(c)
$|\Delta\nu_{-}|$ as a function of $G$. The parameters used are
$\varepsilon={10}^{-2}\omega_{r},\omega_{r}=2\times{10}^{9}Hz$. (b) deals with
three spheres with different densities, which are made of stainless steel,
lead, and tungsten respectively. On the contrary, (c) only deal with a
Tungsten sphere. Note that in (c), $G$ takes values of
$6.67408(31)\times{10}^{-11}m^{3}\bullet{kg}^{-1}\bullet s^{-2}$.

In Fig.6(a), the adopted values of $G$ in CODATA-2014 adjustment [5, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33] are illustrated according
to [2]. To visualize our result (Eq.(23)), we plot $|\Delta\nu_{-}|$ as a
function of $G$ in Fig.6(b)-(c), where the parameters used are
$\varepsilon={10}^{-2}\omega_{r}$, and $\omega_{r}=2\times{10}^{9}Hz$. In (b),
we utilize a Stainless steel sphere with $\rho=9.8\times{10}^{3}Kg/m^{3}$, a
Lead one with $\rho=11.3\times{10}^{3}Kg/m^{3}$ and a Tungsten one with
$\rho=19.35\times{10}^{3}Kg/m^{3}$. On the contrary, only a Tungsten sphere is
considered in (c), where $G$ takes values of
$6.67408(31)\times{10}^{-11}m^{3}\bullet{kg}^{-1}\bullet s^{-2}$ which is the
2014 CODATA recommended value of $G$.

##  IV Summary and outlook

In sum, this paper presents a novel method which can be considered in the
field of precision measurement of $G$. In this method, two cavities and two
membrane resonators constitute an optomechanical system. Gravitational force
gradient perturbs the resonant frequencies of two membranes, resulting the
shifts of the eigenfrequencies of the supermodes. Based on EP enhanced shifts,
we can perform the measurement of $G$. Compared to traditional measurements
methods [34, 35] such as torsion balance, atom interferometry, etc., our
method possesses two distinct characters. They are: a. The proposed setup is a
minute-sized optomechanical system. b. EP effect plays an important role in
the measurement of $G$.

Based on this work where a second-order EP feature emerges in the proposed
system, we can attempt to construct a new system which can be operated at
higher-order EP, and a higher sensitivity may be achieved [36, 13, 37]. So
far, EP effect in optomechanical systems have been observed [38, 39],
indicating that our method can be put into consideration. Moreover, [40] can
be referred to for the detection of the shift of the supermode eigenfrequency.
Finally, we expect our proposal can promote the application of EP-based
sensors to the measurement of $G$.

###### Acknowledgements.

This work was supported by Natural Science Foundation of Shanghai (No.
20ZR1429900).

## References

  * Xue  _et al._ [2020] C. Xue, J. P. Liu, Q. Li, J. F. Wu, S. Q. Yang, Q. Liu, C. G. Shao, L. C. Tu, Z. K. Hu, and J. Luo, Precision measurement of the newtonian gravitational constant, National Science Review (2020).
  * Liu  _et al._ [2018] J.-P. Liu, J.-F. Wu, Q. Li, C. Xue, D.-K. Mao, S.-Q. Yang, C.-G. Shao, L.-C. Tu, Z.-K. Hu, and J. Luo, Experimental progresses on precision measurement of gravitational constant g, Acta Physica Sinica 67, 160603 (2018).
  * Gillies [1987] G. T. Gillies, The newtonian gravitational constant: an index of measurements, Metrologia 24, 1 (1987).
  * Li  _et al._ [2018] Q. Li, C. Xue, J.-P. Liu, J.-F. Wu, S.-Q. Yang, C.-G. Shao, L.-D. Quan, W.-H. Tan, L.-C. Tu, Q. Liu, _et al._ , Measurements of the gravitational constant using two independent methods, Nature 560, 582 (2018).
  * Mohr  _et al._ [2016] P. J. Mohr, D. B. Newell, and B. N. Taylor, Codata recommended values of the fundamental physical constants: 2014, Rev. Mod. Phys. 88, 035009 (2016).
  * Křen [2019] P. Křen, On the inconsistency in experimental results of the gravitational constant, Measurement 136, 622 (2019).
  * Quinn [2014] T. Quinn, Don’t stop the quest to measure big g, Nature 505, 455 (2014).
  * Heiss [2004] W. D. Heiss, Exceptional points of non-hermitian operators, Journal of Physics A: Mathematical and General 37, 2455 (2004).
  * Heiss [2012] W. D. Heiss, The physics of exceptional points, JPhA 45 (2012).
  * Wiersig [2020] J. Wiersig, Review of exceptional point-based sensors, Photon. Res. 8, 1457 (2020).
  * Wiersig [2016] J. Wiersig, Sensors operating at exceptional points: General theory, Phys. Rev. A 93, 033809 (2016).
  * Wiersig [2014] J. Wiersig, Enhancing the sensitivity of frequency and energy splitting detection by using exceptional points: Application to microcavity sensors for single-particle detection, Phys. Rev. Lett. 112, 203901 (2014).
  * Hodaei  _et al._ [2017] H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Enhanced sensitivity at higher-order exceptional points, Nature 548, 187 (2017).
  * Chen  _et al._ [2017] W. Chen, S. K. Ozdemir, G. Zhao, J. Wiersig, and L. Yang, Exceptional points enhance sensing in an optical microcavity, Nature 548, 192 (2017).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Reviews of Modern Physics 86, 1391 (2014).
  * Djorwe  _et al._ [2019] P. Djorwe, Y. Pennec, and B. Djafari-Rouhani, Exceptional point enhances sensitivity of optomechanical mass sensors, Physical Review Applied 12, 024002 (2019).
  * Djorwe  _et al._ [2018] P. Djorwe, Y. Pennec, and B. Djafari-Rouhani, Frequency locking and controllable chaos through exceptional points in optomechanics, Physical Review E 98, 032201 (2018).
  * Harber  _et al._ [2005] D. M. Harber, J. M. Obrecht, J. M. McGuirk, and E. A. Cornell, Measurement of the casimir-polder force through center-of-mass oscillations of a bose-einstein condensate, Phys. Rev. A 72, 033610 (2005).
  * Luther and Towler [1982] G. G. Luther and W. R. Towler, Redetermination of the newtonian gravitational constant g, Physical Review Letters 48, 121 (1982).
  * Karagioz and Izmailov [1996] O. Karagioz and V. Izmailov, Measurement of the gravitational constant with a torsion balance, Measurement Techniques 39, 979 (1996).
  * Bagley and Luther [1997] C. H. Bagley and G. G. Luther, Preliminary results of a determination of the newtonian constant of gravitation: a test of the kuroda hypothesis, Physical Review Letters 78, 3047 (1997).
  * Gundlach and Merkowitz [2000] J. H. Gundlach and S. M. Merkowitz, Measurement of newton’s constant using a torsion balance with angular acceleration feedback, Physical Review Letters 85, 2869 (2000).
  * Quinn  _et al._ [2001] T. Quinn, C. Speake, S. Richman, R. Davis, and A. Picard, A new determination of g using two methods, Physical Review Letters 87, 111101 (2001).
  * Kleinevoß [2002] U. Kleinevoß,  _Bestimmung der Newtonschen gravitationskonstanten G_ , Ph.D. thesis, Universität Wuppertal, Fakultät für Mathematik und Naturwissenschaften … (2002).
  * Armstrong and Fitzgerald [2003] T. Armstrong and M. Fitzgerald, New measurements of g using the measurement standards laboratory torsion balance, Physical review letters 91, 201101 (2003).
  * Hu  _et al._ [2005] Z.-K. Hu, J.-Q. Guo, and J. Luo, Correction of source mass effects in the hust-99 measurement of g, Physical Review D 71, 127505 (2005).
  * Holzschuh  _et al._ [2006] E. Holzschuh, W. Kündig, F. Nolting, R. Pixley, J. Schurr, U. Straumann, _et al._ , Measurement of newton’s gravitational constant, Physical Review D 74, 082001 (2006).
  * Luo  _et al._ [2009] J. Luo, Q. Liu, L.-C. Tu, C.-G. Shao, L.-X. Liu, S.-Q. Yang, Q. Li, and Y.-T. Zhang, Determination of the newtonian gravitational constant g with time-of-swing method, Physical Review Letters 102, 240801 (2009).
  * Tu  _et al._ [2010] L.-C. Tu, Q. Li, Q.-L. Wang, C.-G. Shao, S.-Q. Yang, L.-X. Liu, Q. Liu, and J. Luo, New determination of the gravitational constant g with time-of-swing method, Physical Review D 82, 022001 (2010).
  * Parks and Faller [2010] H. V. Parks and J. E. Faller, Simple pendulum determination of the gravitational constant, Physical Review Letters 105, 110801 (2010).
  * Quinn  _et al._ [2013] T. Quinn, H. Parks, C. Speake, and R. Davis, Improved determination of g using two methods, Physical Review Letters 111, 101102 (2013).
  * Rosi  _et al._ [2014] G. Rosi, F. Sorrentino, L. Cacciapuoti, M. Prevedelli, and G. Tino, Precision measurement of the newtonian gravitational constant using cold atoms, Nature 510, 518 (2014).
  * Newman  _et al._ [2014] R. Newman, M. Bantel, E. Berg, and W. Cross, A measurement of g with a cryogenic torsion pendulum, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20140025 (2014).
  * Rothleitner and Schlamminger [2017] C. Rothleitner and S. Schlamminger, Invited review article: Measurements of the newtonian constant of gravitation, g, Review of Scientific Instruments 88, 111101 (2017).
  * Rosi [2016] G. Rosi, Challenging the ‘big g’measurement with atoms and light, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 202002 (2016).
  * Jing  _et al._ [2017] H. Jing, Ş. Özdemir, H. Lü, and F. Nori, High-order exceptional points in optomechanics, Scientific reports 7, 1 (2017).
  * Miri and Alù [2019] M.-A. Miri and A. Alù, Exceptional points in optics and photonics, Science 363, eaar7709 (2019).
  * Zhang  _et al._ [2018] J. Zhang, B. Peng, s. K. özdemir, K. Pichler, D. O. Krimer, G. Zhao, F. Nori, Y.-x. Liu, S. Rotter, and L. Yang, A phonon laser operating at an exceptional point, Nature Photonics 12, 479 (2018).
  * Xu  _et al._ [2016] H. Xu, D. Mason, L. Jiang, and J. Harris, Topological energy transfer in an optomechanical system with exceptional points, Nature 537, 80 (2016).
  * Peng  _et al._ [2014] B. Peng, A. K. Zdemir, F. Lei, F. Monifi, M. Gianfreda, G. L. Long, S. Fan, F. Nori, C. M. Bender, and L. Yang, Parity-time-symmetric whispering-gallery microcavities, Nature Physics 10, 394 (2014).


# The galaxy 3-point correlation function: a methodological analysis.

A. Veropalumbo,11footnotetext: Corresponding author. A. Binetti  E. Branchini
M. Moresco  P. Monaco  A. Oddo  A. G. Sánchez  E. Sefusatti

###### Abstract

Upcoming galaxy surveys will provide us with an unprecedented view of the
Large-Scale Structure of the Universe and the realistic chance to extract
valuable astrophysical and cosmological information from higher-order
clustering statistics. This perspective poses new challenges, requiring both
accurate and efficient estimators and a renewed assessment of possible
systematic errors in the theoretical models and likelihood assumptions. This
work investigates these issues in relation to the analysis of the 3-point
correlation function (3PCF) in configuration space. We measure the 3PCF of 300
halo catalogs from the Minerva simulations covering a total volume of
$~{}1000h^{-3}\mathrm{Gpc}^{3}$. Each 3PCF measurement includes all possible
triangular configurations with sides between 20 and $130h^{-1}\,\mathrm{Mpc}$.
In the first place, we test different estimates of the covariance matrix, a
crucial aspect of the analysis. We compare the covariance computed numerically
from the limited but accurate benchmark simulations set to the one obtained
from $10000$ approximate halo catalogs generated with the Pinocchio code. We
demonstrate that the two numerically-estimated covariance matrices largely
match, confirming the validity of approximate methods based on Lagrangian
Perturbation Theory for generating mocks suitable for covariance estimation.
We also compare the numerical covariance with a theoretical prediction in the
Gaussian approximation. We find a good match between the two for separations
above 40$h^{-1}\,\mathrm{Mpc}$. We test the 3PCF tree-level model in
Perturbation Theory. The model is adopted in a likelihood analysis aimed at
the determination of bias parameters. We find that, for our sample of halos at
redshift $z=1$, the tree-level model performs well for separations $r\geq
40h^{-1}\,\mathrm{Mpc}$. Results obtained with this scale cut are robust
against different choices of covariance matrix. We compare to the analogous
analysis of the halo bispectrum already presented in a previous publication,
finding a remarkable agreement between the two statistics. We notice that such
comparison relies, to the best of our knowledge for the first time, on a
robust and consistent covariance estimate and on the inclusion of essentially
all measurable configurations in Fourier as in configuration space. We then
test different assumptions to build the model defining a robust combination of
hypotheses that lead to unbiased parameter estimates. Our results confirm the
importance of 3PCF, supplying a solid recipe for its inclusion in likelihood
analyses. Moreover, it opens the path for further improvements, especially in
modelling, to extract information from non-linear regimes.

##  1 Introduction

It has long been realized that the analysis of the large-scale structure (LSS)
of the Universe traced by the spatial distribution of galaxies provides a
wealth of information complementary to other cosmological probes as weak
gravitational lensing or the Cosmic Microwave Background. As a result, several
galaxy redshift surveys have been carried out and planned to probe an
increasingly larger fraction of the observable Universe, mapping the spatial
distributions of tens of millions of galaxies (e.g. DESI [1, 2], the Euclid
spectroscopic survey [3] and the Roman Telescope high-latitude survey [4]). So
far, measurements of the of two-point statistics have been the main tool for
the analysis of the LSS. Such statistics can be either the power spectrum in
Fourier space, or the 2-point correlation function in configuration space.
Yet, larger survey volumes will allow more precise measurements of higher-
order clustering statistics, able to quantify the non-Gaussian nature of
cosmological perturbations at low redshift. Such measurements can therefore
help to disentangle different sources of non-Gaussianity, either of intrinsic
interest if it’s of primordial origin or possibly responsible for relevant
degeneracies with cosmological parameters in the case of non-linear galaxy
bias.

Indeed, the analysis of higher-order statistics in recent datasets is becoming
routine, particularly in the case of the galaxy bispectrum, the 3-point
correlation function in Fourier space [5, 6, 7, 8, 9, 10]. The bispectrum, in
fact, presents some advantages with respect to its configuration-space
counterpart. In the first place, the computational cost of measuring the
bispectrum of the largest galaxy catalogs is moderate, allowing one to
accurately compute its covariance matrix if an adequate number of simulated
catalog is available [11, 12, 5, 13, 14, 15, 16]. In the second place,
theoretical models for the bispectrum, in some cases even up to one loop in
Perturbation Theory, have been extensively validated with large, simulated
data-sets [15, 17, 18, 19, 20, 21, 22, 23].

Similar developments have taken place for three-point correlation function
(3PCF) in configuration space at a relatively slower pace due to the higher
computational cost of its estimation in large galaxy catalogs and simulations
([24, 25, 26, 27, 28] and references therein). In this respect, the situation
has changed dramatically with the advent of a new class of efficient 3PCF
estimators [29, 30, 31, 32] that reduced the computational scaling cost from
$N^{3}$ to $N^{2}$, where $N$ is the number of objects in the catalog. This
allowed the analysis of large datasets and for the first detection of baryonic
acoustic oscillation features in the 3PCF of galaxies [33] and clusters [34],
and to perform a joint two- and three-point correlation analysis [35]
paralleling those performed in Fourier space [5]. Such advances triggered a
development on the theoretical side that pushed the 3PCF model beyond the
analysis of the sole monopole, taking advantage of its redshift-space
anisotropy [33, 36, 37].

The main advantage of the 3PCF over the bispectrum is its ability to deal with
survey geometry. The footprint of spectroscopic surveys is characterized by
irregular boundaries and holes in correspondence to regions that could not be
observed. A complex survey geometry can be naturally dealt with in
configuration space by using appropriate estimators as, in the 3PCF case, the
Szapudi-Szalay one [38]. This is not the case in Fourier space, where the
survey footprint induces mode coupling that needs to be accounted for by
either de-convolving the measured correlation signal or by convolving the
theory prediction with the window [31, 39, 40].

All these considerations justify the renewed interest in the 3PCF statistics
and the urgency to carry out a dedicated study aimed at investigating how
various methodological choices in the 3PCF analysis can affect the inference
of model parameters. In this work we will assess the effect of these choices
on the determination of galaxy bias parameters, one of the traditional
motivations to study the 3PCF.

The most relevant one is the estimation of the covariance matrix. Numerical
estimates of the 3PCF covariance matrix that rely on individual measurements
from a set of mock catalogs are time-consuming, both because of mocks
production and 3PCF measurements, despite the fast estimator mentioned above.
A number of workarounds have been proposed to reduce the computational cost
while keeping the accuracy at an acceptable level. Computational efficiency
can be increased either by speeding up the creation of mocks by means of
approximate halo catalogs (see [41] for a recent review) or by approximating
the covariance matrix itself. In this work we will explore both possibilities.
First of all, we will compare the 3PCF covariance matrix estimated from
catalogs of dark matter halos obtained from the Minerva set of full N-body
simulations to the one estimated from a set of halo catalogs obtained with the
approximated method implemented in the Pinocchio code, based on Lagrangian
Perturbation Theory [42, 43]. Then we shall compare the performances of
different 3PCF covariance estimators, considering the brute-force numerical
estimate from mocks, the analytical Gaussian model of [29] and the shrinkage
method [44], i.e. a hybrid combination of numerical and analytical estimates.
To evaluate the goodness of each method we will perform the same likelihood
analysis with different covariance matrices and compare the recovered best-fit
estimates of the halo bias parameters.

We will explore as well additional aspects of the standard 3PCF analysis. One
is the selection of the triangular configurations included in the analysis. A
related one is the choice of the separation bin size, defining such triangles.
In both cases we will evaluate the impact on the determination of constraints
of the bias parameters. Finally, we will consider different relations among
bias parameters, either theoretically-motivated or based on fits to numerical
simulations, as a way to reduce the parameter space.

Our results will ultimately provide an extensive assessment of the validity of
the tree-level PT model of the 3PCF under a variety of assumptions. The
posterior constraints on bias parameters will be compared to the similar
results obtained in [15] from the halo bispectrum measured from the same data
set. This constitutes, to the best of our knowledge, the most direct and fair
comparison between the two statistics presented so far in the literature,
based, on both sides, on all measurable configurations, as a function of the
smallest scale included, and taking advantage of a robust estimate of the
covariance properties from the same set of approximate mocks.

The paper is structured as follows: in Sec. 2 we briefly describe the two halo
catalogs used in our analysis and their parent simulations. In Section 3 we
describe the estimators used to measure the 3PCF and its covariance matrix,
then we compare the covariance matrices numerically estimated from the two
different sets of mock catalogs obtained from different types of simulations.
In Sec. 4 we present the model for the halo 3PCF adopted in this work and
used, along with the measured 3PCF and its covariance matrix, to perform the
likelihood analysis described in Sec. 5 and to estimate the halo bias
parameters. The results of this analysis, aimed at assessing the impact of the
various methodological options described above are presented in Sec. 6 and
discussed in Sec. 7, in which we draw our main conclusions.

##  2 Data sets

All results presented in this work are based on two different sets of halo
catalogs. The first is constituted by 298 catalogs of dark matter halos
obtained from the Minerva set of full N-body simulations. The second is a much
larger set of 10,000 approximate halo catalogs produced with the Pinocchio
code. Both sets, sharing the same simulation box and cosmological parameters,
are described below.

###  2.1 Minerva halo catalogs

The Minerva simulations [45, 46] have been performed with the GADGET-III code
[47] evolving $1000^{3}$ dark matter particles of mass $m_{p}\simeq 2.67\cdot
10^{10}\,h^{-1}\mathrm{M}_{\odot}$ in a cubic box of side
$L=1500\,h^{-1}\,\mathrm{Mpc}$.

Dark matter halos were identified with a Friends-of-Friends algorithm of
linking length equal to $0.2$ times the mean inter-particle distance.
Following [15], we only consider the $z=1$ output and halos with masses larger
than $1.12\cdot 10^{13}\,h^{-1}\mathrm{M}_{\odot}$, corresponding to a mean
number density ($\bar{n}=2.13\cdot 10^{-4}\,h^{3}\mathrm{Mpc}^{-3}$), roughly
matching that of H-$\alpha$ line emission galaxies that will be observed in
wide surveys like Euclid [3].

###  2.2 Pinocchio halo catalogs

The second set of halo catalogs has been obtained with the Pinocchio code [42,
43]. Unlike the Minerva set, here the dynamical evolution of the cosmic
density field is approximated using third order Lagrangian Perturbation
Theory. Dark matter particles are grouped in halos using a procedure based on
the ellipsoidal collapse model. Since Pinocchio realisations are
computationally cheaper than N-body simulations, a larger number of mock halo
catalogs can be generated. Here we consider the same $10,000$ catalogs used by
[15] and adopt the same threshold mass, originally chosen to match the
amplitude of the halo power spectrum measured in the Minerva simulations at
large scales. It has been shown that this matching minimizes the difference
between the power spectrum and bispectrum variance estimated from the
Pinocchio mocks and the one from the Minerva simulations, both in real [15,
17] as in redshift-space [23]. We highlight that the first 300 Pinocchio mocks
share the same random seed for the initial condition of the Minerva
simulations. This allows a direct comparison that is not affected by sample
variance.

##  3 Measuring the 3PCF

###  3.1 3PCF estimator

To measure the 3PCF of the DM halos in the Pinocchio and Minerva catalogs we
use the Szapudi & Szalay estimator [38] as implemented in [29]. That estimator
relies on the local Spherical Harmonics Decomposition (SHD) of the 3PCF

$\zeta(r_{12},r_{13},\mu)\,=\sum_{\ell}^{\ell_{max}}\zeta_{\ell}(r_{12},r_{13})\,\mathcal{L}_{\ell}(\mu)\,,$
(3.1)

where $r_{12}$ and $r_{13}$ are two sides of a triangle while $\mu$ is the
cosine of the angle between them. The length of the third side is set by the
vector relation $\mathbf{r}_{12}+\mathbf{r}_{13}+\mathbf{r}_{23}=0$.
$\mathcal{L}_{\ell}(\mu)$ is the $\ell$-th Legendre polynomial. The SHD
algorithm measures the coefficients of the expansion
$\zeta_{\ell}(r_{12},r_{13})$ up to a maximum multipole $\ell_{max}$.
Therefore eq. 3.1 provides an approximated estimate of the halo 3PCF. However,
for a typical galaxy catalog, $\ell_{max}=10$ is sufficient to converge to the
correct solution [35].

For the multipoles measurements $\hat{\zeta}_{\ell}(r_{12},r_{23})$, we
consider linear bins of size $\Delta r=10\,h^{-1}\,\mathrm{Mpc}$ up to a
maximum separation $r_{max}=150\,h^{-1}\,\mathrm{Mpc}$. The standard 3PCF is
then obtained as

$\hat{\zeta}(r_{12},r_{13};r_{23})\,=\sum_{\ell}^{\ell_{max}}\hat{\zeta}_{\ell}(r_{12},r_{13})\,\widetilde{\mathcal{L}}_{\ell}(r_{12},r_{13},r_{23})\,,$
(3.2)

where $\widetilde{\mathcal{L}}_{\ell}(r_{12},r_{13},r_{23})$ is the Legendre
polynomial of order $\ell$ weighted over the triangle

$\widetilde{\mathcal{L}}_{\ell}(r_{12},r_{13},r_{23})=\frac{\displaystyle\int_{r_{12}-\Delta_{r}/2}^{r_{12}+\Delta_{r}/2}\mathrm{d}pp^{2}\displaystyle\int_{r_{13}-\Delta_{r}/2}^{r_{13}+\Delta_{r}/2}\mathrm{d}qq^{2}\displaystyle\int_{r_{23}-\Delta_{r}/2}^{r_{23}+\Delta_{r}/2}\mathrm{d}ss^{2}\Theta(p,q,s)\mathcal{L}_{\ell}(p,q,s)}{\displaystyle\int_{r_{12}-\Delta_{r}/2}^{r_{12}+\Delta_{r}/2}\mathrm{d}pp^{2}\displaystyle\int_{r_{13}-\Delta_{r}/2}^{r_{13}+\Delta_{r}/2}\mathrm{d}qq^{2}\displaystyle\int_{r_{23}-\Delta_{r}/2}^{r_{23}+\Delta_{r}/2}\mathrm{d}ss^{2}\Theta(p,q,s)}$
(3.3)

with

$\Theta(r_{12},r_{13},r_{23})=\left\\{\begin{array}[]{ c l
}1&\quad\textrm{if}\,\left|\frac{r_{12}^{2}+r_{13}^{2}-r_{23}^{2}}{2r_{12}r_{13}}\right|\leq
1\\\ 0&\quad\textrm{otherwise}\end{array}\right.$ (3.4)

and it accounts for the binning of the third side [35]. The estimator cannot
properly recover the 3PCF corresponding to isosceles configurations
($r_{12}=r_{13}$) with a finite and manageable number of multipoles. In fact,
this intrinsic limitation extends to nearly isosceles triangles where
$r_{13}-r_{12}\simeq\Delta r$. For this reason, we characterise a given set of
measurements not only in terms of the minimum separation $r_{min}$, but also
by a minimum value $\eta_{min}$ for the relative difference

$\eta\equiv\frac{r_{13}-r_{12}}{\Delta r}\,.$ (3.5)

In other terms, fixed $\Delta r$, $r_{min}$ and $r_{max}$, setting
$\eta_{min}=0$ would correspond to measure all available triangles with sides
within the given range in steps of $\Delta r$. $\eta_{min}=1$ would exclude
configurations with $r_{13}=r_{12}$, while higher values would further reduce
the triangle set. As we will see our ability to properly model the 3PCF
measured with the SHD estimator depends significantly on this parameter. In
Sec. 5 we perform several tests to assess the sensitivity of the 3PCF analysis
to the choice of both $r_{min}$ and $\eta_{min}$.

Input to the estimator are a catalog of halos and a catalog of random objects,
that we created, which are distributed in the same volume but lacking
clustering properties. In the Minerva 3PCF measurements the number of random
objects is 50 times larger than that of the halos. In the Pinocchio case we
use 10 times as many objects. In both cases a different random catalog is
generated for each measurement. For this reason, even for the low density
randoms, we do not expect a biased 3PCF estimation but only an additional
contribution to its variance. The effect on the total error budget, however,
is negligible.

Figure 1: Top panel. Blue curve: ratio between the 3PCF rms scatter measured
in the Minerva, $\sigma_{M}$, and measured in the first 298 Pinocchio mock
catalogs, $\sigma_{P}^{298}$, shown as a function of the triangle type,
identified by an index. Larger indexes identify larger triangles. Grey band:
$2-\sigma$ scatter of $\sigma_{M}/\sigma_{P}^{298}$, obtained propagating the
variance over single rms. Bottom panel. Blue curve: difference between the
3PCF rms scatter measured the Minerva, $\sigma_{M}$, and measured in the
Pinocchio mock catalogs, $\sigma_{P}$, shown as a function of the triangle
type, identified by an index. Grey band: the reference $2-\sigma$ scatter was
obtained from subsets of 298 realisations extracted from the 10,000 Pinocchio
mocks. See text for details. Figure 2: 3PCF correlation matrix measured from
the Minerva mock catalogs (upper triangular matrix) and from the 10000
Pinocchio mocks (lower triangular matrix). Indexes identify different triangle
configurations. Larger indexes indicate larger scales. The correlation
amplitude is color-coded as indicated in the color bar.

Figure 3: Top panel: correlation amplitude of the elements along a cut through
correlation matrix in correspondence of a triangle with sides
$25\,h^{-1}\,\mathrm{Mpc},95\,h^{-1}\,\mathrm{Mpc},105\,h^{-1}\,\mathrm{Mpc}$.
Black dots: Minerva estimates. Yellow curve: Pinocchio estimates. Error bars:
jackknife errors from the Minerva. Bottom panel: same quatity estimated for a
larger triangle with sides

###  3.2 Covariance matrix

The numerical covariance matrix is estimated from the 3PCFs measured in each
mock catalog as follows:

$\widehat{C}_{i,j}=\frac{1}{N_{mocks}-1}\sum_{n=1}^{N_{mocks}}\left(\hat{\zeta}_{i}^{n}-\overline{\zeta}_{i}\right)\left(\hat{\zeta}_{j}^{n}-\overline{\zeta}_{j}\right)$
(3.6)

where the indexes $i$ and $j$ identify any two triangles in the sample with
sides $\\{r_{12},r_{13},r_{23}\\}$ and
$\\{r_{12}^{\prime},r_{13}^{\prime},r_{23}^{\prime}\\}$. $\hat{\zeta}^{n}$
indicates the 3PCF measured in the $n$-th mock and $\overline{\zeta}$ the
average among the mocks

$\overline{\zeta}_{i}=\frac{1}{N_{mocks}}\sum_{n=1}^{N_{mocks}}\hat{\zeta}_{i}^{n}.$
(3.7)

Then, to compare the Minerva and Pinocchio estimates we consider as well the
correlation matrix, defined as

$c_{ij}\equiv\frac{\widehat{C}_{ij}}{\sqrt{\sigma_{i}\sigma_{j}}}$ (3.8)

where $\widehat{C}_{i,j}$ is the covariance matrix element and
$\sigma_{i}\equiv\sqrt{\widehat{C}_{i,i}}$ represents the dispersion of the
measurements for triangle $i$.

###  3.3 3PCF covariance: Minerva vs Pinocchio

In this section we provide a first comparison between the numerical covariance
estimated from the Minerva N-body simulations and the one estimated from the
Pinocchio approximate mocks.

The top panel of Fig. 1 shows the ratio between the Minerva variance and the
variance from the first 298 Pinocchio mocks sharing the same initial
conditions. Each measurement corresponds to all triangles with
$r_{min}=20h^{-1}\,\mathrm{Mpc},r_{max}=130h^{-1}\,\mathrm{Mpc},\eta=1$
ordered assuming $r_{12}<r_{13}\leq r_{23}$. The shaded area shows the
2-$\sigma$ scatter obtained from subsets of 298 realisations extracted from
the 10,000 Pinocchio mocks. The bottom panel shows instead the difference
between the Minerva variance and the Pinocchio variance this time estimated
from the full 10,000 mocks set. We notice that the difference between the two
estimates is of the order of 5% and it is well within the estimated scatter
for the variance based on 298 measurements.

In Fig. 2 we compare instead the two correlation matrices $c_{i,j}$, defined
in eq. (3.8), with the Minerva one in the upper-left triangle and the
Pinocchio one from the 10,000 mocks in the lower-right triangle. Small
differences in the most off-diagonal elements are barely visible.

A more quantitative comparison is shown in Fig. 3 in terms of cuts across the
correlation matrix: the black dots show a row of the Minerva correlation
matrix, with the error bars estimated from the Jackknife error, while the
yellow curve connects the elements from the same row in the Pinocchio matrix.
The 3PCF covariance matrix estimated from the 10000 Pinocchio mocks and the
one from the 298 Minerva simulations agree with each other within the
uncertainty estimated for the latter. We shall consider the Pinocchio
covariance as a reference and will use it to perform all likelihood analyses
presented in what follows.

###  3.4 3PCF covariance: Minerva vs theory

Figure 4: Difference of the scatter from the $298$ Minerva simulations,
$\sigma_{M}$ and the theoretical one, $\sigma_{T}$, blue continuous line, as a
function of the triangles, from small to large size (left to right). The grey
band represent the $2\sigma$ error on the Minerva scatter, measured via
Jackknife. The two estimates are in agreement when the triangle side are
large; for small size triangles the theoretical error is systematically
smaller than the one estimated numerically.

We now consider an analytical model for the covariance matrix, in the Gaussian
approximation, as described in [29]. An analytic expression has the obvious
advantage of being noiseless and flexible to changes in the fiducial
cosmology. Most importantly it does not require large sets of mock catalogs
and the numerically expensive measurements of the 3PCF. On the other hand, the
simple analytic model we consider here makes assumption, as the Gaussianity of
the halo distribution, that could lead to potential systematic errors and does
not take into account the specific geometry of the sample (but see [48] for
possible improvements).

To compute the analytic covariance matrix we need to assume a fiducial, non-
linear power spectrum, the covered volume and a shot noise contribution to the
power spectrum. The former is provided by the best-fit model to the Minerva
power spectrum measurements222The best-fit power spectrum has been derived and
kindly provided by Kevin Pardede using the PT model and assumptions of [17]..
The latter is set equal to the inverse of the mean halo number density in the
mocks, i.e. $1./\bar{n}=4500\,h^{-3}\mathrm{Mpc}^{3}$.

To check if the analytical, Gaussian errors are consistent with those obtained
from the mock catalogs, we show in fig. 4 the difference between the two sets
of errors (blue curve), compared with the 2-$\sigma$ error measured in the
Minerva mocks (grey band) obtained as in Fig. 1. We notice that the theory
errors are systematically smaller than those measured from the mocks, with the
mismatch decreasing with the size of the triangles. This is expected since the
adequacy of the Gaussian hypothesis increases with the scale, while it is
expected to break down on small scales due to the non-linear evolution
inducing non-Gaussian contributions to the covariance.

##  4 The halo 3PCF Model

Since one of the goals of this work is to compare the results of the 3PCF
analysis to the one of the bispectrum analyses, we adopt the model of [15] and
obtain the 3PCF model via inverse-FFT as for example in [24, 25, 49]. The
model corresponds to the simple tree-level expression in Perturbation Theory
including the contribution from gravitational instability plus the one from
non-linear, local and tidal bias [50, 51, 52]. Here we follow [49] to obtain
the pre-cyclic 3PCF part of the prediction as

$\displaystyle\zeta_{pc}(r_{12},r_{13},\mu_{23})=$ $\displaystyle
b_{1}^{3}\left(\frac{34}{21}+\frac{b_{2}}{b_{1}}-\frac{4}{3}\frac{\gamma_{2}}{b_{1}}\right)\xi_{0}(r_{12})\xi_{0}(r_{13})+$
(4.1) $\displaystyle-
b_{1}^{3}\left[\xi_{1}^{+}(r_{12})\xi_{1}^{-}(r_{13})+\xi_{1}^{-}(r_{12})\xi_{1}^{+}(r_{13})\right]\mathcal{L}_{1}(\mu_{23})$
$\displaystyle+b_{1}^{3}\left(\frac{8}{21}+\frac{4}{3}\frac{\gamma_{2}}{b_{1}}\right)\xi_{2}(r_{12})\xi_{2}(r_{13})\mathcal{L}_{2}(\mu_{23})\,,$

where the triangle sides $r_{12},\,r_{13},\,r_{23}$ and the cosine angle
$\mu_{23}$ are related through
$r_{23}^{2}=r_{12}^{2}+r_{13}^{2}-2\mu_{23}r_{12}r_{13}$ and $j_{\ell}$ are
the spherical Bessel functions.

The two-point correlation function moments in eq. (4.1)

$\displaystyle\xi_{\ell}(r)$
$\displaystyle=\frac{1}{2\pi^{2}}\int_{0}^{\infty}\mathrm{d}kk^{2}P(k)j_{\ell}\left(kr\right),$
$\displaystyle\xi_{\ell}^{[\pm]}(r)$
$\displaystyle=\frac{1}{2\pi^{2}}\int_{0}^{\infty}\mathrm{d}kk^{2}k^{\pm
1}P(k)j_{\ell}\left(kr\right)$ (4.2)

depend on the linear matter power spectrum $P(k)$, that we model using CAMB
[53].

The full galaxy 3PCF model is then obtained by summing permutations of the
$\zeta_{pc}$ function:

$\zeta_{h}(r_{12},r_{13},r_{23})=\zeta_{pc}(r_{12},r_{13},\mu_{23})+\zeta_{pc}(r_{13},r_{23},\mu_{12})+\zeta_{pc}(r_{23},r_{12},\mu_{13}).$
(4.3)

To speed up the sampling of the parameter space in the likelihood analysis of
Sec. 5 we first numerically evaluate and tabulate the integrals in eq. (4) in
the reference $\Lambda$CDM cosmology and then maximize the likelihood with
respect to the bias parameters.

This model, however, cannot be compared directly to the quantities measured by
our estimator based on the Spherical Harmonics Decomposition. We should
therefore, in the first place, expand it in multipoles. In addition, we need
to account for the finite size of separation bins and we do that at the
multipoles level, computing

$\widetilde{\zeta}_{\ell}(r_{12},r_{13})=\frac{2l+1}{2}\int_{-1}^{1}\mathrm{d}\mu\mathcal{L}_{\ell}(\mu)\int_{r_{12}-\Delta
r/2}^{r_{12}+\Delta
r/2}\\!\\!\\!\\!\\!\\!\\!\mathrm{d}q\,q^{2}\int_{r_{13}-\Delta
r/2}^{r_{13}+\Delta
r/2}\\!\\!\\!\\!\\!\\!\mathrm{d}p\,p^{2}\,\zeta(p,q,\mu)\,,$ (4.4)

where $r_{12},r_{13}$ represent the center of the bins, $\Delta r$ is the bin
size and $\zeta(p,q,\mu)$ is the model, eq. (4.3). In the end, the binned 3PCF
model that we compare with the data is

$\widetilde{\zeta}(r_{12},r_{13},r_{23})=\sum_{l=0}^{{\ell}_{max}}\widetilde{\zeta}_{\ell}(r_{12},r_{13})\,\widetilde{\mathcal{L}}_{\ell}(r_{12},r_{13},r_{23}).$
(4.5)

We notice that the comparison of the model to an accurate estimation of the
full 3PCF takes advantage of a clear dependence on the physical scales
represented by the three sides. This allows in turn to a proper definition of
the minimal separation $r_{min}$ at which the model fails. The accuracy of the
estimation of the full 3PCF $\zeta(r_{12},r_{13},r_{23})$ from the equation
above, however, is limited by the finite value of $\ell_{max}=10$, required if
we want to preserve the efficiency of the implementation of [29]. In section 6
we will explore in detail the subset of triangular configurations affected by
this approximation and therefore not properly described by the tree-level
expression.

##  5 Parameter inference

We validate the galaxy 3PCF model, in terms of the posterior distribution on
the three bias parameters, testing various assumptions on the likelihood
analysis.

For each catalog, identified by the $\alpha$ index, we define the $\chi^{2}$
function

$\chi^{2}_{\alpha}(\mathbf{\mu},\mathbb{C})=(\mathbf{d}_{\alpha}-\boldsymbol{\mu})^{T}\mathbb{C}^{-1}(\mathbf{d}_{\alpha}-\boldsymbol{\mu}),$
(5.1)

where $\mathbf{\mu}$ and $\mathbf{d}_{\alpha}$ represent, respectively, the
model and data vector while $\mathbb{C}$ is the covariance matrix with
elements $C_{i,j}$. In the case of the (noiseless) analytic covariance matrix
we adopt, for the single realization, the Gaussian likelihood function

$\mathcal{L}_{\alpha}\propto\exp\left(-\frac{\chi^{2}_{\alpha}}{2}\right)\,.$
(5.2)

We will also consider covariance matrices, $\widehat{\mathbb{C}}$, estimated
numerically from a limited number of $N_{m}$ halo catalogs, whose statistical
noise can bias the determination of the precision matrix $\mathbb{C}^{-1}$. We
then include the correction suggested in [54] and given by

$\mathbb{C}^{-1}=\frac{N_{m}-n_{d}-2}{N_{m}-1}\,\widehat{\mathbb{C}}^{-1}\,$
(5.3)

where $N_{m}$ is the number of independent measurements used to estimate
$\widehat{\mathbb{C}}$ while $n_{d}$ is the size of the data vector. The
correction inflates the errors in order to encompass any possible systematics
on $\mathbb{C}^{-1}$. In the case of the covariance estimated from the full
set of 10,000 Pinocchio mocks, given a data vector comprising at most 184
triangular configurations, the correction is rather small and negligible. This
is not the case, naturally, for the covariance estimated from the more limited
set of Minerva simulations, where the factor can be as low as 0.35. We
considered as well the prescription of [55] finding marginal differences in
our applications.

The total likelihood adopted for each evaluation is obtained as the product of
all individual likelihood, or

$\log{\mathcal{L}_{tot}}=\sum_{\alpha=1}^{N_{m}}\log{\mathcal{L}_{\alpha}}\,.$
(5.4)

Finally, we use the correction factor in [56] to fully take into account the
propagation of covariance matrix uncertainty on the parameter posterior
distributions. The factor is

$m_{1}=\frac{1+B(n_{d}-n_{p})}{1+A+B(n_{p}+1)};$ (5.5)

with $n_{p}$ the size of the model parameter vector, and

$\displaystyle A$ $\displaystyle=\frac{2}{(N_{m}-n_{d}-1)(N_{m}-n_{d}-4)},$
(5.6) $\displaystyle B$
$\displaystyle=\frac{N_{m}-n_{d}-2}{(N_{m}-n_{d}-1)(N_{m}-n_{d}-4)}\,.$ (5.7)

We adopt flat, uninformative priors, specified in Table 1. To sample the
likelihood function we use a Monte Carlo Markov Chain (MCMC) approach as
implemented in the emcee software [57].

##  6 Results

Name | Binning | Covariance | $b_{1}$ | $b_{2}$ | $\gamma_{2}$ | Sec.
---|---|---|---|---|---|---
Reference | Integrated | Pinocchio | $\mathcal{U}(0.5,5)$ | $\mathcal{U}(-10,10)$ | $\mathcal{U}(-10,10)$ | 6.1
Benchmark | Integrated | Minerva | $\mathcal{U}(0.5,5)$ | $\mathcal{U}(-10,10)$ | $\mathcal{U}(-10,10)$ | 6.2
Gaussian | Integrated | Theoretical | $\mathcal{U}(0.5,5)$ | $\mathcal{U}(-10,10)$ | $\mathcal{U}(-10,10)$ | 6.3
Shrinkage | Integrated | Shrinkage | $\mathcal{U}(0.5,5)$ | $\mathcal{U}(-10,10)$ | $\mathcal{U}(-10,10)$ | 6.4
Bin-center | Center | Pinocchio | $\mathcal{U}(0.5,5)$ | $\mathcal{U}(-10,10)$ | $\mathcal{U}(-10,10)$ | 6.5
$\gamma_{2}(b_{1})$ | Integrated | Pinocchio | $\mathcal{U}(0.5,5)$ | $\mathcal{U}(-10,10)$ | Eq. 6.7 | 6.6
$b_{2}(b_{1})$ | Integrated | Pinocchio | $\mathcal{U}(0.5,5)$ | Eq. 6.8 | $\mathcal{U}(-10,10)$ | 6.6
Table 1: Summary of test case and their main characteristics.

In this section we will consider the effect of different triangle selections,
covariance matrices and binning schemes, as well as the inclusion of relations
among bias parameters to reduce the parameter space. Table 1 lists the
characteristics of the various tests and summarizes the methodological choices
adopted in each of them. All runs share the same Minerva halo catalogs at
$z=1$ where the 3PCF is measured with the SHD estimator using $\ell_{max}=10$,
$\Delta r=10h^{-1}\,\mathrm{Mpc}$ and $\eta_{min}=1$. The Gaussian likelihood
is maximized with respect to the three halo bias parameters $b_{1},\,b_{2}$
and $\gamma_{2}$. All other parameters are kept fixed.

###  6.1 Sensitivity to triangle size and shape

|  | $r_{min}\,[h^{-1}\,\mathrm{Mpc}]$
---|---|---
|  | 10 | 20 | 30 | 40 | 50 | 60 | 70
$\eta_{min}$ | 1 | 216 | 186 | 152 | 117 | 84 | 56 | 35
2 | 170 | 143 | 113 | 83 | 56 | 35 | 20
3 | 131 | 107 | 81 | 56 | 35 | 20 | 10
Table 2: Number of triangles as a function of $\eta_{min}$ and $r_{min}$.
$r_{max}$ is fixed to $130\,h^{-1}\,\mathrm{Mpc}$.

Figure 5: Inferred bias parameters as a function of the minimum triangle side
$r_{min}$ (thick red curve) with their $1\sigma$ uncertainty (shaded area) for
the Reference case. Dashed black lines with grey shaded region represent the
bias values and their $1\sigma$ errors estimated from the bispectrum analysis
of [15]. Top, central and bottom panels: all results are for the reference
case using triangles with $\eta\geq 1,\,2$ and 3, respectively. Figure 6:
$1D-2D$ posterior probability distributions of the bias parameters inferred
from the likelihood analysis in the reference case. The green, blue and red
contours (curves) are used for the cases the cases
$r_{min}=20,30,40\,h^{-1}\,\mathrm{Mpc}$ and $\eta_{min}=3$. Reference dark
grey contours (curves) are from the bispectrum analysis of [15].

To assess the sensitivity of the 3PCF analysis to the size and shapes of the
triangles we performed different likelihood analyses in which we fix
$\eta_{min}$, hence selecting the triangle types, and progressively increase
the value of $r_{min}$. The number of triangles that remain after a specific
selection is shown in Tab. 2.

The results are shown in Fig. 5. In each panel we plot the best fit bias
parameters ($b_{1}$, $b_{2}$, $\gamma_{2}$, from top to bottom) as a function
of $r_{min}$. Shaded areas represent 1-$\sigma$ uncertainty. We also plot, for
reference, the best fit values (and uncertainties) obtained by [15] from the
bispectrum analysis of the same dataset with Fourier space triangles smaller
than $k_{max}=0.082\,h\,{\rm Mpc}^{-1}$ (dashed black lines and grey shaded
areas). Each panel shows the results obtained with $\eta_{min}=1,2$ and $3$
(red, green and blue curves, respectively), corresponding to excluding those
triangles with no or a small difference between the sides $r_{12}$ and
$r_{13}$, see eq. (3.5).

Let us focus first on the cases with $\eta_{min}=2$ and 3. In these cases,
both parameters $b_{2}$ and $\gamma_{2}$, for $r_{min}\geq
20\,h^{-1}\,\mathrm{Mpc}$, converge to a constant value consistent with the
one obtained in the bispectrum analysis. Only the linear bias parameters
$b_{1}$ exhibits a scale dependence, particularly significant for $r_{min}$ up
to $\sim 40h^{-1}\,\mathrm{Mpc}$ and for $\eta_{min}=2$.

The case of $\eta_{min}=1$ is different. In this case $\gamma_{2}$ and $b_{1}$
show some dependence on $r_{min}$; the former converges to the correct value,
but we should notice how the behavior of $b_{2}$ shows no clear trend and it
is inconsistent with the bispectrum results for all values of $r_{min}$. This
can be explained recalling that the case $\eta_{min}=1$ includes the quasi-
isosceles configurations, that is those with $r_{12}\simeq r_{13}$, with the
third side $r_{23}$ approaching zero. When this happens, the 3PCF tree-level
model that we adopted in our analysis, specific to the fast estimator
described in section 3.1, fails [29, 35]. In fact, both the estimator and the
model are defined by a maximum value of $\ell_{max}=10$ in the Legendre
expansion, eq. (3.1), and probably a larger value of $\ell_{max}$ would have
reduced the dependence on $\eta_{min}$. In the rest of this work we adopt
$\eta_{min}=3$ for all our results.

To inspect the degeneracy between the bias parameters and how this depends on
$r_{min}$, we plot their 1D and 2D marginalized distributions in the triangle
plot of Fig. 6. We compare the three cases of $r_{min}=20,30$ and
$40h^{-1}\,\mathrm{Mpc}$ (respectively in blue green and red) with the result
of the bispectrum analysis (grey) from [15], where a scale cut of
$k_{max}=0.09\,h\,{\rm Mpc}^{-1}$ is assumed.

For $r_{min}\leq 30h^{-1}\,\mathrm{Mpc}$ the constraints from the 3PCF
analysis are slightly tighter and less degenerate than in the bispectrum case,
with the $b_{1}$ value recovered from the 3PCF 1-$\sigma$ lower than the one
from the bispectrum. On the other hand, the error on $b_{2}$ inferred from the
bispectrum is smaller. Yet, the overall agreement is an indication that our
understanding of halo 3-point correlation in terms of local and tidal bias
along with the control of systematic in both the Fourier-space as the
configuration-space estimations and analyses are on rather solid grounds.

###  6.2 Noisy estimates of the numerical covariance

In Section 3.3 we found no significant difference between the 3PCF covariance
matrices estimated from the Minerva simulations and the Pinocchio approximate
mocks. We can nevertheless expect a bias, particularly on the inverse
covariance matrix, when the numerical estimation is based on a small number of
catalogs. As described in Sec. 5, we account for this enlarging the error bars
to encompass possible systematic effects.

Figure 7: Marginalized posterior distribution of three bias parameters
assuming $\eta_{min}=3$, $r_{min}=20$ (top panels) and $r_{min}=40$ (bottom
panels), obtained from different covariance matrices: numerical estimates from
the full Pinocchio set (black continuous curve), small set of 298 catalogs
from the Minerva simulations (red dot-dashed), 298 catalogs from the Pinocchio
mocks (green dotted) and the analytic Gaussian model (blue dashed).

In Fig. 7 we compare the marginalized posteriors on the bias parameters
obtained assuming the covariance matrix estimated from the small sets of 298
Minerva simulations (red curves) and the one estimated from the same number of
Pinocchio mocks (green) to the reference case of the covariance from the full
Pinocchio set. (black). The top panels assume
$r_{min}=20h^{-1}\,\mathrm{Mpc}$, while for the lower ones the data vector is
limited to $r_{min}=40h^{-1}\,\mathrm{Mpc}$. For
$r_{min}=20h^{-1}\,\mathrm{Mpc}$ the main difference is clearly due to the
correction factor applied to the poorly estimated covariances. In this case we
notice, in addition, a shift of about 1-$\sigma$ in the mean of the $b_{2}$,
sign of a systematic error induced by a poor estimation of the covariance. We
also notice a shift in the mean values between the 298 Minerva and 298
Pinocchio cases, despite they share the same initial conditions. This suggests
that noisy covariance obtained with approximated methods should be taken with
care. This discrepancy disappears when excluding triangles with sides smaller
that 40 $h^{-1}\,\mathrm{Mpc}$ from the analysis. Similar results are obtained
with different values for $\eta_{min}$ values and for larger $r_{min}$ sizes
and for this reason they are not shown here.

###  6.3 Analytical covariance model

Figure 8: Value of the shrinkage intensity $\lambda^{*}$ as a function of the
minimum triangle scale. Different curves correspond to different choices of
the numerical matrix $C_{i,j}$ estimated from 298 Minerva simulations (black,
continuous) and 100, 1000 and 3500 Pinocchio mocks (respectively red dashed,
blue dot-dashed and green dotted). For this figure we select only triangles
with $\eta_{min}=3$.

We turn now to the effect of the analytical covariance on the on the
likelihood analysis. The resulting posterior distributions of the bias
parameters are shown as blue dashed curves in Fig. 7. Not surprisingly, in the
$r_{min}=20\,h^{-1}\,\mathrm{Mpc}$ case (top panels) the marginalized
uncertainties are underestimated w.r.t. the case of the reference numerical
covariance. However, the difference is not very large and we do not notice a
significant bias on the central values of the parameters. The comparison
improves significantly for $r_{min}=40\,h^{-1}\,\mathrm{Mpc}$ (bottom panels),
indicating that the Gaussian hypothesis is rather adequate for 3PCF analyses
above this scale.

###  6.4 Shrinkage covariance matrix

Several methods have been proposed to estimate the covariance matrix with a
limited number of mocks while keeping a good accuracy [44, 58]. We consider
here the shrinkage technique of [44], based on the combination of a noisy but
unbiased estimate of the covariance (e.g. from N-body simulations) with a
possibly biased but more precise estimate or noiseless model (e.g. from
theory). The benefit of this combination is that of reducing the number of
mocks but also, in our case, the non-negligible computational cost of
estimating their 3PCFs.

In practice, we consider in the linear combination the (biased but noiseless)
theoretical, Gaussian covariance $T_{i,j}$ with the numerical covariance,
$\widehat{C}_{i,j}$, that is

$S_{i,j}=\lambda^{*}T_{i,j}+(1-\lambda^{*})\widehat{C}_{i,j}\,,$ (6.1)

where the parameter $\lambda^{*}$ represent the shrinkage intensity. For
$\lambda^{*}=0$ one recover the numerical covariance while for $\lambda^{*}=1$
we have the theory model. The optimal choice for $\lambda$ is found by
minimizing a mean-square difference cost function giving, in our case

$\lambda^{\star}=\frac{\sum_{i,j}\widehat{\operatorname{Var}}\left(\widehat{C}_{ij}\right)}{\sum_{i,j}\left(T_{ij}-\widehat{C}_{ij}\right)^{2}}\,,$
(6.2)

where $\widehat{\operatorname{Var}}\left(\widehat{C}_{ij}\right)$ is the
variance associated with each $\widehat{C}_{i,j}$, that we compute as [44]

$\widehat{\operatorname{Var}}\left(\widehat{C}_{ij}\right)=\frac{N_{m}}{(N_{m}-1)^{3}}\sum_{\alpha=1}^{N_{m}}\left(W_{ij}^{(\alpha)}-\widehat{W}_{ij}\right)^{2}\,,$
(6.3)

with

$W_{ij}^{(\alpha)}=\left(\hat{\zeta}_{i}^{(\alpha)}-\overline{\zeta}_{i}\right)\left(\hat{\zeta}_{j}^{(\alpha)}-\overline{\zeta}_{j}\right)\,$
(6.4)

computed in terms of the 3PCF measurement $\hat{\zeta}_{i}^{(\alpha)}$ for the
$i$-th triangle in the $\alpha$-th realization and the average
$\overline{\zeta}_{i}$ at triangle $i$, and with

$\widehat{W}_{ij}=\frac{1}{N_{m}}\sum_{\alpha=1}^{N_{m}}W_{ij}^{(\alpha)}\,.$
(6.5)

We show our estimates of $\lambda^{*}$ as a function of $r_{min}$, assuming
$\eta_{min}=3$, in fig. 8. Different curves correspond to different choices of
the numerical matrix $\widehat{C}_{i,j}$. As expected, on large scales
$\lambda^{*}\rightarrow 1$ since the covariance is approximately Gaussian,
although when 1000 or more mocks are available, the parameter $\lambda^{*}$
never reaches unity, showing relevant differences between the numerical
estimate and analytical predictions at all scale. At small scales, where
deviations from Gaussianity are larger, the numerical covariance is always
favoured. In the Minerva case, we also run the likelihood analysis with the
shrinkage covariance. Unsurprisingly, the same trend found for $\lambda^{*}$
reflects on the constraints on the bias parameters, which stay close to the
constraints shown in Fig. 7.

###  6.5 Sensitivity to binning

Figure 9: Percentage difference between bias constraints obtained with and
without averaging over the triangle side bin, as a function of $r_{min}$ for
$\eta\geq 3$. The black dashed line drawn at the 0 value (perfect match)
indicates when the two hypothesis give consistent constraints.

In all the results presented so far we have used the 3PCF multipole model
predictions obtained by performing the three-dimensional integral in eq.
(4.4). This step significantly contributes to the total computing budget and
becomes particularly burdensome when performed at each likelihood evaluation
in a Markov Chain exploring cosmological parameters. An alternative, cheaper
option could be the simple angle integration

$\zeta_{\ell}(r_{12},r_{13})=\frac{2\ell+1}{2}\int_{-1}^{1}\mathrm{d}\mu\mathcal{L}_{\ell}(\mu)\zeta(r_{12},r_{13},r_{23}),$
(6.6)

with $r_{12}$ and $r_{13}$ evaluated at the center of the separation bin.

In Fig. 9 we compare the results of the likelihood analysis obtained in these
two cases. We show, for each bias parameters, the percent difference relative
to the reference case in which the binning effect has been fully modelled,
along with the relative marginalized uncertainty. The simplification given by
the evaluation at the bin center can induce significant deviations from the
reference values, as large as 50% in the case of $b_{2}$, and for all values
of $r_{min}$ up to 70$h^{-1}\,\mathrm{Mpc}$.

We acknowledge the fact that the effect may be particularly large due to the
bin size of $\Delta r=10\,h^{-1}\,\mathrm{Mpc}$, and that using finer binning
will reduce the mismatch. That said, our results show that the analysis is
sensitive to the binning adopted and that an alternative definition eq. (6.6)
to speed up the analysis should be considered carefully.

###  6.6 Bias relations

Figure 10: Best fit values of the bias parameters as a function of $r_{min}$
obtained in the reference case (red continuous curve), assuming the LLB
hypothesis (dashed green curves) and using the N-body calibrated fitting
function (blue, dot-dashed curve). Dashed areas represent the 1-$\sigma$
uncertainty interval. All results are obtained for triangles with $\eta\geq
3$.

So far we assume all three bias parameters free to vary. However, these
parameters are not fully independent, and several relations among them can be
found either from theoretical considerations or as fits to numerical
simulations. Here we consider specifically the local Lagrangian bias (LLB)
prescription for the tidal parameter [51, 52]

$\gamma_{2}(b_{1})=-\frac{2}{7}(b_{1}-1)\,,$ (6.7)

and the numerical fit [59]:

$b_{2}(b_{1})=0.412-2.143\,b_{1}+0.929\,b_{1}^{2}+0.008\,b_{1}^{3}+\frac{4}{3}\gamma_{2}\,.$
(6.8)

Fig. 10 shows the mean over marginalized posterior of bias parameters of the
best fit and their 1-$\sigma$ uncertainties as a function of $r_{min}$ for the
case of $\eta_{min}=3$. The results obtained assuming the LLB prescription
(dashed green curves) are not compatible with those obtained in the reference
case (red continuous curves), especially for the derived $\gamma_{2}$
parameter, confirming previous results [60, 61, 15]. On the other hand, the
posteriors obtained assuming the N-body-calibrated relation in eq. (6.8) (dot-
dashed blue curves) are in very good agreement with the reference ones for all
$r_{min}$ values. A similarly good agreement has been found in the bispectrum
analysis of [15].

##  7 Discussion and conclusions

This work extends to the 3-point correlation function the detailed analysis,
previously presented for the power spectrum and bispectrum in [15, 21, 62, 17,
40, 23], of a large set of 298 halo catalogs from the Minerva simulations
[45], corresponding to a total volume of about 1,000 $\,h^{-3}\,{\rm
Gpc}^{3}$. These are supplemented with an even larger set of 10,000
approximate halo catalogs obtained with the Pinocchio code on the same box and
cosmology and reproducing the large-scale amplitude of the simulations’ halo
power spectrum. While not directly corresponding to any galaxy population from
actual surveys, these large data-sets allow for rigorous test of several
methodological choices in the analysis of galaxy clustering, taking advantage
of a robust estimate of their covariance properties, even in the case of
higher-order statistics.

We measure the real-space 3PCF from all the Minerva halo catalogs as well as
from the Pinocchio mocks. Each measurement includes all measurable triangular
configurations with sides between 20 and 130$h^{-1}\,\mathrm{Mpc}$, with only
a small subset excluded due to the limitations of the fast 3PCF estimator of
[29]. This is crucial to explore the potential of the 3PCF in constraining
bias and cosmological parameters. In our tests, we focus on the determination
of the three bias parameters required by the tree-level model in Perturbation
Theory.

  * •

We explored how the recovered bias parameter depend on the value of the
smallest separation $r_{min}$ included in the analysis as well as on the
minimal, relative difference $\eta_{min}$ between the sides $r_{12}$ and
$r_{13}$, defined in eq. (3.5). For $\eta_{min}=1$ the sum of eq. (4.5)
limited to $\ell_{max}=10$ does not provide an accurate estimate of the full
3PCF $\zeta$ and consequently is not properly described by the tree-level
model. We find that setting $\eta_{min}=3$ leads to a consistent determination
of the linear bias $b_{1}$ and the non-linear parameters $b_{2}$ and
$\gamma_{2}$, with values independent of $r_{min}$ for $r_{min}\geq
20\,h^{-1}\,\mathrm{Mpc}$. Some dependence of $b_{1}$ on $r_{min}$ is evident
instead for $\eta_{min}=2$, while $\eta_{min}=1$ leads to largely inconsistent
results for all bias parameters.

  * •

We compare our constraints on the bias parameters for $\eta_{min}=3$ (assumed
as reference value in all that follows) to the same constraints obtained from
the bispectrum analysis of [15]. We find a remarkable agreement, particularly
for the scale cuts given by $r_{min}=40h^{-1}\,\mathrm{Mpc}$ and
$k_{max}=0.09\,h\,{\rm Mpc}^{-1}$ for the 3PCF and bispectrum, respectively.
This validates both the theoretical description of the large-scale 3-point
statistics in terms of local and tidal bias contributions as well as our
accurate estimation and modelling. This is, as far as we are aware, the most
balanced comparison between the two statistics, based, on both sides, on all
measurable configurations and on a robust numerical estimate of the covariance
from a very large set of mocks.

  * •

We compare the results obtained from covariance matrix estimated from the
Pinocchio mocks to those assuming a covariance obtained from more limited set
of Minerva N-body simulations (which also provide our data-vector), accounting
for possible systematic errors due to poor statistics in terms of the factor
proposed by [54]. The results are very similar for
$r_{min}=40h^{-1}\,\mathrm{Mpc}$, while we notice some systematic shift in the
determination of $b_{2}$ for $r_{min}=20h^{-1}\,\mathrm{Mpc}$, possibly due to
small-scale suppression of power characteristic of the Lagrangian PT on which
Pinocchio, as many other similar approximate methods, are based.

  * •

We show a comparison of the same results with those obtained using an
analytical model for the covariance matrix, assuming the Gaussian
approximation [29]. Again we find largely consistent results for $r_{min}\geq
40\,h^{-1}\,\mathrm{Mpc}$, while including smaller scales we can notice a
slight underestimate of the marginalized uncertainty for all parameters,
w.r.t. the reference numerical covariance. This is expected also from a direct
comparison in terms of the variance alone, with the analytical model providing
lower values at small scales.

  * •

We investigate how the shrinkage technique [44] can provide an hybrid estimate
of the covariance based on the possibly biased analytical model and on a noisy
numerical estimate, determining the shrinkage parameter $\lambda^{*}$ as a
function of $r_{min}$ for different choices of the numerical matrix. We find
that in all cases the numerical result is preferred when $r_{min}\lesssim
30h^{-1}\,\mathrm{Mpc}$, but also that a significant numerical component is
retained when the number of available mocks is sufficiently large ($\sim
1000$).

  * •

We compare our reference analysis, fully accounting for the size of the
separation bins of the 3PCF estimator, to a more efficient evaluation that
neglects binning effects, finding large systematic errors on the bias
parameters. Although such errors can be reduced by using a smaller bin size
than the one adopted here ($10h^{-1}\,\mathrm{Mpc}$), a careful assessment has
to be made before adopting this approximation, also considering that a smaller
bin implies a much larger covariance matrix, with all its significant
numerical requirements.

  * •

Finally, we consider relations among the bias parameters in order to reduce
the parameter space. We confirmed that the simple local Lagrangian expression
for the tidal bias leads to significant systematic errors in the 3PCF case as
well. The relation between $b_{2}$ and $b_{1}$ obtained by [59] as a fit to
numerical simulations, instead leads to results consistent with our reference
case where all three parameters are left free.

With this work we have provided a solid test of 3-point correlation function
as a useful tool in the determination of second-order bias parameters, showing
constraints consistent to those from a bispectrum analysis. A fully numerical
estimates of the 3PCF covariance matrix, however, is always a very numerically
demanding task. We discussed here possible strategies to tackle this problem
and provide ideal tests based on periodic-box simulations. We aim, in the
future, to explore in further details prescriptions targeting a specific
survey geometry. Other possible extensions of this work go in the direction of
testing more sophisticated models (e.g. [63]) to extract as much information
as possible on both bias and cosmological parameters.

## Acknowledgments

We want to thank Luigi Guzzo, Elena Sarpa and Massimo Guidi for useful
discussions, Claudio Dalla Vecchia for performing and sharing the Minerva
simulations used in this work. The Minerva simulations have been performed on
the Hydra and Euclid clusters at the Max Planck Computing and Data Facility
(MPCDF) in Garching. The Pinocchio mocks were run on the GALILEO cluster at
CINECA, thanks to an agreement with the University of Trieste. AV, EB and MM
are partially supported by ASI/INAF agreement n. 2018-23-HH.0 “Scientific
activity for Euclid mission, Phase D" EB, AV, PM and ES are also supported by
INFN project “InDark". EB and AV are also supported by MIUR/PRIN 2017 “From
Darklight to Dark Matter: understanding the galaxy-matter connection to
measure the Universe". EB is also supported by ASI/INAF agreement n.
2017-14-H.O “Unveiling Dark Matter and Missing Baryons in the high-energy
sky". ES and PM are also partially supported by PRIN MIUR 2015 "Cosmology and
Fundamental Physics: illuminating the Dark Universe with Euclid." MM
acknowledges support from MIUR, PRIN 2017 (grant 20179ZF5KS).

## References

  * [1] DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen et al., _The DESI Experiment Part I: Science,Targeting, and Survey Design_ , _arXiv e-prints_ (2016) arXiv:1611.00036 [1611.00036].
  * [2] DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen et al., _The DESI Experiment Part II: Instrument Design_ , _arXiv e-prints_ (2016) arXiv:1611.00037 [1611.00037].
  * [3] R. Laureijs, J. Amiaux, S. Arduini, J. L. Auguères, J. Brinchmann, R. Cole et al., _Euclid Definition Study Report_ , _arXiv e-prints_ (2011) arXiv:1110.3193 [1110.3193].
  * [4] R. Akeson, L. Armus, E. Bachelet, V. Bailey, L. Bartusek, A. Bellini et al., _The Wide Field Infrared Survey Telescope: 100 Hubbles for the 2020s_ , _arXiv e-prints_ (2019) arXiv:1902.05569 [1902.05569].
  * [5] H. Gil-Marín, W. J. Percival, L. Verde, J. R. Brownstein, C.-H. Chuang, F.-S. Kitaura et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: RSD measurement from the power spectrum and bispectrum of the DR12 BOSS galaxies_ , _Mon. Not. R. Astron. Soc._ 465 (2017) 1757 [1606.00439].
  * [6] G. d’Amico, J. Gleyzes, N. Kokron, K. Markovic, L. Senatore, P. Zhang et al., _The cosmological analysis of the SDSS/BOSS data from the Effective Field Theory of Large-Scale Structure_ , _Journal of Cosmology and Astro-Particle Physics_ 2020 (2020) 005 [1909.05271].
  * [7] G. D’Amico, M. Lewandowski, L. Senatore and P. Zhang, _Limits on primordial non-Gaussianities from BOSS galaxy-clustering data_ , _arXiv e-prints_ (2022) arXiv:2201.11518 [2201.11518].
  * [8] O. H. E. Philcox and M. M. Ivanov, _BOSS DR12 full-shape cosmology: $\Lambda$ CDM constraints from the large-scale galaxy power spectrum and bispectrum monopole_, _Phys. Rev. D_ 105 (2022) 043517 [2112.04515].
  * [9] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović and M. Zaldarriaga, _Constraints on Single-Field Inflation from the BOSS Galaxy Survey_ , _arXiv e-prints_ (2022) arXiv:2201.07238 [2201.07238].
  * [10] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović and M. Zaldarriaga, _Constraints on Multi-Field Inflation from the BOSS Galaxy Survey_ , _arXiv e-prints_ (2022) arXiv:2204.01781 [2204.01781].
  * [11] E. Sefusatti, M. Crocce, S. Pueblas and R. Scoccimarro, _Cosmology and the bispectrum_ , _Phys. Rev. D_ 74 (2006) 023522 [arXiv: astro-ph/0604505].
  * [12] K. C. Chan and L. Blot, _Assessment of the information content of the power spectrum and bispectrum_ , _Phys. Rev. D_ 96 (2017) 023528 [1610.06585].
  * [13] M. Colavincenzo, E. Sefusatti, P. Monaco, L. Blot, M. Crocce, M. Lippich et al., _Comparing approximate methods for mock catalogues and covariance matrices - III: bispectrum_ , _Mon. Not. R. Astron. Soc._ 482 (2019) 4883 [1806.09499].
  * [14] C. Hahn, F. Villaescusa-Navarro, E. Castorina and R. Scoccimarro, _Constraining M ν with the bispectrum. Part I. Breaking parameter degeneracies_, _Journal of Cosmology and Astro-Particle Physics_ 2020 (2020) 040 [1909.11107].
  * [15] A. Oddo, E. Sefusatti, C. Porciani, P. Monaco and A. G. Sánchez, _Toward a robust inference method for the galaxy bispectrum: likelihood function and model selection_ , _Journal of Cosmology and Astro-Particle Physics_ 2020 (2020) 056 [1908.01774].
  * [16] D. Gualdi, H. Gil-Marín and L. Verde, _Joint analysis of anisotropic power spectrum, bispectrum and trispectrum: application to N-body simulations_ , _Journal of Cosmology and Astro-Particle Physics_ 2021 (2021) 008 [2104.03976].
  * [17] A. Oddo, F. Rizzo, E. Sefusatti, C. Porciani and P. Monaco, _Cosmological parameters from the likelihood analysis of the galaxy power spectrum and bispectrum in real space_ , _Journal of Cosmology and Astro-Particle Physics_ 2021 (2021) 038 [2108.03204].
  * [18] T. Baldauf, M. Garny, P. Taule and T. Steele, _Two-loop bispectrum of large-scale structure_ , _Phys. Rev. D_ 104 (2021) 123551 [2110.13930].
  * [19] A. Moradinezhad Dizgah, M. Biagetti, E. Sefusatti, V. Desjacques and J. Noreña, _Primordial non-Gaussianity from biased tracers: likelihood analysis of real-space power spectrum and bispectrum_ , _Journal of Cosmology and Astro-Particle Physics_ 2021 (2021) 015 [2010.14523].
  * [20] M. M. Ivanov, O. H. E. Philcox, T. Nishimichi, M. Simonović, M. Takada and M. Zaldarriaga, _Precision analysis of the redshift-space galaxy bispectrum_ , _arXiv e-prints_ (2021) arXiv:2110.10161 [2110.10161].
  * [21] D. Alkhanishvili, C. Porciani, E. Sefusatti, M. Biagetti, A. Lazanu, A. Oddo et al., _The reach of next-to-leading-order perturbation theory for the matter bispectrum_ , _arXiv e-prints_ (2021) arXiv:2107.08054 [2107.08054].
  * [22] A. Eggemeier, R. Scoccimarro, R. E. Smith, M. Crocce, A. Pezzotta and A. G. Sánchez, _Testing one-loop galaxy bias: Joint analysis of power spectrum and bispectrum_ , _Phys. Rev. D_ 103 (2021) 123550 [2102.06902].
  * [23] F. Rizzo, C. Moretti, K. Pardede, A. Eggemeier, A. Oddo, E. Sefusatti et al., _The Halo Bispectrum Multipoles in Redshift Space_ , _arXiv e-prints_ (2022) arXiv:2204.13628 [2204.13628].
  * [24] Y. P. Jing and G. Boerner, _Three-point correlation function in the quasilinear regime._ , _Astron. Astrophys._ 318 (1997) 667 [astro-ph/9606122].
  * [25] J. Barriga and E. Gaztañaga, _The three-point function in large-scale structure - I. The weakly non-linear regime in N-body simulations_ , _Mon. Not. R. Astron. Soc._ 333 (2002) 443 [astro-ph/0112278].
  * [26] J. Bel, K. Hoffmann and E. Gaztañaga, _Non-local bias contribution to third-order galaxy correlations_ , _Mon. Not. R. Astron. Soc._ 453 (2015) 259 [1504.02074].
  * [27] K. Hoffmann, E. Gaztañaga, R. Scoccimarro and M. Crocce, _Testing the consistency of three-point halo clustering in Fourier and configuration space_ , _Mon. Not. R. Astron. Soc._ 476 (2018) 814 [1708.08941].
  * [28] J. Kuruvilla and C. Porciani, _The n-point streaming model: how velocities shape correlation functions in redshift space_ , _Journal of Cosmology and Astro-Particle Physics_ 2020 (2020) 043 [2005.05331].
  * [29] Z. Slepian and D. J. Eisenstein, _Computing the three-point correlation function of galaxies in $O(N^{2})$ time_, _Mon. Not. R. Astron. Soc._ 454 (2015) 4142 [1506.02040].
  * [30] Z. Slepian and D. J. Eisenstein, _A practical computational method for the anisotropic redshift-space three-point correlation function_ , _Mon. Not. R. Astron. Soc._ 478 (2018) 1468 [1709.10150].
  * [31] N. S. Sugiyama, S. Saito, F. Beutler and H.-J. Seo, _A complete FFT-based decomposition formalism for the redshift-space bispectrum_ , _Mon. Not. R. Astron. Soc._ 484 (2019) 364 [1803.02132].
  * [32] O. Umeh, _Optimal computation of anisotropic galaxy three point correlation function multipoles using 2DFFTLOG formalism_ , _Journal of Cosmology and Astro-Particle Physics_ 2021 (2021) 035 [2011.05889].
  * [33] Z. Slepian, D. J. Eisenstein, J. R. Brownstein, C.-H. Chuang, H. Gil-Marín, S. Ho et al., _Detection of baryon acoustic oscillation features in the large-scale three-point correlation function of sdss boss dr12 cmass galaxies_ , _Mon. Not. R. Astron. Soc._ 469 (2017) 1738 [1607.06097].
  * [34] M. Moresco, A. Veropalumbo, F. Marulli, L. Moscardini and A. Cimatti, _C 3: Cluster Clustering Cosmology. II. First Detection of the Baryon Acoustic Oscillations Peak in the Three-point Correlation Function of Galaxy Clusters_, _Astrophys. J._ 919 (2021) 144 [2011.04665].
  * [35] A. Veropalumbo, I. Sáez Casares, E. Branchini, B. R. Granett, L. Guzzo, F. Marulli et al., _A joint 2- and 3-point clustering analysis of the VIPERS PDR2 catalogue at z 1: breaking the degeneracy of cosmological parameters_ , _Mon. Not. R. Astron. Soc._ 507 (2021) 1184 [2106.12581].
  * [36] Z. Slepian and D. J. Eisenstein, _Modelling the large-scale redshift-space 3-point correlation function of galaxies_ , _Mon. Not. R. Astron. Soc._ 469 (2017) 2059 [1607.03109].
  * [37] N. S. Sugiyama, S. Saito, F. Beutler and H.-J. Seo, _Towards a self-consistent analysis of the anisotropic galaxy two- and three-point correlation functions on large scales: application to mock galaxy catalogues_ , _Mon. Not. R. Astron. Soc._ 501 (2021) 2862 [2010.06179].
  * [38] I. Szapudi and A. S. Szalay, _A New Class of Estimators for the N-Point Correlations_ , _Astrophys. J. Lett._ 494 (1998) L41.
  * [39] O. H. E. Philcox, _Cosmology without window functions. II. Cubic estimators for the galaxy bispectrum_ , _Phys. Rev. D_ 104 (2021) 123529 [2107.06287].
  * [40] K. Pardede, F. Rizzo, M. Biagetti, E. Castorina, E. Sefusatti and P. Monaco, _Bispectrum-window convolution via Hankel transform_ , _arXiv e-prints_ (2022) arXiv:2203.04174 [2203.04174].
  * [41] P. Monaco, _Approximate Methods for the Generation of Dark Matter Halo Catalogs in the Age of Precision Cosmology_ , _Galaxies_ 4 (2016) 53 [1605.07752].
  * [42] P. Monaco, T. Theuns and G. Taffoni, _The pinocchio algorithm: pinpointing orbit-crossing collapsed hierarchical objects in a linear density field_ , _Mon. Not. R. Astron. Soc._ 331 (2002) 587 [astro-ph/0109323].
  * [43] E. Munari, P. Monaco, E. Sefusatti, E. Castorina, F. G. Mohammad, S. Anselmi et al., _Improving fast generation of halo catalogues with higher order Lagrangian perturbation theory_ , _Mon. Not. R. Astron. Soc._ 465 (2017) 4658 [1605.04788].
  * [44] A. C. Pope and I. Szapudi, _Shrinkage estimation of the power spectrum covariance matrix_ , _Mon. Not. R. Astron. Soc._ 389 (2008) 766 [0711.2509].
  * [45] J. N. Grieb, A. G. Sánchez, S. Salazar-Albornoz and C. Dalla Vecchia, _Gaussian covariance matrices for anisotropic galaxy clustering measurements_ , _Mon. Not. R. Astron. Soc._ 457 (2016) 1577 [1509.04293].
  * [46] M. Lippich, A. G. Sánchez, M. Colavincenzo, E. Sefusatti, P. Monaco, L. Blot et al., _Comparing approximate methods for mock catalogues and covariance matrices - I. Correlation function_ , _Mon. Not. R. Astron. Soc._ 482 (2019) 1786 [1806.09477].
  * [47] V. Springel, _The cosmological simulation code GADGET-2_ , _Mon. Not. R. Astron. Soc._ 364 (2005) 1105 [astro-ph/0505010].
  * [48] O. H. E. Philcox and D. J. Eisenstein, _Estimating covariance matrices for two- and three-point correlation function moments in Arbitrary Survey Geometries_ , _Mon. Not. R. Astron. Soc._ 490 (2019) 5931 [1910.04764].
  * [49] Z. Slepian, D. J. Eisenstein, F. Beutler, C.-H. Chuang, A. J. Cuesta, J. Ge et al., _The large-scale three-point correlation function of the SDSS BOSS DR12 CMASS galaxies_ , _Mon. Not. R. Astron. Soc._ 468 (2017) 1070 [1512.02231].
  * [50] J. N. Fry and E. Gaztañaga, _Biasing and hierarchical statistics in large-scale structure_ , _Astrophys. J._ 413 (1993) 447 [astro-ph/9302009].
  * [51] K. C. Chan, R. Scoccimarro and R. K. Sheth, _Gravity and large-scale nonlocal bias_ , _Phys. Rev. D_ 85 (2012) 083509 [1201.3614].
  * [52] T. Baldauf, U. Seljak, V. Desjacques and P. McDonald, _Evidence for quadratic tidal tensor bias from the halo bispectrum_ , _Phys. Rev. D_ 86 (2012) 083540 [1201.4827].
  * [53] A. Lewis, A. Challinor and A. Lasenby, _Efficient Computation of Cosmic Microwave Background Anisotropies in Closed Friedmann-Robertson-Walker Models_ , _Astrophys. J._ 538 (2000) 473 [astro-ph/9911177].
  * [54] J. Hartlap, P. Simon and P. Schneider, _Why your model parameter confidences might be too optimistic. Unbiased estimation of the inverse covariance matrix_ , _Astron. Astrophys._ 464 (2007) 399 [astro-ph/0608064].
  * [55] E. Sellentin and A. F. Heavens, _Parameter inference with estimated covariance matrices_ , _Mon. Not. R. Astron. Soc._ 456 (2016) L132 [1511.05969].
  * [56] W. J. Percival, A. J. Ross, A. G. Sánchez, L. Samushia, A. Burden, R. Crittenden et al., _The clustering of Galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: including covariance matrix errors_ , _Mon. Not. R. Astron. Soc._ 439 (2014) 2531 [1312.4841].
  * [57] D. Foreman-Mackey, D. W. Hogg, D. Lang and J. Goodman, _emcee: The MCMC Hammer_ , _Publications of the Astronomical Society of the Pacific_ 125 (2013) 306 [1202.3665].
  * [58] D. J. Paz and A. G. Sánchez, _Improving the precision matrix for precision cosmology_ , _Mon. Not. R. Astron. Soc._ 454 (2015) 4326 [1508.03162].
  * [59] T. Lazeyras, C. Wagner, T. Baldauf and F. Schmidt, _Precision measurement of the local bias of dark matter halos_ , _Journal of Cosmology and Astro-Particle Physics_ 2016 (2016) 018 [1511.01096].
  * [60] T. Lazeyras and F. Schmidt, _Beyond LIMD bias: a measurement of the complete set of third-order halo bias parameters_ , _Journal of Cosmology and Astro-Particle Physics_ 2018 (2018) 008 [1712.07531].
  * [61] M. M. Abidi and T. Baldauf, _Cubic halo bias in Eulerian and Lagrangian space_ , _Journal of Cosmology and Astro-Particle Physics_ 2018 (2018) 029 [1802.07622].
  * [62] J. Byun, A. Oddo, C. Porciani and E. Sefusatti, _Towards cosmological constraints from the compressed modal bispectrum: a robust comparison of real-space bispectrum estimators_ , _Journal of Cosmology and Astro-Particle Physics_ 2021 (2021) 105 [2010.09579].
  * [63] A. Eggemeier, R. Scoccimarro and R. E. Smith, _Bias loop corrections to the galaxy bispectrum_ , _Phys. Rev. D_ 99 (2019) 123514 [1812.03208].
  * [64] N. S. Sugiyama, S. Saito, F. Beutler and H.-J. Seo, _A complete FFT-based decomposition formalism for the redshift-space bispectrum_ , _Mon. Not. R. Astron. Soc._ 484 (2019) 364 [1803.02132].
  * [65] D. Gualdi, H. Gil-Marín, M. Manera, B. Joachimi and O. Lahav, _Geometrical compression: a new method to enhance the BOSS galaxy bispectrum monopole constraints_ , _Mon. Not. R. Astron. Soc._ (2019) L1 [1901.00987].
  * [66] F. Bernardeau, S. Colombi, E. Gaztañaga and R. Scoccimarro, _Large-scale structure of the Universe and cosmological perturbation theory_ , _Phys. Rep._ 367 (2002) 1 [astro-ph/0112551].
  * [67] L. Blot, M. Crocce, E. Sefusatti, M. Lippich, A. G. Sánchez, M. Colavincenzo et al., _Comparing approximate methods for mock catalogues and covariance matrices II: power spectrum multipoles_ , _Mon. Not. R. Astron. Soc._ 485 (2019) 2806 [1806.09497].
  * [68] E. Sefusatti, M. Crocce, R. Scoccimarro and H. M. P. Couchman, _Accurate estimators of correlation functions in Fourier space_ , _Mon. Not. R. Astron. Soc._ 460 (2016) 3624 [1512.07295].
  * [69] V. Yankelevich, I. G. McCarthy, J. Kwan, S. G. Stafford and J. Liu, _The halo bispectrum as a sensitive probe of massive neutrinos and baryon physics_ , _arXiv e-prints_ (2022) arXiv:2202.07680 [2202.07680].
  * [70] S. Avila, S. G. Murray, A. Knebe, C. Power, A. S. G. Robotham and J. Garcia-Bellido, _HALOGEN: a tool for fast generation of mock halo catalogues_ , _Mon. Not. R. Astron. Soc._ 450 (2015) 1856 [1412.5228].
  * [71] A. Izard, M. Crocce and P. Fosalba, _ICE-COLA: towards fast and accurate synthetic galaxy catalogues optimizing a quasi-N-body method_ , _Mon. Not. R. Astron. Soc._ 459 (2016) 2327 [1509.04685].
  * [72] F. S. Kitaura, G. Yepes and F. Prada, _Modelling baryon acoustic oscillations with perturbation theory and stochastic halo biasing._ , _Mon. Not. R. Astron. Soc._ 439 (2014) L21 [1307.3285].
  * [73] G. Stein, M. A. Alvarez and J. R. Bond, _The mass-Peak Patch algorithm for fast generation of deep all-sky dark matter halo catalogues and its N-body validation_ , _Mon. Not. R. Astron. Soc._ 483 (2019) 2236 [1810.07727].
  * [74] V. Desjacques, D. Jeong and F. Schmidt, _Large-scale galaxy bias_ , _Phys. Rep._ 733 (2018) 1 [1611.09787].