# Cross-Correlated Force Measurement for Thermal Noise Reduction in Torsion
Pendulum

Yusuke Okuma<EMAIL_ADDRESS>Department of Physics,
University of Tokyo, Bunkyo, Tokyo 113-0033, Japan  Kiwamu Izumi
<EMAIL_ADDRESS>Department of Space Astronomy and Astrophysics,
Institute of Space and Astronautical Science, JAXA  Kentaro Komori  Research
Center for the Early Universe (RESCEU), Graduate School of Science, University
of Tokyo, Bunkyo, Tokyo 113-0033, Japan  Department of Physics, University of
Tokyo, Bunkyo, Tokyo 113-0033, Japan  Masaki Ando  Research Center for the
Early Universe (RESCEU), Graduate School of Science, University of Tokyo,
Bunkyo, Tokyo 113-0033, Japan  Department of Physics, University of Tokyo,
Bunkyo, Tokyo 113-0033, Japan

###### Abstract

The torsion pendulum is a prevailing instrument for measuring small forces
acting on a solid body or those between solid bodies. While it offers powerful
advantages, the measurement precision suffers from thermal noises of the
suspending wires giving rise to stochastic torque noises. This letter proposes
a new scheme to reduce the effect of such noise by employing a double torsion
pendulum and cross-correlation technique based on the theoretical analysis
that the thermal torque noise appears at each end of the suspending wire
differentially. Cross-correlating two synthesized data streams which are
composed of the rotation angles of two torsion stages, it yields the power
spectral density estimate of external forces acting on the lower stage with
the reduced effect from the thermal torque noises. As an example use case, we
discuss the application to the study on the coupling strength of ultra light
dark matter to standard model particles. Our evaluation indicates that the
upper limit may be improved by an order of magnitude than the previous
experiments at 2 mHz, which corresponds to about $8\times
10^{-18}~{}\mathrm{eV}$.

##  I Introduction

Precision measurement of small forces acting on a solid body has played an
indispensable role in exploring modern physics. In these measurements, the
torsion pendulum has been often employed as a core instrument. Since the
epoch-making measurements of the Coulomb force and gravitational constant [1]
in the late 18th century, the torsion pendulum has remained as a powerful
experimental tool. Nowadays, it has been successfully utilized in multiple
disciplines such as Casimir force measurements [2], tests of the equivalence
principle [3], verification of gravitational inverse square law [4] and test
of gravitational coupling between small masses [5]. Moreover, using optical
spring have made torsion pendulum technique more sensitive sensor [6, 7].

The torsion pendulum offers advantages critical for the precision force
measurement, including the provision of a low mechanical resonant frequency
and decoupling of the rotational mode from the gravitational pull from the
Earth [8]. However, it has been known that the measurement precision for
external forces is limited by thermal noise from the suspending wire(s) [9,
10]. A straightforward workaround is to cryogenically cool down the apparatus
to reduce the magnitude of thermal noise [11]. However, in practice, it could
introduce a set of complex systems and therefore the design must be carried
out with a great care. Furthermore, the setup must not lose the desired
surrounding environment relevant to the application.

In this letter, we propose a new measurement concept for reducing the effect
of thermal noise in the torsion pendulum while maintaining the ability to
measure small torques or equivalently forces acting on the torsion pendulum.
The core idea is to synthesize two measurement signal streams such that the
cross correlated spectrum between the two provides an estimate of power
spectral density (PSD) of the forces where the contamination from thermal
noise is integrated away due to the lack of coherence. We determined that the
scheme can be implemented in a double torsion setup as opposed to the single
one. This scheme can be implemented together with a cryogenic system and other
noise attenuation systems.

Albeit the drawback of losing the phase information associated with the forces
being measured, it offers a way to probe the magnitude of the forces in the
form of PSD. A number of experiments may benefit from this scheme, including
test of gravitational inverse square law [12], gravitational coupling between
small masses [5], dark matter search [13] and the ground verification of the
gravitational reference sensors for space gravitational wave antennas [14, 15,
16, 17].

Additionally, this letter discusses the application to the study on the
coupling strength of ultra light dark matter to standard model particles,
$g_{B-L}$[13], as an example. We show that our scheme can improve the
sensitivity level by an order of magnitude than the previous record at 2 mHz.

##  II Setup and definitions

Let us assume the torques or forces being measured to be stationary
stochastic. The quantity we desire to estimate in the measurement is the
torques in an amplitude spectral density in Nm/$\sqrt{\textrm{Hz}}$. The
torque can be subsequently converted into the effective force acting on a mass
by taking the lever arm into account.

As discussed later, thermal noise from the suspending wire would not be
reduced if a single-stage torsion pendulum was employed. For the reason, we
begin the discussion with a two-stage torsion pendulum as shown in Fig. 1. The
lower wire suspends the first torsion mass. The upper wire and second torsion
mass suspends the lower wire underneath. Each stage rotates about the wire
horizontally with natural frequencies at some mHz. For simplicity, the
discussion hereafter is focused on these rotational degrees of freedom only.

The upper and lower stages receive external torques as well as thermal torque
noises from the wires. The rotational angles of the two stages, $\phi_{1}$ and
$\phi_{2}$, are continuously measured and recorded.

Figure 1: Illustration of two stage torsion pendulum.

The equations of motion for the two stages in the rotational degrees of
freedom can be expressed by

$\displaystyle\begin{pmatrix}I_{1}\ddot{\phi_{1}}\\\
I_{2}\ddot{\phi_{2}}\end{pmatrix}=\begin{pmatrix}-I_{1}\omega_{1}^{2}&I_{1}\omega_{1}^{2}\\\
I_{1}\omega_{1}^{2}&-I_{2}\omega_{2}^{2}-I_{1}\omega_{1}^{2}\end{pmatrix}\begin{pmatrix}\phi_{1}\\\
\phi_{2}\end{pmatrix}+\begin{pmatrix}n_{1}\\\ n_{2}\end{pmatrix},$ (1)

where $I_{j}$, $\omega_{j}$ and $n_{j}$ are the moments of inertia, resonant
frequencies and torques applied to the $j$-th stage with $j\in 1,2$. In our
convention, $j=1$ corresponds to the lower stage while $j=2$ for the upper
(see Fig. 1). To reduce clutter, damping effects are implicitly incorporated
in the resonant frequencies, such that $\omega_{j}^{2}$ can be later replaced
in the frequency domain by $\omega_{j}^{2}\left(1+iQ_{j}^{-1}\right)$ with
$Q_{j}$ the $Q$ value.

Fourier transforming Eq. (1) , we obtain

$\displaystyle\begin{pmatrix}\Phi_{1}\\\
\Phi_{2}\end{pmatrix}=\hat{M}\begin{pmatrix}N_{1}\\\ N_{2}\end{pmatrix},$ (2)

where $N_{j}\left(\omega\right)$ and $\Phi_{j}\left(\omega\right)$ are the
Fourier domain representations of the torques $n_{j}(t)$ and rotational angles
$\phi_{j}(t)$, respectively. $\hat{M}\left(\omega\right)$ is a transfer matrix
of the system defined by

$\displaystyle\hat{M}$ $\displaystyle=$ $\displaystyle
a\begin{pmatrix}-I_{2}(\omega_{2}^{2}-\omega^{2})+I_{1}\omega_{1}^{2}&I_{1}\omega_{1}^{2}\\\
I_{1}\omega_{1}^{2}&I_{1}(\omega_{1}^{2}-\omega^{2})\end{pmatrix},$ (3)
$\displaystyle a$ $\displaystyle=$
$\displaystyle\frac{1}{I_{1}I_{2}\omega^{4}-I_{1}(I_{1}\omega_{1}^{2}+I_{2}\omega_{2}^{2}+I_{2}\omega_{1}^{2})\omega^{2}+I_{1}I_{2}\omega_{1}^{2}\omega_{2}^{2}}.$

It is obvious that the rotational motions are consequences of the torques
which are linearly combined.

Now, even if no external forces exist, the stages rotate in stochastic manner
due to the presence of torques generated by thermal noise of the two wires.
Such thermal torque noise in the lower stage can be calculated by directly
applying the fluctuation dissipation theorem [9] as,

$\displaystyle S^{N}_{\rm th}=4k_{B}T\Re\left\\{\frac{1}{i\omega
M_{11}}\right\\}\quad[\rm N^{2}m^{2}/Hz],$ (4)

where $M_{11}$ is the (1,1) element of the transfer matrix $\hat{M}$, $k_{B}$
is the Boltzmann constant and $T$ the temperature.

An important remark is that Eq. (4) is consistent with the physical process
that thermal torque noise of the lower wire drives both upper and lower stages
in differential fashion. For this reason and later convenience, we redefine
the torque vector in Eq. (2) as,

$\vec{N}=N_{1}\begin{pmatrix}1\\\ 0\end{pmatrix}+N_{\rm
th2}\begin{pmatrix}0\\\ 1\end{pmatrix}+N_{\rm th1}\begin{pmatrix}1\\\
-1\end{pmatrix},$ (5)

where we have introduced three independent and uncorrelated torques, $N_{1}$
and $N_{\rm thj}$ ($j\in 1,2$). $N_{1}$ represents an external torque acting
on the lower stage with thermal noise torques removed. This is the quantity we
desire to measure. On the other hand, $N_{\rm th1}$ and $N_{\rm th2}$
represent thermal torque noises from the lower and upper wires, respectively.

One can prove that $N_{\textrm{th1}}$ indeed comes into the system
differentially on the upper and lower stages by propagating the PSD of thermal
torque noises $S^{N}_{\rm thj}=4k_{B}TI_{j}\omega_{j}^{2}/\omega Q_{j}$ via
the equation of motion to the net torques as seen in the lower stage. It
results in the identical expression to Eq. (4).

##  III Cross-correlated force measurement

Plugging torque vector (5) into Eq. (2), one can explicitly write down the
rotational angles as,

$\displaystyle\Phi_{1}$ $\displaystyle=$ $\displaystyle
M_{11}N_{1}+M_{12}N_{\rm th2}+(M_{11}-M_{12})N_{\rm th1},$ (6)
$\displaystyle\Phi_{2}$ $\displaystyle=$ $\displaystyle
M_{21}N_{1}+M_{22}N_{\rm th2}+(M_{21}-M_{22})N_{\rm th1}.$ (7)

There are three unknown quantities i.e., $N_{1}$ and$~{}N_{\rm thj}$, while
the number of equations is only two. Therefore, given the measurement of
angles $\Phi_{j}$, one cannot solve the equation for $N_{1}$ without a
contamination from $N_{\textrm{thj}}$.

Nonetheless, we discovered a generic solution to estimate the PSD of $N_{1}$,
namely $S^{N}_{11}$, independently of the other two. Linearly combining Eqs
(6) and (7), one can obtain two synthesized data streams as,

$\displaystyle N_{1}+N_{\rm th2}$ $\displaystyle=$
$\displaystyle\frac{(M_{22}-M_{21})\Phi_{1}-(M_{12}-M_{11})\Phi_{2}}{\textrm{det}\hat{M}},$
(8) $\displaystyle N_{1}+N_{\rm th1}$ $\displaystyle=$
$\displaystyle\frac{M_{22}\Phi_{1}-M_{12}\Phi_{2}}{\textrm{det}\hat{M}},$ (9)

where $\textrm{det}\hat{M}$ is the determinant of matrix $\hat{M}$.
Subsequently, taking the cross correlated spectrum of two synthesized data
sets, one obtains,

$\displaystyle S^{N}_{11}$ $\displaystyle=$
$\displaystyle\frac{1}{|\textrm{det}\hat{M}|^{2}}\\{M_{22}(M_{22}-M_{21})^{*}S^{\Phi}_{11}$
(10) $\displaystyle-
M_{22}(M_{12}-M{11})^{*}S^{\Phi}_{12}-M_{12}(M_{22}-M_{21})^{*}S^{\Phi}_{21}$
$\displaystyle+M_{12}(M_{12}-M{11})^{*}S^{\Phi}_{22}\\},$

where $N_{1}$ and $N_{\rm thj}$ are assumed to be independent of each other
i.e., $S^{N}_{1\rm thj}=S^{N}_{\rm th1th2}=0$. With measurements for a
sufficiently long duration of time, one can obtain $S^{N}_{11}$ separated from
the contributions of the thermal noises. We note that since an external torque
acting on the upper stage come into the system in the manner identical to
$N_{\rm th2}$, it is also rejected as long as it is incoherent to the others.

Plugging Eq. (3) into (10), one can explicitly write down the PSD of torque
acting on the lower stage as,

$\displaystyle S^{N}_{11}$ $\displaystyle=$ $\displaystyle
I_{1}^{2}\omega^{2}(\omega^{2}-\omega_{1}^{*2})S^{\Phi}_{11}+I_{1}^{2}\omega^{2}\omega_{1}^{*2}S^{\Phi}_{12}$
(11)
$\displaystyle+I_{1}I_{2}(\omega^{2}-\omega_{2}^{2})(\omega^{2}-\omega_{1}^{*2})S^{\Phi}_{21}$
$\displaystyle+I_{1}I_{2}(\omega^{2}-\omega_{2}^{2})\omega_{1}^{*2}S^{\Phi}_{22}.$

This serves as a recipe for obtaining the cross-correlated force estimate. The
rotation angles $\phi_{j}$ must be continuously monitored and recorded to
obtain the auto- and cross-correlated spectra, $S^{\Phi}_{jk}$. Note that, in
real experiment, one must take the real part of Eq. (11) to obtain the PSD.
Taking the absolute value of Eq. (11) is fine for evaluating only
$S_{11}^{N}$. But considering residual thermal noise (see following sections),
taking the real part results in smaller noise. Besides, the mechanical
properties must be characterized to obtain the handful parameters, namely
$I_{j}$ and $\omega_{j}$. An implicit assumption is that the mechanical
properties are time-invariant throughout the measurements.

We remark that the single-stage torsion pendulum cannot adopt this cross-
correlated force estimate. In the case of a single-stage torsion pendulum, it
is equivalent to set $I_{2}\rightarrow\infty$, resulting in
det$\hat{M}\rightarrow 0$ in Eq. (3). This means that $N_{1}$ and
$N_{\textrm{th1}}$ are degenerate unless another information is incorporated.

One drawback in this scheme is that the phase information becomes
inaccessible. The resulting estimation of external torque $n_{1}$ is presented
in the form of PSD. Therefore, it is not possible to bring this quantity back
into a set of time series data representing $n_{1}$.

##  IV Practical limitations

The discussion so far has been restricted to an ideal situation. However, in
practice, the performance of the force estimation would be deteriorated in
three aspects, namely an estimation bias, residual thermal noise torques and
sensing noises. While the first two effects arise due to inaccurate knowledge
on the system parameters, the last adds another noise into the estimation
independently of the accuracy for the system parameters.

Figure 2: (Left) Contribution of residual thermal torque and sensing noises to
the torque sensitivity of lower stage. The colored contours represent the
occurrence intervals of 99.7%, 95% and 68% at given frequency points. (a) Red
line represents the most frequent value. (b) Dashed line shows thermal noise
without cross-correlated technique using Eq (4) for comparison. (c) Dash-
dotted line shows sensing noise. (Right) Relative bias in the estimate of
$S^{N}_{11}$. Filled contours obey the same color scheme as the left panel
based on the occurrence. (d) Red line represents the most frequent value.

As written in Eq. (11), the system parameters $I_{j}$, $\omega_{j}$ and
$Q_{j}$ have to be known accurately. Any errors in the system parameter values
would lead to an incomplete result. In fact, the resulting PSD estimate would
underestimate or overestimate the true value $S^{N}_{11}$ and simultaneously
leave residual thermal noise torques. When the system parameters are
inaccurate, the right hand side in Eq. (11) can be expressed by,

$\displaystyle({\rm RHS(\ref{eq:sol})})=S^{N}_{11}+\Delta S^{N}_{11}+\Delta
S^{N}_{\rm th1}+\Delta S^{N}_{\rm th2},$ (12)

where $\Delta X$ denotes the error associated with quantity $X$. While the
first term represents the true value, the second term acts as an estimation
bias against the true value. The last two terms are residual thermal torque
noises which would mask the true value in the spectrum.

Table 1 summarizes the contributions of various errors to the estimation. The
errors were expanded up to the first order and propagated to the resulting
estimate in two distinct frequency regimes, namely $\omega\ll\omega_{j}$ and
$\omega_{j}\ll\omega$. As seen in the table, the relative signs between the
errors matter.

To quantitatively highlight such effects, we ran a numerical simulation with a
concrete set of system parameters as $T=300$ K, $I_{j}=3\times 10^{-3}$ kg m2,
$\omega_{j}/(2\pi)=10^{-3}$ Hz and $Q=10$. In the simulation, the system
parameters were deliberately changed by $\alpha+\delta\alpha$ ($\alpha\in
I_{j},\,\omega_{j},\,Q_{j}$), with $\delta\alpha$ an error with respect to the
true vale $\alpha$. We adopted the Monte Carlo analysis [18], where we assumed
$\delta\alpha$ to obey a gaussian distribution with zero-mean and fractional
standard deviation of 0.01. In total, 100,000 realizations were computed. The
system parameters were randomly sampled at each time.

The left panel in Fig. 2 shows the effect from the residual thermal torque
noises in the PSD. As summarized in Table 1, the dominant contributors at
frequencies higher than the resonances were identified to be errors in the
moment of inertia or $\delta I_{1,2}/I_{1,2}$. This is because the mechanical
responses of the system resemble a free mass response in this regime. In
contrast, residual thermal torque noise was confirmed to scale as $\omega^{2}$
in PSD at frequencies lower than the resonances. Therefore, thermal noise is
further suppressed at low frequencies.

Table 1: Summary of various contributions to the resulting torque estimate. | Low freq., $\omega\ll\omega_{1,2}$ | High freq., $\rm\omega_{1,2}\ll\omega$
---|---|---
Residual thermal | $\displaystyle\frac{\omega^{2}}{\omega_{1}^{*2}}\left(-\frac{\delta I_{1}}{I_{1}}+\frac{\delta I_{2}}{I_{2}}+2\frac{\delta\omega_{2}^{*}}{\omega_{2}^{*}}\right)S^{N}_{\rm th1}+\frac{2I_{1}\omega^{2}}{I_{2}\omega_{2}^{2}}\frac{\delta\omega_{1}}{\omega_{1}}S^{N}_{\rm th2}$ | $\displaystyle\left(\frac{\delta I_{1}}{I_{1}}-\frac{\delta I_{2}}{I_{2}}\right)S^{N}_{\rm th1}$
Relative uncertainty | $\displaystyle\left(\frac{\delta I_{1}}{I_{1}}+\frac{\delta I_{2}}{I_{2}}+\frac{2\delta\omega_{1}}{\omega_{1}}+\frac{2\delta\omega^{*}_{2}}{\omega^{*}_{2}}\right)S^{N}_{11}$ | $\displaystyle\frac{2\delta I_{1}}{I_{1}}S^{N}_{11}$
Sensing noise | $-I_{1}I_{2}\omega_{1}^{*2}\omega_{2}^{2}S^{\Phi}_{n2}$ | $I_{1}^{2}\omega^{4}S^{\Phi}_{n1}$

The right panel in Fig. 2 shows the simulation results for the estimation bias
against the true value in $S^{N}_{1,1}$. The relative bias normalized by the
true value is almost flat across the frequency band with amplification
structures at around the resonant frequencies.

Finally, the effect of sensing noise is studied. We add moderate sensing noise
of $\sqrt{S^{\Phi}_{nj}}=10^{-8}$ rad/$\sqrt{\textrm{Hz}}$ to the measurements
of the rotational degrees of freedom. As shown in Table 1, the lower moment of
inertia and the lower resonant frequency, the smaller level sensing noise
limit lie. This is because a small moment of inertia and resonant frequencies
lead to mechanically soft system which gives a large response in rotation
angles for given torques. At frequencies above $\omega_{j}$, sensing noise
reaches the level as high as those without the cross-correlation scheme. In
this region, the sensing noise scales as $\omega^{4}$ in PSD.

The system parameters may be fine-tuned offline to minimize the residual
thermal torque noises and bias in the estimation. Even if such an optimization
is implemented, sensing noise would persist as a major source limiting the
torque sensitivity. Nonetheless, sensing noise can be also reduced if a pair
of independent angular sensors are introduced for each stage. In fact, such a
scheme has been used in the development of the gravitational reference sensor
[19]. Therefore both thermal and sensing noises would be eliminated with the
measurement for infinitely long time. In reality, the measurement time of
approximately a day is required for the cross-correlated spectrum to reach the
level limited by the inaccuracy of system parameters as shown in Fig(2) above
1 mHz.

##  V Discussion

Figure 3: Expected upper limits on the coupling strength of ULDM to standard
model particle. (a) Upper limits on the amplitude of field along X axis in
geocentric celestial coordinates [13]. (b) Upper limit by MICROSCOPE [20, 21].
(c) Upper limit by Eöt-Wash equivalence principle experiment [22]. (d) Thermal
noise without cross-correlation. (e) Sum of the most frequent value and
sensing noise.

One of the experiments that may benefit from the cross-correlated force
measurement is the dark matter search using a torsion pendulum. Ultra light
dark matter (ULDM), a candidate for dark matter, has the possibility to be
composed of vector field [23] and could couple to B-L with B and L being
baryon and lepton numbers, respectively. Using the torsion pendulum, one can
experimentally study the coupling strength of ULDM to standard model
particles, $g_{\rm B-L}$. The Eöt-Wash equivalence principle experiment [22]
limits $g_{\rm B-L}\leq 1\times 10^{-24}$ below tens of mHz. At 2 mHz, the
upper limit for $g_{\rm B-L}$ is set to be less than $1\times 10^{-25}$ [13],
which is the most stringent limit to our best knowledge.

The cross-correlated force measurement is capable of improving these limits
below a few mHz including the one at 2 mHz by an order of magnitude as shown
in Fig. 3. Our expectation is derived based on [24]. We ran another set of
Monte-Carlo simulation using the system parameters $T=300$ K,
$I_{j}=3.78\times 10^{-5}$ kg m2, $\omega_{j}/(2\pi)=0.5\times 10^{-3}$ Hz,
$Q_{j}=200$, $l=0.05$ m, $m=19.4\times 10^{-3}$ kg, $\tau=91.4$ days and
$\Delta_{\rm B-L}$ = 0.0359, where $l$ is the arm length of torsion pendulum,
$m$ is the mass of test bodies, $\tau$ is the data measurement time and
$\Delta_{\rm B-L}$ the difference in the charge to-mass ratios of the two
atomic species in the dipole. Also in this simulation, $I_{j},~{}\omega_{j}$
and $Q_{j}$ are assumed to obey a gaussian distribution with zero-mean and
fractional standard deviation of 0.01. We chose small $Q$ values comparing to
[13] in which $Q=460,000$ is achieved. This is because a small $Q$ value makes
it easy to commission the apparatus in practice.

##  VI Conclusion

This letter proposes a new scheme, the cross-correlated force measurement. The
scheme utilizes the cross-correlation to overcome the sensitivity limit of
thermal torque noises which would otherwise limit the sensitivity of precision
force measurement in the torsion pendulum. Both qualitative and quantitative
analyses for the sensitivity are presented. Inaccurate knowledge on the system
parameters would limit the practical sensitivity, which invites a concrete
idea for offline optimization scheme to further reduce residual thermal torque
noises. The scheme can be implemented with the sensing noise reduction method
[19], enabling an apparatus to be significantly less sensitive to both thermal
and sensing noises. As a use case, the application to the Ultra-light dark
matter search is discussed. The cross-correlated force measurement is capable
of improving the upper limits of $g_{\rm B-L}$ at frequencies below tens of
mHz. The scheme, however, shall be experimentally verified.

## References

  * [1] H. Cavendish, Philosophical Transactions of the Royal Society of London, 469 (1798)
  * [2] S. K. Lamoreaux, Physical Review Letters 78, 5 (1997)
  * [3] E. G. Adelberger, C. Stubbs, B. R. Heckel, Y. Su, H. Swanson, G. Smith, J. Gundlach, and W. Rogers, Physical Review D 42, 3267 (1990)
  * [4] J. Hoskins, R. Newman, R. Spero, and J. Schultz, Physical Review D 32, 3084 (1985)
  * [5] T. Westphal, H. Hepach, J. Pfaff, and M. Aspelmeyer, Nature 591, 225 (2021)
  * [6] Q.-L. Wang, H.-C. Yeh, Z.-B. Zhou, and J. Luo, Physical Review A 80, 043811 (2009)
  * [7] K. Komori, Y. Enomoto, C. P. Ooi, Y. Miyazaki, N. Matsumoto, V. Sudhir, Y. Michimura, and M. Ando, Physical Review A 101, 011802 (2020)
  * [8] G. Gillies and R. Ritter, Review of scientific instruments 64, 283 (1993)
  * [9] P. R. Saulson, Physical Review D 42, 2437 (1990)
  * [10] G. I. González and P. R. Saulson Physics Letters A 201, 12 (1995)
  * [11] R. Newman, M. Bantel, E. Berg, and W. Cross, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20140025 (2014)
  * [12] W.-H. Tan, S.-Q. Yang, C.-G. Shao, J. Li, A.-B. Du, B.-F. Zhan, Q.-L. Wang, P.-S. Luo, L.-C. Tu, and J. Luo, Physical Review Letters 116, 131101 (2016)
  * [13] E. Shaw, M. Ross, C. Hagedorn, E. Adelberger, and J. Gundlach, Physical Review D 105, 042007 (2022)
  * [14] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, et al., arXiv preprint arXiv:1702.00786 (2017)
  * [15] J. Luo, L.-S. Chen, H.-Z. Duan, Y.-G. Gong, S. Hu, J. Ji, Q. Liu, J. Mei, V. Milyukov, M. Sazhin, et al., Classical and Quantum Gravity 33, 035010 (2016)
  * [16] W.-R. Hu and Y.-L. Wu, The taiji program in space for gravitational wave physics and the nature of gravity (2017)
  * [17] S. Kawamura, M. Ando, N. Seto, S. Sato, T. Nakamura, K. Tsubono, N. Kanda, T. Tanaka, J. Yokoyama, I. Funaki, et al., Classical and Quantum Gravity 28, 094011 (2011)
  * [18] N. Metropolis and S. Ulam, Journal of the American statistical association 44, 335 (1949)
  * [19] L. Carbone, G. Ciani, R. Dolesi, M. Hueller, D. Tombolato, S. Vitale, W. J. Weber, and A. Cavalleri, Physical Review D 75, 042001 (2007)
  * [20] P. Touboul, G. Métris, M. Rodrigues, Y. André, Q. Baghi, J. Bergé, D. Boulanger, S. Bremer, P. Carle,R. Chhun, et al., Physical Review Letters 119, 231101 (2017)
  * [21] P. Fayet, Physical Review D 97, 055039 (2018)
  * [22] T. A. Wagner, S. Schlamminger, J. Gundlach, and E. G. Adelberger, Classical and Quantum Gravity 29, 184002 (2012)
  * [23] A. E. Nelson and J. Scholtz, Physical Review D 84, 103501 (2011)
  * [24] J. Manley, M. D. Chowdhury, D. Grin, S. Singh, and D. J. Wilson, Physical Review Letters 126, 061301 (2021)


# The monopole and quadrupole moments of the Epoch of Reionization (EoR) 21-cm
bispectrum

Sukhdeep Singh Gill1, Suman Pramanick1, Somnath Bharadwaj1, Abinash Kumar
Shaw2 and Suman Majumdar3,4
1Department of Physics & Centre for Theoretical Studies, Indian Institute of
Technology Kharagpur, Kharagpur 721302, India.
2Astrophysics Research Center of the Open University (ARCO), The Open
University of Israel, 1 University Road, Ra’anana 4353701, Israel.
3Department of Astronomy, Astrophysics and Space Engineering, Indian Institute
of Technology, Indore, India
4Department of Physics, Blackett Laboratory, Imperial College, London SW7 2AZ,
U.K E-mail<EMAIL_ADDRESS>
###### Abstract

We study the monopole ($\bar{B}^{0}_{0}$) and quadrupole ($\bar{B}^{0}_{2}$)
moments of the 21-cm bispectrum (BS) from EoR simulations and present results
for squeezed and stretched triangles. Both $\bar{B}^{0}_{0}$ and
$\bar{B}^{0}_{2}$ are positive at the early stage of EoR where the mean
neutral hydrogen (H i) density fraction $\bar{x}_{\rm\text{H\,{i}}}\approx
0.99$. The subsequent evolution of $\bar{B}^{0}_{0}$ and $\bar{B}^{0}_{2}$ at
large and intermediate scales $(k=0.29$ and $0.56\,{\rm Mpc}^{-1}$
respectively) is punctuated by two sign changes which mark transitions in the
H i distribution. The first sign flip where $\bar{B}^{0}_{0}$ becomes negative
occurs in the intermediate stages of EoR $(\bar{x}_{\rm\text{H\,{i}}}>0.5)$,
at large scale first followed by the intermediate scale. This marks the
emergence of distinct ionized bubbles in the neutral background.
$\bar{B}^{0}_{2}$ is relatively less affected by this transition, and it
mostly remains positive even when $\bar{B}^{0}_{0}$ becomes negative. The
second sign flip, which affects both $\bar{B}^{0}_{0}$ and $\bar{B}^{0}_{2}$,
occurs at the late stage of EoR $(\bar{x}_{\rm\text{H\,{i}}}<0.5)$. This marks
a transition in the topology of the H i distribution, after which we have
distinct H i islands in an ionized background. This causes $\bar{B}^{0}_{0}$
to become positive. The negative $\bar{B}^{0}_{2}$ is a definite indication
that the H i islands survive only in under-dense regions.

###### keywords:

cosmology: reionization, diffuse radiation –– simulations - methods:
statistical

††pubyear: 2023††pagerange: The monopole and quadrupole moments of the Epoch
of Reionization (EoR) 21-cm bispectrum–References

##  1 Introduction

The Epoch of Reionization (EoR), when the first luminous objects ionized the
neutral hydrogen (H i ) gas in the intergalactic medium (IGM), is one of the
least understood era in the evolution of our Universe. Observations of the
Cosmic Microwave Background (CMB) (e.g. Planck Collaboration et al. 2016), the
absorption spectra of quasars (e.g. Becker et al. 2001), and the luminosity
function of Lyman-$\alpha$ emitters (e.g. Trenti et al. 2010) all indirectly
probe the EoR, however, these provide rather limited information about the
EoR. Observations of the redshifted H i 21-cm radiation are a very promising
probe of the EoR (e.g. Furlanetto et al. 2006). Despite continued efforts
using several telescopes, a detection of the EoR 21-cm power spectrum (PS) is
still forthcoming. The main impediment for detecting the 21-cm PS comes from
foregrounds, which are $\sim 4-5$ order of magnitude stronger (Ali et al.,
2008). At present, the best upper limits on the EoR 21-cm PS come from the
HERA experiment (Abdurashidova et al., 2022). It is anticipated that more
sensitive observations with existing instruments and also the forthcoming
SKA1-Low (Koopmans et al., 2015) will result in a detection.

Much of the effort to quantify the EoR 21-cm signal has focused on the 21-cm
PS. While this completely describes Gaussian fields, an early paper (Bharadwaj
& Pandey, 2005) predicts the EoR 21-cm signal to be highly non-Gaussian and
also provides predictions for the 21-cm bispectrum (BS), which is the lowest
order statistics that captures this non-Gaussianity. The 21-cm BS holds the
potential to offer valuable insights into the evolution of the H i
distribution during the EoR (Yoshiura et al., 2015; Majumdar et al., 2018,
2020; Hutter et al., 2020), tighten the constraints on EoR models (Watkinson
et al., 2022; Tiwari et al., 2022) and probe the IGM physics (Watkinson et
al., 2021; Kamran et al., 2021a, b, 2022).

Figure 1: Evolution of the redshift space EoR 21-cm brightness temperature
maps. This shows slices from a single realization of the simulated signal
during different stages of EoR.

Redshift space distortion (RSD) is an important effect imprinted in the 21-cm
signal (Bharadwaj et al., 2001; Bharadwaj & Ali, 2004), causing an anisotropy
in the 21-cm PS along the line of sight (LoS) direction. Majumdar et al.
(2013) have quantified this anisotropy using the quadrupole moment of the
21-cm PS estimated from semi-numerical simulations of the EoR. Considering the
quadrupole to monopole ratio, they show that at large scales this becomes
negative when the mean H i density fraction $\bar{x}_{\text{H\,{i}}~{}}$ drops
to $\bar{x}_{\text{H\,{i}}~{}}\leq 0.7$. They propose this as a direct
signature of the inside-out reionization scenario. Analytical predictions (eg.
Scoccimarro et al. 1999; Bharadwaj et al. 2020) suggest that the anisotropy in
the BS can similarly be quantified by decomposing it into angular multipoles
using spherical harmonics. In a recent study with the simulated signal,
Majumdar et al. (2020) have shown that RSD can change the amplitude of the
monopole 21-cm BS by $50-200\%$, and can also cause it to flip sign. To the
best of our knowledge, previous studies of the 21-cm BS have all been
restricted to the monopole moment. In this paper, we present predictions for
the quadrupole moment of the EoR 21-cm BS using a large set of statistically
independent simulated signal realizations. As the quadrupole moment is
entirely due to the line of sight effects, it can potentially measure and
quantify the impact of these effects on the non-Gaussian nature of the H i
distribution. Here, we study the quadrupole moment due to RSD. Moreover, the
quadrupole moment of the BS used in conjunction with the monopole can impose
tighter constraints on the EoR model (and its parameters) in comparison to
those achievable using the monopole alone.

A brief outline of the paper follows. Section 2 presents the Methodology,
Section 3 presents the Results, and Section 4 presents the Discussion and
Conclusions. This paper has used the Plank+WP best-fit values of the
cosmological parameters (Planck Collaboration et al., 2020).

##  2 Methodology

We have simulated the EoR 21-cm signal using the semi-numerical code
ReionYuga111https://github.com/rajeshmondal18/ReionYuga (Mondal et al.,
2021a), whose details have also been presented in Majumdar et al. (2013) and
which closely follows the homogeneous reionization scheme of Choudhury et al.
(2009). The dark matter only PM N-body simulation considers a $(215.04~{}{\rm
Mpc})^{3}$ comoving volume with a $3072^{3}$ grids and $1536^{3}$ particles
each of mass $1.09\times 10^{8}\>{\rm M}_{\sun}$. Dark matter halos were
identified as having a minimum mass of $1.09\times 10^{9}{\rm M}_{\sun}$ (10
particles) using FoF algorithm with linking length $0.2$ grid units. The H i
reionization field was calculated on a coarser $384^{3}$ grids of spacing
$0.56\,{\rm Mpc}$ (comoving). The reionization model has three free parameters
namely $(M_{\rm min}$, $N_{\rm ion}$, $R_{\rm mfp})$ for which we have used
the values $(1.09\times 10^{9}~{}{\rm M}_{\sun},23.21,20\,{\rm Mpc})$
respectively. Here $M_{\rm min}$ is the minimum halo mass assumed to host a
reionization source (UV luminous galaxy), $N_{\rm ion}$ is a proportionality
constant that linearly relates the number of ionizing photons emitted by an
ionizing source to the amount of hydrogen contained in the host halo and
$R_{\rm mfp}$ is the mean free path of the ionizing photons. This choice of
simulation parameters corresponds to a model of EoR where reionization starts
around $z=13$ and ends around $z=6$, marking $50\%$ ionization at $z\approx
8$. The integrated Thomson scattering optical depth computed for our fiducial
model is $\tau=0.057$, which is consistent with $\tau=0.054\pm 0.007$ obtained
from CMB measurements (Planck Collaboration et al., 2020). These fiducial
simulations have been used in many previous studies (Mondal et al., 2015;
Mondal et al., 2017a, b; Majumdar et al., 2018, 2020; Shaw et al., 2019;
Pramanick et al., 2023). Varying the model parameters will change the
evolution history of the mean neutral hydrogen fraction
$\bar{x}_{\rm\text{H\,{i}}}$. However, we do not expect the overall topology
of the 21-cm signal at different stages of reionization to change drastically.
The reader is referred to Shaw et al. (2020) for a detailed discussion on the
effect of varying these parameters. We have generated snapshots of the
redshift space 21-cm signal at six different redshifts $z=13,11,10,9,8$ and
$7$. We have run $45$ statistically independent realizations of the
simulations which provide the mean values and the error bars for the results
presented here.

The EoR 21-cm redshift space BS
$B^{s}(\textbf{{k}}_{1},\textbf{{k}}_{2},\textbf{{k}}_{3})$ has been
decomposed into angular multipoles $\bar{B}_{\ell}^{m}(k,\mu,t)$ using
spherical harmonics $Y_{\ell}^{m}(\hat{\textbf{{p}}})$ where

$\bar{B}^{m}_{\ell}(k,\mu,t)=\sqrt{\frac{(2\ell+1)}{4\pi}}\int[Y^{m}_{\ell}(\hat{\textbf{{p}}})]^{*}B^{s}(k,\mu,t,\hat{\textbf{{p}}})\,d\Omega_{\hat{\textbf{{p}}}}\,.$
(1)

With reference to a triangle
$(\textbf{{k}}_{1},\textbf{{k}}_{2},\textbf{{k}}_{3})$, $k$ (the largest side)
quantifies its size, the two dimensionless parameters $(\mu,t)$ with
$0.5\leq\mu,t\leq 1$ and $2\mu t\geq 1$ quantify its shape and the unit vector
$\hat{\textbf{{p}}}$ quantifies its orientation with respect to the LoS
(Bharadwaj et al., 2020). Our entire analysis is restricted to the monopole
$(\ell=0)$ and quadrupole $(\ell=2)$ with $m=0$. Shaw et al. (2021) have
recently proposed a fast estimator to quantify the shape dependence of the
bispectrum monopole. Here, we have slightly modified this and used it to
estimate the binned bispectrum monopole and quadrupole moments for triangles
of all possible unique shapes. The monopole (Majumdar et al., 2020; Mondal et
al., 2021b; Tiwari et al., 2022) and the quadrupole are both found to have the
largest values and the highest statistical significance (signal-to-noise
ratio) for linear triangles, which corresponds to the situation where the two
largest sides of the triangle are nearly aligned in the same direction and the
cosine of the angle between these two sides $\mu\approx 1$. Results are
presented here only for linear triangles of two specific shapes
$(\mu,t)\approx(1.0,1.0)$ and $(1.0,0.5)$ which respectively correspond to
squeezed triangles $(k_{1}\approx k_{2},k_{3}\rightarrow 0)$ and stretched
triangles $(k_{1}/2\approx k_{2}\approx k_{3})$ both of which are illustrated
in Figure 2.

Figure 2:  This shows a squeezed triangle and a stretched triangle, the two
triangle shapes for which the bispectrum results are presented in Figure 3.
Figure 3: The solid lines show the redshift and scale dependence of the
simulated bispectrum multipole moments. We show monopole
$[\Delta^{s}_{0}]^{3}$ (top row) and quadrupole $[\Delta^{s}_{2}]^{3}$ (bottom
row) in terms of mean cubed brightness temperature fluctuations. The left and
right columns correspond to squeezed and stretched triangles, respectively.
The analytical predictions from second-order standard perturbation theory
(2PT) are shown by dashed lines. The yellow-shaded region indicates the range
where the $y$ scale is linear.

##  3 Results

Figure 1 visually depicts the redshift evolution of the 21-cm brightness
temperature distribution. At the early stage ($z=13$), the H i is nearly
completely neutral with $\bar{x}_{\rm\text{H\,{i}}}=0.99$. There are a few
small isolated ionized (H ii ) bubbles that cannot be visually identified.
These isolated H ii bubbles grow and merge during the intermediate stages of
EoR ($z=11,10,9$) where $\bar{x}_{\rm\text{H\,{i}}}<0.5$. The mid-stage
($\bar{x}_{\rm\text{H\,{i}}}\approx 0.5$), which occurs at $z=8$, marks a
transition in the H i distribution. Subsequently, at the late stage $(z=7)$,
the bulk of the H i is ionized ($\bar{x}_{\rm\text{H\,{i}}}=0.15$), and we are
left with isolated H i islands. The emergence of large H ii bubbles during the
intermediate and mid stages and H i islands at the late stage makes the EoR
21-cm signal highly non-Gaussian.

The top row of Figure 3 shows the redshift evolution of
$[\Delta^{s}_{0}]^{3}=k^{6}\bar{B}^{0}_{0}(k,\mu,t)/(2\pi^{2})^{2}$ (redshift
space monopole) which we may interpret as the mean cubed brightness
temperature fluctuation, and the lower row shows $[\Delta^{s}_{2}]^{3}$
(redshift space quadrupole), which has been similarly defined using
$\bar{B}^{0}_{2}(k,\mu,t)$. The left and right columns correspond to squeezed
and stretched triangles, respectively. The different curves in each panel
correspond to $k=0.29,0.56,1.18$ and $2.48\,{\rm Mpc}^{-1}$, which we refer to
as large, intermediate, small, and smallest scales, respectively. We have also
shown (dashed lines) the predictions from second-order perturbation theory
(2PT) using the formulae from Mazumdar et al. (2020) with $P^{r}(k)$ being the
real space 21-cm PS estimated from the simulations as input. The 2PT
predictions assume that the H i traces the underlying matter distribution. Due
to the limited simulation volume, we do not obtain a statistically significant
(SNR $>3$) estimate of $[\Delta^{s}_{2}]^{3}$ for stretched triangles at large
scales, and we have not shown this here. Considering the overall behaviour, we
see that in most cases, the results for squeezed and stretched triangles are
very similar, and we discuss these together.

The early stage of EoR $(z=13)$. Here, we see that both $[\Delta^{s}_{0}]^{3}$
and $[\Delta^{s}_{2}]^{3}$ are positive, and the value of the quadrupole is
comparable to that of the monopole. In some cases, the quadrupole exceeds the
monopole by a factor in the range of $2$ to $3$. We further see that at all
scales $[\Delta^{s}_{0}]^{3}$ and $[\Delta^{s}_{2}]^{3}$ closely follow the
2PT predictions $[\Delta^{s}_{0}]_{\rm 2PT}^{3}$ and $[\Delta^{s}_{2}]_{\rm
2PT}^{3}$ respectively. This indicates that it is reasonable to assume that at
this stage the H i traces the underlying matter ($M$) distribution, whereby
$[\Delta^{s}_{0}]^{3}\approx[\Delta^{s}_{0}]_{M}^{3}$ and
$[\Delta^{s}_{2}]^{3}\approx[\Delta^{s}_{2}]_{M}^{3}$ which arise nearly
entirely from the non-Gaussianity in the underlying matter distribution. Note
that we have also used $[\Delta^{s}_{0}]_{\rm 2PT}^{3}$ and
$[\Delta^{s}_{2}]_{\rm 2PT}^{3}$ as proxies for $[\Delta^{s}_{0}]_{M}^{3}$ and
$[\Delta^{s}_{2}]_{M}^{3}$ respectively at lower redshift. While this is not
strictly correct once there is significant reionization, we still expect this
to correctly capture some of the salient features. In particular, we see that
both $[\Delta^{s}_{0}]_{M}^{3}$ and $[\Delta^{s}_{2}]_{M}^{3}$ increase by
several orders of magnitude as we go from large scale to the smallest scale
(increasing $k$), a behaviour we may also expect to hold at lower $z$. In
addition, we also expect both $[\Delta^{s}_{0}]_{M}^{3}$ and
$[\Delta^{s}_{2}]_{M}^{3}$ to scale as $\bar{x}_{\rm\text{H\,{i}}}^{3}$ as
reionization proceeds.

The intermediate and mid stages of EoR $(z=11-8)$. We find that both
$[\Delta^{s}_{0}]^{3}$ and $[\Delta^{s}_{2}]^{3}$ deviate from the 2PT
predictions during these stages. We attribute these deviations to the
emergence of “extra ionization ($I$) features" in the H i distribution and use
$[\Delta^{s}_{0}]^{3}\approx[\Delta^{s}_{0}]_{M}^{3}+[\Delta^{s}_{0}]_{I}^{3}$
to qualitatively interpret the results for the monopole. We note that this
simple model ignores several ‘matter-ionization’ cross-terms. However, earlier
studies (Majumdar et al., 2018, 2020) have shown that these cross-terms make a
sub-dominant contribution to the monopole, and the two dominant terms
considered here provide a reasonable qualitative understanding. As
reionization proceeds, in all cases, $[\Delta^{s}_{0}]^{3}$ falls below
$[\Delta^{s}_{0}]_{\rm 2PT}^{3}$, and it continues to decline to become
negative at large and intermediate scales. The sign-flip to negative values
occurs at the large scale first, followed by the intermediate scale, and for
both these length-scales and for both squeezed and stretched triangles,
$[\Delta^{s}_{0}]^{3}$ has particularly large negative values in the redshift
range $z=8-9$. For the small and smallest scales, $[\Delta^{s}_{0}]^{3}$
remains positive throughout.

The 21-cm signal during the intermediate and mid stages of EoR may be thought
to originate from a distribution of distinct ionized bubbles embedded in a
neutral medium (Figure 1). Bharadwaj & Pandey (2005) have calculated the
bispectrum for the distribution of non-overlapping ionized spheres in a
neutral medium, and this is predicted to be negative. Further,
$\mid[\Delta^{s}_{0}]_{I}^{3}\mid$ is also predicted to increase as $k^{6}$
for small $k$, and remain constant or decline rapidly for $k\geq\pi/R$, where
$R$ is the typical bubble radius (which varies with $z$). During these stages
where $\bar{x}_{\rm\text{H\,{i}}}\geq 0.5$, we also expect
$\mid[\Delta^{s}_{0}]_{I}^{3}\mid$ to increase as reionization proceeds. These
considerations lead to a picture where the negative values of
$[\Delta^{s}_{0}]^{3}$ arise due to the emergence of the ionized bubbles
whereby $\mid[\Delta^{s}_{0}]_{I}^{3}\mid>[\Delta^{s}_{0}]_{M}^{3}$ at large
and mediate scales. As noted earlier, the value of $[\Delta^{s}_{0}]_{M}^{3}$
increases with $k$ and it exceeds that of $\mid[\Delta^{s}_{0}]_{I}^{3}\mid$
at the small and smallest scales where $[\Delta^{s}_{0}]^{3}$ is positive
throughout. Several earlier works (Majumdar et al., 2018, 2020; Hutter et al.,
2020) have proposed that the sign-flip from positive to negative values seen
at large and intermediate scales is a distinct signature of the onset of
reionization.

We now consider $[\Delta^{s}_{2}]^{3}$ during the intermediate and mid stages
of EoR. Similar to $[\Delta^{s}_{0}]^{3}$, the values of
$[\Delta^{s}_{2}]^{3}$ also deviate from the 2PT predictions, and these
deviations mostly increase as reionization proceeds. However, unlike
$[\Delta^{s}_{0}]^{3}$, we cannot use
$[\Delta^{s}_{2}]^{3}\approx[\Delta^{s}_{2}]_{M}^{3}+[\Delta^{s}_{2}]_{I}^{3}$
to qualitatively interpret the evolution. It is important to note that
$[\Delta^{s}_{2}]^{3}$ arises due to RSD, and $[\Delta^{s}_{2}]^{3}$ is
predicted to be zero in the absence of RSD. The peculiar velocities
responsible for RSD are sourced by fluctuations in the underlying matter
distribution, and a reference to the matter fluctuations is unavoidable for
modeling $[\Delta^{s}_{2}]^{3}$. Here we use
$[\Delta^{s}_{2}]^{3}\approx[\Delta^{s}_{2}]_{M}^{3}+[\Delta^{s}_{2}]_{IM}^{3}$
to qualitatively interpret the results for the quadrupole. The term
$[\Delta^{s}_{2}]_{IM}^{3}$ encapsulates many possible matter-ionization
cross-correlations which arise due to the effect of peculiar velocities on the
ionization field. Notable among these is a term which arises from the
rearrangement of the ionization field due to the peculiar velocities, a
similar term for the matter field also arises in the 2PT predictions
(Scoccimarro et al., 1999).

We first consider the large scale for squeezed triangles. Here we see that the
value of $[\Delta^{s}_{2}]^{3}$ matches the 2PT prediction
$[\Delta^{s}_{2}]_{\rm 2PT}^{3}$ at $z=11$. We also notice that
$[\Delta^{s}_{2}]^{3}\approx 2~{}|[\Delta^{s}_{0}]^{3}|$ where
$[\Delta^{s}_{0}]^{3}$ is negative at $z=11$. This indicates that at $z=11$ we
already have ionized bubbles which cause $[\Delta^{s}_{0}]_{I}^{3}$ to become
negative and exceed $[\Delta^{s}_{0}]_{M}^{3}$ in magnitude. However, this
does not appear to have much impact on $[\Delta^{s}_{2}]^{3}$ implying that
$[\Delta^{s}_{2}]_{IM}^{3}\ll[\Delta^{s}_{2}]_{M}^{3}$. At lower $z$ $(10\leq
z\leq 8)$, the values of $[\Delta^{s}_{2}]^{3}$ drop below those of
$[\Delta^{s}_{2}]_{\rm 2PT}^{3}$ and they are also small in comparison to the
monopole $([\Delta^{s}_{2}]^{3}\leq 0.3\mid[\Delta^{s}_{0}]^{3}\mid)$ which is
negative.

We next consider the intermediate scales. Here $[\Delta^{s}_{0}]^{3}$ remains
positive at $z=11$ and $10$ for both squeezed and stretched triangles. For
squeezed triangles $[\Delta^{s}_{2}]^{3}$ is below $[\Delta^{s}_{2}]_{\rm
2PT}^{3}$, however the values are quite large with
$[\Delta^{s}_{2}]^{3}\approx[\Delta^{s}_{0}]^{3}$ and $\approx
2\,[\Delta^{s}_{0}]^{3}$ at $z=11$ and $10$ respectively. For stretched
triangles, $[\Delta^{s}_{2}]^{3}$ matches $[\Delta^{s}_{2}]_{\rm 2PT}^{3}$ and
$\approx[\Delta^{s}_{0}]^{3}$ at $z=11$, whereas $[\Delta^{s}_{2}]^{3}$
slightly exceeds $[\Delta^{s}_{2}]_{\rm 2PT}^{3}$ and $\approx
10\,[\Delta^{s}_{0}]^{3}$ at $z=10$. These results seem to indicate that we
have $[\Delta^{s}_{0}]_{I}^{3}<[\Delta^{s}_{0}]_{M}^{3}$ and
$[\Delta^{s}_{2}]_{IM}^{3}<[\Delta^{s}_{2}]_{M}^{3}$ at $z=11$ and $10$,
although the reason why $[\Delta^{s}_{2}]^{3}$ slightly exceeds
$[\Delta^{s}_{2}]_{\rm 2PT}^{3}$ for stretched triangles at $z=10$ is not
clear at present. We see that $[\Delta^{s}_{0}]^{3}$ drops sharply and becomes
negative at lower redshifts $(z=9,8)$ for both squeezed and stretched
triangles. In contrast, $[\Delta^{s}_{2}]^{3}$ declines gradually and mostly
remain positive in this $z$ range. Here, the values of $[\Delta^{s}_{2}]^{3}$
are relatively small in comparison to the monopole, and they are in the range
$[\Delta^{s}_{2}]^{3}\leq 0.4\mid[\Delta^{s}_{0}]^{3}\mid$.

Considering the large and intermediate scales together, the monopole
$[\Delta^{s}_{0}]^{3}$ exhibits a transition from positive to negative as
reionization proceeds. In contrast , the quadrupole $[\Delta^{s}_{2}]^{3}$
mostly remains positive. Although $[\Delta^{s}_{2}]^{3}$ does become negative
at $z=8$ for two cases, the error bars are large, and these values are
consistent with zero.

For the small and smallest scales, $[\Delta^{s}_{2}]^{3}$ roughly matches the
2PT predictions at $z=11$ and deviates from these as EoR progresses. At these
scales the values of $[\Delta^{s}_{2}]^{3}$ remain positive throughout. For
squeezed triangles, the ratio
$\mid[\Delta^{s}_{2}]^{3}/[\Delta^{s}_{0}]^{3}\mid$ has values in the range
$0.5-1.5$, with the exception of $z=9$ where it shoots up to $\approx 2.5$ for
small length scales. For stretched triangles, the ratio
$\mid[\Delta^{s}_{2}]^{3}/[\Delta^{s}_{0}]^{3}\mid$ has values in the range
$1.2-2.5$, with the exception of $z=8$ where it drops to $\approx 0.8$ for
small length scales. Considering all the length scales together, we note that
the values of $\mid[\Delta^{s}_{2}]^{3}\mid$ are larger for stretched
triangles in comparison to the squeezed triangles.

We reiterate that the quadrupole moment arises entirely due to RSD, and this
is predicted to be zero in the absence of RSD. The 2PT predictions illustrate
that RSD has a strong effect on the underlying matter distribution, as the
induced $[\Delta^{s}_{2}]^{3}_{M}$ is positive and is comparable in magnitude
to $[\Delta^{s}_{0}]^{3}_{M}$. In contrast, RSD is expected to have a less
pronounced effect on the ionized bubbles. It is possible to estimate this by
considering the divergence of the peculiar velocities (in suitable units, eq.
(10) of Bharadwaj & Ali 2004), which is comparable to the density contrast in
the underlying matter distribution. This will be small in comparison to the H
i fluctuations produced by the ionized bubbles, and we expect
$\mid[\Delta^{s}_{2}]_{IM}^{3}\mid\ll\mid[\Delta^{s}_{0}]_{I}^{3}\mid$. This
provides a qualitative interpretation for the fact that $[\Delta^{s}_{2}]^{3}$
matches the 2PT predictions at $z=11$ whereas $[\Delta^{s}_{0}]^{3}$ shows
considerable deviations from $[\Delta^{s}_{0}]_{\rm 2PT}^{3}$ at the same
redshift. This further explains why $[\Delta^{s}_{2}]^{3}$ remains positive at
large and intermediate scales, in contrast to $[\Delta^{s}_{0}]^{3}$, which
becomes negative. This also leads us to expect that
$[\Delta^{s}_{2}]^{3}\approx[\Delta^{s}_{2}]^{3}_{M}$ at the small and
smallest scales.

The late stage of EoR $(z=7)$. Considering the large and intermediate scales,
we see that both $[\Delta^{s}_{0}]^{3}$ and $[\Delta^{s}_{2}]^{3}$ flip signs
to respectively become positive and negative at this stage. Here the magnitude
of $[\Delta^{s}_{2}]^{3}$ is around $0.3\,[\Delta^{s}_{0}]^{3}$. At the small
and smallest scales, both $[\Delta^{s}_{0}]^{3}$ and $[\Delta^{s}_{2}]^{3}$
remain positive at this stage. Here the magnitude of $[\Delta^{s}_{2}]^{3}$ is
comparable to that of $[\Delta^{s}_{0}]^{3}$ at the smallest scale, whereas it
is approximately an order of magnitude smaller than $[\Delta^{s}_{0}]^{3}$ at
the small scale.

During the late stage of EoR, the 21-cm signal originates from individual H i
islands embedded in an ionized background (Figure 1). This is practically a
negative image of the H i distribution during the intermediate stages of EoR,
where we have individual ionized bubbles in a neutral medium, and we expect
$[\Delta^{s}_{0}]^{3}_{I}$ to be similar but with an opposite sign.
Considering the inside-out reionization scenario implemented here, we expect
the over-dense region to be ionized first with the H i islands surviving
mainly in the under-dense regions of the underlying matter distribution.
Considering a hypothetical situation where the H i islands reside in the over-
dense regions, RSD would then produce a positive $[\Delta^{s}_{2}]_{IM}^{3}$
somewhat similar to $[\Delta^{s}_{2}]_{M}^{3}$. However, the peculiar
velocities responsible for the RSD are reversed in the under-dense regions
where the H i islands actually reside, and this leads to a negative quadrupole
$[\Delta^{s}_{2}]_{IM}^{3}$ during the late stage of EoR. For the large and
intermediate scales, both $[\Delta^{s}_{0}]_{M}^{3}$ and
$[\Delta^{s}_{2}]_{M}^{3}$ are expected to be quite small, and the
considerations outlined above for $[\Delta^{s}_{0}]_{I}^{3}$ and
$[\Delta^{s}_{2}]_{IM}^{3}$ provide a qualitative interpretation for the
positive and negative values of $[\Delta^{s}_{0}]^{3}$ and
$[\Delta^{s}_{2}]^{3}$ respectively. The magnitude of
$[\Delta^{s}_{0}]_{M}^{3}$ and $[\Delta^{s}_{2}]_{M}^{3}$ both increase with
$k$, and we expect these to dominate at the small and smallest scales where
$[\Delta^{s}_{0}]^{3}$ and $[\Delta^{s}_{2}]^{3}$ remain positive throughout
EoR.

##  4 Discussion and Conclusions

We have used $45$ independent realizations of the simulated EoR 21-cm signal
to obtain statistically reliable estimates for both the monopole
$(\bar{B}^{0}_{0})$ and quadrupole $(\bar{B}^{0}_{2})$ moments of the BS, and
also the $1\sigma$ error-bars for these. Results are presented only for
squeezed and stretched triangles where the BS is expected to have the largest
magnitude. To the best of our knowledge, this is the first study of the
quadrupole moment of the EoR 21-cm BS. The quadrupole moment arises entirely
due to RSD, and this is predicted to be zero in the absence of RSD.

We first consider the large and intermediate scales ($k=0.29,0.56\,{\rm
Mpc}^{-1}$ respectively) which are particularly sensitive to reionization. At
the early stage, both $\bar{B}^{0}_{0}$ and $\bar{B}^{0}_{2}$ are positive and
they match the 2PT predictions indicating that they trace the non-Gaussianity
of the underlying matter distribution. The subsequent evolution of
$\bar{B}^{0}_{0}$ and $\bar{B}^{0}_{2}$ is punctuated by two changes in sign,
each of which marks a transition in the H i distribution. The first sign-flip
where $\bar{B}^{0}_{0}$ becomes negative occurs during the intermediate stages
of EoR $(\bar{x}_{\rm\text{H\,{i}}}>0.5)$. This sign-flip occurs first at
large scales and subsequently at intermediate scales, and it marks the
emergence of distinct ionized bubbles embedded in the neutral background. It
has been proposed that the negative value of $\bar{B}^{0}_{0}$ is a distinct
signature of the onset of reionization (Majumdar et al., 2018, 2020; Hutter et
al., 2020). The ionized bubbles grow and merge as reionization proceeds. This
is reflected in the negative values of $\bar{B}^{0}_{0}$, which also grows as
reionization proceeds, and we find the largest negative values near the mid-
stage of EoR $(\bar{x}_{\rm\text{H\,{i}}}\approx 0.5)$. In comparison to
$\bar{B}^{0}_{0}$, the quadrupole $\bar{B}^{0}_{2}$ is much less affected by
the onset of reionization. In fact, in some cases, $\bar{B}^{0}_{2}$ continues
to match the 2PT predictions even when $\bar{B}^{0}_{0}$ has become negative.
Further, the values of $\bar{B}^{0}_{2}$ mostly remain positive in contrast to
$\bar{B}^{0}_{0}$, which becomes negative after the onset of reionization.

As reionization proceeds, the ionized regions merge, and we are left with a
distribution of distinct H i islands embedded in an ionized background. The
second sign-flip, when both $\bar{B}^{0}_{0}$ and $\bar{B}^{0}_{2}$ change
signs, marks this change in the topology of the H i distribution at the late
stage of EoR $(\bar{x}_{\rm\text{H\,{i}}}<0.5)$. The H i islands are somewhat
like a negative image of the ionized bubbles. This causes the sign of
$\bar{B}^{0}_{0}$ to change from negative to positive at the late stage. Here
we note a very recent work (Raste et al., 2023)222This appeared on arxiv while
the present manuscript was under preparation. which has reported a similar
sign-change arising from H i islands at the late stage of EoR. The work
presented there consider two different reionization models and a variety of
triangle shapes, it is however restricted to the monopole $\bar{B}^{0}_{0}$.
Considering the quadrupole $\bar{B}^{0}_{2}$, we see that this also exhibits a
sign-flip and changes from positive to negative at the late stage of EoR. The
H i islands are expected to survive only in the under-dense regions of the
underlying matter distribution. This causes the peculiar velocities to be
anti-correlated with the H i distribution, resulting in a negative
$\bar{B}^{0}_{2}$. Note that $\bar{B}^{0}_{2}$ is negative only if the H i
islands reside in under-dense regions, and it would be positive if the H i
islands were to reside in over-dense regions instead. Hence, unlike
$\bar{B}^{0}_{0}$ which only captures the transition in topology from the
emergence of H ii bubbles to H i islands, $\bar{B}^{0}_{2}$ can potentially
distinguish between inside-out and outside-in models of the EoR. We note that
the position-dependent power spectrum studied by Giri et al. (2019) is also a
useful tool to distinguish between such EoR models.

The present work shows that both $\bar{B}^{0}_{0}$ and $\bar{B}^{0}_{2}$ are
sensitive to changes in the topology of the H i distribution. Considering
simulations which exhibit a reionization history very similar to the present
work, Bag et al. (2019) have identified a percolation transition in the H ii
distribution at $z=9$ where $\bar{x}_{\rm\text{H\,{i}}}=0.728$ and barely
$12\%$ of the volume is ionized. It is interesting to note that the BS seems
to indicate a change in topology (or at least a transition) at a lower $z$,
somewhere between $z=7$ and $z=8$ where $\bar{x}_{\rm\text{H\,{i}}}<0.5$. The
quadrupole moment of the large scale 21-cm PS, however, becomes negative for
$\bar{x}_{\rm\text{H\,{i}}}\leq 0.7$ which is close to the percolation
threshold. However, it is not clear if this behaviour of the 21-cm PS is in
anyway related to the topology of the H i distribution.

The quadrupole $\bar{B}^{0}_{2}$ quantifies aspects of the non-Gaussian EoR
21-cm signal that are not captured by the monopole $\bar{B}^{0}_{0}$. A joint
detection of $\bar{B}^{0}_{0}$ and $\bar{B}^{0}_{2}$ would provide a strong
confirmation of the non-Gaussian behavior of the EoR 21-cm signal.
Furthermore, these two independent statistics used together would
substantially enhance our understanding of reionization and impose tighter
constraints on EoR model parameters as compared to the situation where we use
$\bar{B}^{0}_{0}$ alone. However, the exploration of EoR model parameter
estimation falls outside the scope of this work. It would also be interesting
to investigate how $\bar{B}^{0}_{0}$ and $\bar{B}^{0}_{2}$ evolve under
different EoR models where the topology and evolution of the 21-cm signal
differ significantly. We plan to explore this in the future.

In conclusion, we note that $\bar{B}^{0}_{2}$ complements $\bar{B}^{0}_{0}$ as
a probe of the H i distribution. The value of $\bar{B}^{0}_{2}$ can be
comparable to that of $\bar{B}^{0}_{0}$, and it may also exceed it by a factor
in the range $2$ to $3$ – particularly at the small and smallest scales where
both $\bar{B}^{0}_{0}$ and $\bar{B}^{0}_{2}$ are predicted to be positive
throughout EoR. At large and intermediate scales, changes in the signs of
$\bar{B}^{0}_{0}$ and $\bar{B}^{0}_{2}$ mark transitions in the H i
distribution. In particular, we propose that a positive $\bar{B}^{0}_{0}$
accompanied by a negative $\bar{B}^{0}_{2}$ are a distinct and definite
signature of H i islands residing in under-dense regions at the late stage of
reionization.

## Acknowledgements

SSG and SP acknowledge the support of the Prime Minister’s Research Fellowship
(PMRF). AKS acknowledges support by the Israel Science Foundation (grant no.
255/18). SM acknowledges support from the Science and Engineering Research
Board (SERB), Department of Science and Technology, GoI through the Core
Research Grant No. CRG/2021/004025. We acknowledge the National Supercomputing
Mission (NSM) for providing computing resources of ‘PARAM Shakti’ at IIT
Kharagpur, which is implemented by C-DAC and supported by the Ministry of
Electronics and Information Technology (MeitY) and Department of Science and
Technology (DST), Government of India.

## Data availability

The simulated data and the package involved in this work will be shared on
reasonable request to the corresponding author.

## References

  * Abdurashidova et al. (2022) Abdurashidova Z., et al., 2022, ApJ, 925, 221
  * Ali et al. (2008) Ali S. S., Bharadwaj S., Chengalur J. N., 2008, MNRAS, 385, 2166
  * Bag et al. (2019) Bag S., Mondal R., Sarkar P., Bharadwaj S., Choudhury T. R., Sahni V., 2019, MNRAS, 485, 2235
  * Becker et al. (2001) Becker R. H., et al., 2001, AJ, 122, 2850
  * Bharadwaj & Ali (2004) Bharadwaj S., Ali S. S., 2004, MNRAS, 352, 142
  * Bharadwaj & Pandey (2005) Bharadwaj S., Pandey S. K., 2005, MNRAS, 358, 968
  * Bharadwaj et al. (2001) Bharadwaj S., Nath B. B., Sethi S. K., 2001, JoAA, 22, 21
  * Bharadwaj et al. (2020) Bharadwaj S., Mazumdar A., Sarkar D., 2020, MNRAS, 493, 594
  * Choudhury et al. (2009) Choudhury T. R., Haehnelt M. G., Regan J., 2009, MNRAS, 394, 960
  * Furlanetto et al. (2006) Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433, 181
  * Giri et al. (2019) Giri S. K., D’Aloisio A., Mellema G., Komatsu E., Ghara R., Majumdar S., 2019, J. Cosmology Astropart. Phys., 2019, 058
  * Hutter et al. (2020) Hutter A., Watkinson C. A., Seiler J., Dayal P., Sinha M., Croton D. J., 2020, MNRAS, 492, 653
  * Kamran et al. (2021a) Kamran M., Majumdar S., Ghara R., Mellema G., Bharadwaj S., Pritchard J. R., Mondal R., Iliev I. T., 2021a, arXiv e-prints, p. arXiv:2108.08201
  * Kamran et al. (2021b) Kamran M., Ghara R., Majumdar S., Mondal R., Mellema G., Bharadwaj S., Pritchard J. R., Iliev I. T., 2021b, MNRAS, 502, 3800
  * Kamran et al. (2022) Kamran M., Ghara R., Majumdar S., Mellema G., Bharadwaj S., Pritchard J. R., Mondal R., Iliev I. T., 2022, JCAP, 2022, 001
  * Koopmans et al. (2015) Koopmans L., et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 1 (arXiv:1505.07568), doi:10.22323/1.215.0001
  * Majumdar et al. (2013) Majumdar S., Bharadwaj S., Choudhury T. R., 2013, MNRAS, 434, 1978
  * Majumdar et al. (2018) Majumdar S., Pritchard J. R., Mondal R., Watkinson C. A., Bharadwaj S., Mellema G., 2018, MNRAS, 476, 4007
  * Majumdar et al. (2020) Majumdar S., Kamran M., Pritchard J. R., Mondal R., Mazumdar A., Bharadwaj S., Mellema G., 2020, MNRAS, 499, 5090
  * Mazumdar et al. (2020) Mazumdar A., Bharadwaj S., Sarkar D., 2020, MNRAS, 498, 3975
  * Mondal et al. (2015) Mondal R., Bharadwaj S., Majumdar S., Bera A., Acharyya A., 2015, MNRAS: Letters, 449, L41
  * Mondal et al. (2017a) Mondal R., Bharadwaj S., Majumdar S., 2017a, MNRAS, 464, 2992
  * Mondal et al. (2017b) Mondal R., Bharadwaj S., Datta K. K., 2017b, MNRAS, 474, 1390
  * Mondal et al. (2021a) Mondal R., Bharadwaj S., Majumdar S., 2021a, ReionYuga: Epoch of Reionization neutral Hydrogen field generator, Astrophysics Source Code Library, record ascl:2107.005 (ascl:2107.005)
  * Mondal et al. (2021b) Mondal R., Mellema G., Shaw A. K., Kamran M., Majumdar S., 2021b, MNRAS, 508, 3848
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 596, A108
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Pramanick et al. (2023) Pramanick S., Mondal R., Bharadwaj S., 2023, arXiv preprint arXiv:2304.14171
  * Raste et al. (2023) Raste J., Kulkarni G., Watkinson C. A., Keating L. C., Haehnelt M. G., 2023, arXiv e-prints, p. arXiv:2308.09744
  * Scoccimarro et al. (1999) Scoccimarro R., Couchman H. M. P., Frieman J. A., 1999, Astrophys. J., 517, 531
  * Shaw et al. (2019) Shaw A. K., Bharadwaj S., Mondal R., 2019, MNRAS, 487, 4951
  * Shaw et al. (2020) Shaw A. K., Bharadwaj S., Mondal R., 2020, MNRAS, 498, 1480
  * Shaw et al. (2021) Shaw A. K., Bharadwaj S., Sarkar D., Mazumdar A., Singh S., Majumdar S., 2021, JCAP, 2021, 024
  * Tiwari et al. (2022) Tiwari H., Shaw A. K., Majumdar S., Kamran M., Choudhury M., 2022, JCAP, 2022, 045
  * Trenti et al. (2010) Trenti M., Stiavelli M., Bouwens R. J., Oesch P., Shull J. M., Illingworth G. D., Bradley L. D., Carollo C. M., 2010, ApJ, 714, L202
  * Watkinson et al. (2021) Watkinson C. A., Trott C. M., Hothi I., 2021, MNRAS, 501, 367
  * Watkinson et al. (2022) Watkinson C. A., Greig B., Mesinger A., 2022, MNRAS, 510, 3838
  * Yoshiura et al. (2015) Yoshiura S., Shimabukuro H., Takahashi K., Momose R., Nakanishi H., Imai H., 2015, MNRAS, 451, 266