# Velocity recostruction with graph neural networks

Hideki Tanimura<EMAIL_ADDRESS>Center for Data-Driven Discovery,
Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583,
Japan  Albert Bonnefous  École Normale Supérieure de Lyon, 15 parvis René-
Descartes 69342 Lyon Cedex 07, France  Jia Liu  Center for Data-Driven
Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba
277-8583, Japan  Sanmay Ganguly  ICEPP, The University of Tokyo. Hongo,
Bunkyo-ku, Tokyo 113-0033, Japan  Indian Institute of Technology Kanpur,
Kalyanpur, Kanpur 208016, Uttar Pradesh, India

###### Abstract

In this work, we seek to improve the velocity reconstruction of clusters by
using Graph Neural Networks—a type of deep neural network designed to analyze
sparse, unstructured data. In comparison to the Convolutional Neural Network
(CNN) which is built for structured data such as regular grids, GNN is
particularly suitable for analyzing galaxy catalogs. In our GNNs, galaxies as
represented as nodes that are connected with edges. The galaxy positions and
properties—stellar mass, star formation rate, and total number of galaxies
within 100 $h^{-1}\mathrm{Mpc}$—are combined to predict the line-of-sight
velocity of the clusters. To train our networks, we use mock SDSS galaxies and
clusters constructed from the Magneticum hydrodynamic simulations. Our GNNs
reach a precision in reconstructed line-of-sight velocity of $\Delta\varv$=163
km/s, outperforming by $\approx$10% the perturbation theory ($\Delta\varv$=181
km/s) or the CNN ($\Delta\varv$=179 km/s). The stellar mass provides
additional information, improving the precision by $\approx$6% beyond the
position-only GNN, while other properties add little information. Our GNNs
remain capable of reconstructing the velocity field when redshift-space
distortion is included, with $\Delta\varv$=210 km/s which is again 10% better
than CNN with RSD. Finally, we find that even with an impressive, nearly 70%
increase in galaxy number density from SDSS to DESI, our GNNs only show an
underwhelming 2% improvement, in line with previous works using other methods.
Our work demonstrates that, while the efficiency in velocity reconstruction
may have plateaued already at SDSS number density, further improvements are
still hopeful with new reconstruction models such as the GNNs studied here.

##  I Introduction

The velocity field of the large-scale structure is a useful probe of
cosmological and astrophysical parameters. Reconstructing the velocity field
was an essential step in previous works to e.g. sharpen the baryon acoustic
oscillations (BAO) peaks [1, 2, 3] and to extract the kinematic Sunyaev-
Zel’dovich (kSZ) [4] via cross-correlations with the Cosmic Microwave
Background (CMB) measurements [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17].
The BAO signal originates from the acoustic density waves in the primordial
plasma. It is used as a standard ruler to measure the expansion history of the
universe [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36]. The kSZ signal originates from the scattering of CMB photons off
coherently moving electrons. Because the kSZ effect is a particularly powerful
probe of the electron distribution in under-dense regions, it was used to
investigate the “missing baryon problem” [37, 38, 39, 12] and to constrain
baryonic feedback [40, 41, 15, 42, 9, 43].

In the past, the velocity field was typically reconstructed by solving the
linearized continuity equation on the galaxy density field [44, 45, 46, 47,
48, 14]. While the method works well on large scales where density fluctuation
is small (“linear”), three issues complicate the procedure on small scales
($\lesssim$ 30 $h^{-1}\mathrm{Mpc}$). First, structure formation becomes
nonlinear on these scales and hence is not fully captured in the linearized
continuity equation [49, 50]. Second, to obtain the underlying matter field,
we need to correct the galaxy bias, which becomes nonlinear on small scales
and often requires numerical simulations to calibrate [51, 52, 53, 54, 55,
56]. Finally, the line-of-sight (LOS) positions of the galaxies are measured
in redshift, and hence we need to correct for redshift-space distortion (RSD)
to obtain the comoving distance [57, 58].

In this work, we use the Graph Neural Network (GNN)—a deep neural network
particularly suitable to sparse, heterogeneous data—to reconstruct the
velocity field. The application of machine learning to velocity reconstruction
is not new. Previously, [59, 60, 15, 61] performed such tasks using the
Convolutional Neural Network (CNN). However, there are several limitations
with CNNs in handling galaxy catalogs that motivated us to experiment with the
GNN architecture [62]. First, CNNs are optimized for grid-like data but less
so for irregularly structured data such as galaxy catalogs. Second, CNNs focus
mainly on local features and hence likely to overlook global patterns. Third,
CNNs are trained on a single data type, i.e. galaxy positions, and hence
unable to incorporate other potentially useful galaxy features such as color,
shape, stellar mass, and star formation rate. In comparison, GNNs can handle
unstructured data of arbitrary size and complex topology, capture global
patterns, and allow us to incorporate multiple data types such as galaxy
positions and properties. Thanks to these advantages, GNNs with galaxies have
been recently used to predict halo masses [63] and the cosmic matter density
$\Omega_{m}$ [64, 65]. In this work, to quantify the performance of GNNs on
realistic observations, we reconstruct the velocity fields using GNNs trained
on mock galaxies of the Sloan Digital Sky Survey (SDSS) and the Dark Energy
Spectroscopic Instrument (DESI), and compare them with the true velocity
fields.

This paper is organized as follows. Section II describes the mock simulations
used in this study, the general GNN architecture, the specific network design
for our galaxy catalogs, and our training procedure. The results are presented
in Section III. We conclude in Section IV.

##  II Methodology

In this section, we describe the SDSS and DESI [66] mock catalogs used in this
work, the basics of GNN, the specific GNN architecture we adopted to model the
galaxy distribution, and our training procedure.

###  II.1 SDSS and DESI mocks from the Magneticum simulations

We construct SDSS and DESI mock galaxy catalogs from the Magneticum
simulations111http://www.magneticum.org/simulations.html, which are among the
largest cosmological hydrodynamic simulations [67, 68]. They are based on the
standard $\Lambda$CDM cosmology with {$\Omega_{\rm m}$, $\Omega_{\rm b}$,
$h$}={0.272, 0.046, 0.704} [69]. Several simulation boxes with different sizes
and resolutions are available. In this work, we use their Box2b with a box
size of 640 $h^{-1}\mathrm{Mpc}$ and a total number of particles
$N_{p}=2\times 2880^{3}$, with equal numbers of dark matter and baryon
particles. The simulations include post-processed data of galaxy and cluster
catalogs.

To obtain realistic SDSS-like galaxies, we removed clusters with
$M_{500}<10^{13}$ $h^{-1}{\rm M}_{\odot}$ and galaxies with $M_{*}<2.0\times
10^{11}$ $h^{-1}{\rm M}_{\odot}$ from the main catalogs, leaving us with
47,132 clusters and 77,101 galaxies. In addition, to investigate the effect of
varying galaxy number density, we create DESI-like mock catalogs by removing
galaxies with $M_{*}<1.4\times 10^{11}$ $h^{-1}{\rm M}_{\odot}$, leaving
133,721 galaxies in our sample. Throughout this paper, we focus on simulations
at $z\approx$ 0.47222“snapshot 26” for Box2b. to compare our results to those
using CNN as done in [15] (hereafter T22, for more details, see Sec. 3 in
T22).

###  II.2 Graph neural network

A galaxy catalog can be naturally represented by a graph, where galaxies as
nodes are connected by edges. Several features can be added as node
attributes, e.g., 3-dimensional (3D) position, stellar mass, and star
formation rate. An example graph centered at a galaxy cluster is shown in Fig.
1.

With a set of $N$ nodes, a graph can be described with an $N\times N$
adjacency matrix $G_{ij}$, where $G_{ij}=1$ if nodes $i,j$=$1,2,...N$ are
connected by an edge and $G_{ij}=0$ otherwise. The neighborhood of a node $i$
is defined as

$\mathcal{N}_{i}=\\{j\>|\>G_{ij}=1\\}.$ (1)

Each node receives information from its neighborhood via the message-passing
scheme. In this scheme, the neighborhood’s features (“messages”) are passed
through a multilayer perceptron (MLP)—a fully connected neural network—and
aggregated to form hidden feature vectors as

$\mathbf{h}_{i}=\bigoplus_{j\in\mathcal{N}_{i}}\psi(\mathbf{x}_{i},\mathbf{x}_{j}),$
(2)

where $\bigoplus$ is the differentiable, permutation invariant aggregation
function (e.g., the maximum, the mean, or the sum), $\psi$ is the message-
passing MLP, and $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ are the feature vectors
of the node itself and its neighbor $j$, respectively.

To obtain the output quantities $\mathbf{f}(\mathcal{G})$—cluster velocities
in our case—the hidden feature vectors $\mathbf{h}_{i}$ and the global feature
vector $\mathbf{u}$ are passed through a second MLP $\phi$,

$\mathbf{f}(\mathcal{G})=\phi\left(\bigoplus_{i\in\mathcal{G}}\mathbf{h}_{i}\,,\mathbf{u}\right),$
(3)

where $\mathbf{u}$ describes global features of the graph such as the total
number of nodes and that of edges.

Finally, the network is trained to minimize the loss between the output (e.g.,
the peculiar velocity of a cluster) and its true value.

###  II.3 Galaxy distribution as a graph

In this work, we focus on predicting the LOS velocities of clusters, instead
of the full 3D velocities. Galaxy clusters are the largest gravitationally
bound structures in the Universe. The LOS velocity is particularly relevant to
the kSZ observations for upcoming high-resolution CMB experiments such as the
Simons Observatory [70] and CMB-S4 [71] to constrain cosmology and
astrophysics. The kSZ signal is sensitive to the LOS electron momentum,

$\frac{\Delta{T}_{\rm kSZ}}{T_{\rm CMB}}=-\sigma_{\rm T}\int n_{\rm
e}\,\left(\frac{\bm{\varv}\cdot\bm{\hat{n}}}{c}\right)\,{\rm
d}{l}\approx-\tau\left(\frac{\bm{\varv}\cdot\bm{\hat{n}}}{c}\right),$ (4)

where $\sigma_{\rm T}$ is the Thomson scattering cross section, $c$ is the
speed of light, $n_{\rm e}$ is the electron number density, and
$\bm{\varv}\cdot\bm{\hat{n}}$ is the peculiar velocity of electrons along the
LOS. The integral $\tau=\sigma_{\rm T}\int n_{\rm e}{\rm d}{l}$ can be
performed given that the typical correlation length of electron velocity
(given by $\bm{\varv}\cdot\bm{\hat{n}}$), 80–100 $h^{-1}\mathrm{Mpc}$ [10], is
much larger than that of the gas density ($\approx$5 $h^{-1}\mathrm{Mpc}$).
Therefore, the reconstructed LOS velocity field can help break its degeneracy
with the electron density in kSZ observations.

To construct the graph for each cluster, we define the nodes to be galaxies
within 100 $h^{-1}\mathrm{Mpc}$ of the cluster. The scale was chosen to match
the typical correlation length of peculiar velocities of 80–100
$h^{-1}\mathrm{Mpc}$ [10]. Galaxy pairs within 15 $h^{-1}\mathrm{Mpc}$ are
connected by edges333While theoretically, all galaxies can be fully connected
in the graph without the 15 $h^{-1}\mathrm{Mpc}$ limit, realistically the
graph quickly becomes too large to store in the GPU memory..

Figure 1: An example graph centered at one cluster, projected onto a 2D plane.
Galaxies within 100 $h^{-1}$Mpc of the cluster are represented as nodes.
Galaxy pairs within 15 $h^{-1}\mathrm{Mpc}$ are connected by edges. Figure 2:
The GNN architecture used in our work. It consists of a message-passing phase
in red and a global pooling phase in blue (see Sec. II.3 for details). In the
message-passing phase, node $i$ receives messages (3D position, stellar mass,
star formation rate in our training) $\mathbf{x}$ from its neighbors
($a,b,c,d$) and passes them through an MLP $\psi$. These messages are
aggregated to form the hidden feature vector $\mathbf{h}_{i}$. Next, a global
feature $\mathbf{u}$ (number of galaxies in a graph) is added to the hidden
feature vector and passed through a second MLP $\phi$. The network outputs the
predicted LOS galaxy cluster velocity $\mathbf{f}(\mathcal{G})$.

Fig. 1 shows an example of a graph around one cluster. We illustrate the GNN
architecture used in our work in Fig. 2. We adopt edge convolutional layers
[72, 63], where the relative vectors $\mathbf{x}_{i}$-$\mathbf{x}_{j}$ are
combined with the feature vector $\mathbf{x}_{i}$, replacing Eq. 2 with

$\mathbf{h}_{i}=\bigoplus_{j\in
N_{i}}\psi(\mathbf{x}_{i},\mathbf{x}_{i}-\mathbf{x}_{j}).$ (5)

This choice allows us to combine the global shape structure ($\mathbf{x}_{i}$)
with the local information ($\mathbf{x}_{i}-\mathbf{x}_{j}$). Table 1
summarizes the features included for each node and in the global-pooling
phase.

Type | Feature | Symbol
---|---|---
Node | 3D position | $p$
Node | Stellar mass | $m_{*}$
Node | Star formation rate | $m_{SFR}$
Global | Number of galaxies in a graph | $N_{g}$
Table 1: Galaxy properties used in our GNN.

###  II.4 GNN Training

For the training, we first split the Magneticum simulation box of 640
($h^{-1}\mathrm{Mpc}$)3 into eight independent regions and use seven for
training and validation and one for test. The input features such as the 3D
positions and galaxy properties (Table 1) are provided as part of the
simulation. Our network outputs the LOS velocity. We adopt an L2 loss,

${\it Loss}=\dfrac{1}{N}\sum\limits_{i}^{N}(\varv_{\rm LOS,i}^{\rm
truth}-\varv_{\rm LOS,i}^{\rm pred})^{2},$ (6)

where $\varv_{\rm LOS,i}^{\rm truth}$ is the true LOS velocity of i-th galaxy
cluster and $\varv_{\rm LOS,i}^{\rm pred}$ is its predicted value by our GNN.

To optimize the hyperparameters of our network constituting of one aggregation
function $\bigoplus$ and two MLPs, $\psi$ and $\phi$, we use the automatic
hyperparameter optimization software Optuna [73]. For the aggregation function
$\bigoplus$, we use the mean as it outperforms the sum or the maximum. Our
massage-passing MLP ($\phi$) and global-pooling MLP ($\psi$) both consist of
three hidden layers with 200, 200, and 100 channels, with ReLU activation
functions. The output of the network is the LOS velocity of one galaxy
cluster. We adopt the optimization algorithm Adam with a learning rate of
$10^{-3}$ and weight decay of $10^{-6}$.

##  III Results

In this section, we compare the LOS velocities predicted by our GNN to
theoretical predictions based on linear perturbation theory and those obtained
with CNN in T22. We also investigate the impact of RSD and galaxy number
density.

Figure 3: Predicted LOS velocities of mock SDSS clusters by our GNN compared
to the true values. without RSD (left panel) and with RSD (right panel). The
dashed lines indicate when the two values are equal.

###  III.1 Reconstructed LOS velocity

We show in the left panel of Fig. 3 the predicted LOS velocities on the test
set ($\approx$12.5% of the total data), compared with the true velocities. We
note that the predicted LOS velocities by our GNN were $\approx$11% lower than
the true values, and this bias is corrected using the procedure in the
Appendix. A. The predicted and true LOS velocities are positively correlated.
We quantify the precision of our GNN predictions by

$\Delta\varv=\sqrt{\frac{1}{N}\sum\limits_{i}^{N}(\varv_{\rm LOS,i}^{\rm
truth}-\varv_{\rm LOS,i}^{\rm pred})^{2}.}$ (7)

We obtain an uncertainty of $\Delta\varv$=163 km/s. Our GNN outperforms the
CNN applied to the SDSS galaxy mocks ($\Delta\varv$=179 km/s, see T22). This
improvement is owing to the arbitrarily high spatial resolution that GNN can
achieve, compared to the fixed grid size in CNN, and the inclusion of galaxy
properties such as the stellar mass (see below).

###  III.2 Comparison with theory

To evaluate the performance of our GNN, we compare our results to those from
theoretical predictions. The peculiar velocity of a galaxy cluster can be
derived in the linear regime by

$\bm{\varv}(\bm{k})=-f(\Omega)aH(a)\,\frac{i\bm{k}}{k^{2}}\,\delta(\bm{k})\,,$
(8)

where $\bm{\varv}(\bm{k})$ and $\delta(\bm{k})$ are the velocity and matter
density fields in Fourier space, respectively, and $a$ is the scale factor,
$H(a)$ is the Hubble parameter, $f(\Omega)$ is the linear velocity growth rate
given by $f(\Omega)\approx\Omega_{\rm m}^{0.545}$ [74, 75, 76, 77, 78].

We first compute the galaxy density field in 2003 $h^{3}\mathrm{Mpc}^{-3}$
cubes with a grid size of 5 $h^{-1}\mathrm{Mpc}$ centered at each cluster. We
smooth the density field by 15 $h^{-1}\mathrm{Mpc}$ [79] to remove shot noise.
We then compute the 3D peculiar velocities of clusters using Eq. 8. The galaxy
bias is found to be $\approx$ 1.3, which is used to translate from the galaxy
density field to the matter density field. We take the z direction of the 3D
velocity in the simulation to be the LOS velocity. The estimated LOS
velocities are positively correlated with the true values, with an uncertainty
of $\Delta\varv$=181 km/s.

Our GNN again outperforms the linear theory, as it implicitly models smaller
scales and individual galaxy biases without the need to smooth the grid or
take an average bias value.

###  III.3 Contribution of galaxy properties

To estimate the contribution from galaxy properties in our velocity
reconstruction, we remove galaxy features one by one, such as the total number
of galaxies in a graph (within 100 $h^{-1}\mathrm{Mpc}$), the star formation
rate, and the stellar mass (Table 1). The results are presented in Table 2.
Adding the stellar mass reduces the velocity uncertainty by 6%, compared to
the position-only GNN. Intuitively, galaxy mass contains information on the
gravitational potential, e.g. through the stellar-to-halo mass relation.
However, negligible improvements are seen from other features. In our final
reported results, we use all the features listed in Table 1.

Features |  $\Delta\varv$ (km/s) |  $\Delta\varv$ (km/s)
---|---|---
| without RSD | with RSD
$p$ \+ $m_{*}$ \+ $m_{SFR}$ \+ $N_{g}$ | 163 | 210
$p$ \+ $m_{*}$ \+ $m_{SFR}$ | 164 | 211
$p$ \+ $m_{*}$ | 163 | 210
$p$ | 174 | 217
Table 2: Uncertainties in GNN predicted cluster LOS velocities.

###  III.4 Effect of redshift-space distortion

So far, we have input into GNN the galaxy positions in real space. In
observations, however, galaxy distances are measured in redshift space.
Galaxies’ peculiar velocity causes Doppler shifts in addition to the
cosmological expansion, resulting in RSD and therefore introducing
uncertainties in distance measurement.

To investigate the effect of RSD in our model, we retrain our network by
replacing the comoving positions with redshifts. The redshifts are computed
using the LOS velocities of galaxies and clusters. The results are shown in
the right panel of Fig. 3. The predicted LOS velocities using redshifts show a
positive correlation with the true LOS velocities but with a larger
uncertainty of $\Delta\varv$=210 km/s compared to that without RSD (163 km/s).
Even with RSD, our GNN outperforms the CNN results with RSD in T22
($\Delta\varv$=232 km/s) by $\approx$10%.

Figure 4: Relative uncertainties of the predicted LOS velocities of clusters
by our GNN as a function of galaxy density, with RSD included. The blue and
red dashed lines represent the galaxy densities of SDSS and DESI,
respectively.

###  III.5 Effect of galaxy number density

Finally, we investigated the performance of our GNN when the galaxy number
density is varied. The number density of SDSS galaxies is about $3\times
10^{-4}$ $h^{3}\mathrm{Mpc}^{-3}$, or a mean separation distance of
$\approx$15 $h^{-1}\mathrm{Mpc}$ [80]. DESI will see an increase to
$\approx$$5\times 10^{-4}$ $h^{3}\mathrm{Mpc}^{-3}$, or a mean separation
distance of $\approx$12.5 $h^{-1}\mathrm{Mpc}$. We expect an improved GNN
performance as the galaxy number density increases. We test this assumption
with our DESI-like mock (Sec. II.1).

We show in Fig. 4 the resulting uncertainties in LOS velocity predictions. In
addition to the DESI-like galaxies, we vary the galaxy number density from
$1\times 10^{-4}$ to $6\times 10^{-4}$ $h^{3}\mathrm{Mpc}^{-3}$. Improvements
from increasing the galaxy number density are very mild. For example, for a
DESI-like survey, only a $\approx$ 2 % improvement (from $\Delta\varv$=210 to
205 km/s) is seen, despite an almost 70% increase in the number density. This
shows that using the DESI galaxies instead of SDSS galaxies is unlikely to
improve significantly the LOS velocity estimation of clusters.

We compared our results with similar studies done in [8, 9]. They used the
Baryon Oscillation Spectroscopic Survey (BOSS) constant-mass (CMASS) galaxies
DR10 [81] and reconstructed velocities of the galaxies by solving the
linearized continuity equation in redshift space [82, 79]. The performance of
the velocity reconstruction was quantified as a correlation coefficient
between the true and reconstructed galaxy velocities,

$r_{\varv}=\frac{\langle\varv_{\rm true}\varv_{\rm
rec}\rangle}{\sigma_{\rm\varv}^{\rm true}\sigma_{\rm\varv}^{\rm rec}},$ (9)

where $\sigma_{\rm\varv}^{true}$ and $\sigma_{\rm\varv}^{\rm rec}$ are the
standard deviations of the true and reconstructed radial velocities,
respectively. The value of $r_{\varv}$ was found to be $\approx$0.7 by Refs.
[8, 9] when a Wiener filter was applied. When replaced with a Gaussian filter
(as done in [79]), [9] found that $r_{\varv}$ worsened to $\approx$0.5. Our
GNN result shows $r_{\varv}\approx$0.74, slightly outperforming these previous
works. Ref. [83] performed the velocity reconstruction for the mock DESI
galaxies from the AbacusSummit simulation [84] based on the method in [82],
predicting $r_{\varv}\approx$0.73. The authors also found little improvement
when increasing the galaxy number density beyond that of BOSS (see their
Fig.6). Our GNN predicts $r_{\varv}\approx$0.76 for the DESI galaxies, and
hence is in line with previous findings.

##  IV Conclusion

In this work, we investigated possible improvements in reconstructed LOS
velocities of clusters from using Graph Neural Networks. In this study, we
trained our GNNs to predict the LOS velocities of clusters from the 3D
positions and properties of their surrounding galaxies, using realistic SDSS
and DESI mock galaxies.

Our GNNs outperformed both theoretical predictions and CNN by $\approx$10%
(Sec. III), with an uncertainty of $\Delta\varv$=165 km/s. When redshift-space
distortion is included, the predictions degrade to $\Delta\varv$=210 km/s ,
though remain 10% better than CNN. We found that increasing the number density
of galaxies from SDSS to DESI does not significantly improve the LOS velocity
reconstruction, unfortunately, in line with previous findings. Therefore, a
significant improvement in velocity reconstruction may be more hopeful with
improved methods such as the GNNs studied here, rather than a higher galaxy
number density.

GNNs appear to be an exceptionally suitable deep learning architecture to
reconstruct the velocity field, owing to several of its advantages:

  * •

Galaxy properties can be included in addition to their positions. A 6%
improvement was seen from adding galaxy masses (Sec. III.3). However, other
features (Table. 1) did not appear to provide extra information.

  * •

They do not require explicit modeling of the galaxy bias, which is required in
traditional methods based on perturbation theory.

  * •

With or without RSD, they outperform CNNs, likely owing to their ability to
capture arbitrarily small scales, in comparison to CNN which is limited by the
grid size.

In summary, our work demonstrated the potential of GNNs in improving velocity
reconstruction with existing and upcoming galaxy observations such as DESI
[66] and PFS [85]. Improved velocity measurements will drive our understanding
of the gas distribution in our universe when combined with high-precision kSZ
measurements, such as those expected from the upcoming CMB experiments Simons
Observatory [70] and CMB-S4 [71].

###### Acknowledgements.

The authors thank Klaus Dolag and Antonio Ragagnin for providing the
Magneticum simulations. This work was supported by JSPS KAKENHI Grants
23K13095 and 23H00107 (to JL). This work was supported by ILANCE, CNRS -
University of Tokyo, International Research Laboratory, Kashiwa, Chiba
277-8582, Japan. This research used computing resources at Kavli IPMU. SG was
partially supported by Beyond AI and ICEPP at University of Tokyo.

## References

  * [1] J. A. Blazek, J. E. McEwen, and C. M. Hirata,  Phys. Rev. Lett. 116, 121303 (2016), 1510.03554.
  * [2] F. Schmidt,  Phys. Rev. D94, 063508 (2016), 1602.09059.
  * [3] Z. Slepian et al.,  MNRAS 474, 2109 (2018), 1607.06098.
  * [4] R. A. Sunyaev and I. B. Zeldovich,  ARA&A 18, 537 (1980).
  * [5] A. Cooray,  Phys. Rev. D64, 063514 (2001), astro-ph/0105063.
  * [6] S. DeDeo, D. N. Spergel, and H. Trac,  arXiv e-prints , astro (2005), astro-ph/0511060.
  * [7] J. Shao, P. Zhang, W. Lin, Y. Jing, and J. Pan,  MNRAS 413, 628 (2011), 1004.1301.
  * [8] E. Schaan et al.,  Phys. Rev. D93, 082002 (2016), 1510.06442.
  * [9] E. Schaan et al.,  Phys. Rev. D103, 063513 (2021), 2009.05557.
  * [10] Planck Collaboration,  A&A 586, A140 (2016), 1504.03339.
  * [11] M. Münchmeyer, M. S. Madhavacheril, S. Ferraro, M. C. Johnson, and K. M. Smith,  Phys. Rev. D100, 083508 (2019), 1810.13424.
  * [12] S. H. Lim, H. J. Mo, H. Wang, and X. Yang,  ApJ 889, 48 (2020), 1912.10152.
  * [13] N.-M. Nguyen, J. Jasche, G. Lavaux, and F. Schmidt,  JCAP 2020, 011 (2020), 2007.13721.
  * [14] H. Tanimura, S. Zaroubi, and N. Aghanim,  A&A 645, A112 (2021), 2007.02952.
  * [15] H. Tanimura, N. Aghanim, V. Bonjean, and S. Zaroubi,  A&A 662, A48 (2022), 2201.01643.
  * [16] B. Bolliet, J. Colin Hill, S. Ferraro, A. Kusiak, and A. Krolewski,  JCAP 2023, 039 (2023), 2208.07847.
  * [17] J. Cayuso, R. Bloch, S. C. Hotinli, M. C. Johnson, and F. McCarthy,  JCAP 2023, 051 (2023), 2111.11526.
  * [18] W. J. Percival et al.,  MNRAS 381, 1053 (2007), 0705.3323.
  * [19] W. J. Percival et al.,  MNRAS 401, 2148 (2010), 0907.1660.
  * [20] F. Beutler et al.,  MNRAS 416, 3017 (2011), 1106.3366.
  * [21] C. Blake et al.,  MNRAS 418, 1707 (2011), 1108.2635.
  * [22] L. Anderson et al.,  MNRAS 427, 3435 (2012), 1203.6594.
  * [23] L. Anderson et al.,  MNRAS 441, 24 (2014), 1312.4877.
  * [24] K. S. Dawson et al.,  AJ 145, 10 (2013), 1208.0022.
  * [25] A. J. Ross et al.,  MNRAS 449, 835 (2015), 1409.3242.
  * [26] S. Alam et al.,  MNRAS 470, 2617 (2017), 1607.03155.
  * [27] T. M. C. Abbott et al.,  MNRAS 480, 3879 (2018), 1711.00403.
  * [28] T. M. C. Abbott et al.,  MNRAS 483, 4866 (2019), 1712.06209.
  * [29] M. Ata et al.,  MNRAS 473, 4773 (2018), 1705.06373.
  * [30] J. E. Bautista et al.,  ApJ 863, 110 (2018), 1712.08064.
  * [31] J. E. Bautista et al.,  MNRAS 500, 736 (2021), 2007.08993.
  * [32] V. de Sainte Agathe et al.,  A&A 629, A85 (2019), 1904.03400.
  * [33] M. Blomqvist et al.,  A&A 629, A86 (2019), 1904.03430.
  * [34] H. du Mas des Bourboux et al.,  ApJ 901, 153 (2020), 2007.08995.
  * [35] H. Gil-Marín et al.,  MNRAS 498, 2492 (2020), 2007.08994.
  * [36] J. Hou et al.,  MNRAS 500, 1201 (2021), 2007.08998.
  * [37] J. C. Hill, S. Ferraro, N. Battaglia, J. Liu, and D. N. Spergel,  Phys. Rev. Lett. 117, 051301 (2016), 1603.01608.
  * [38] S. Ferraro, J. C. Hill, N. Battaglia, J. Liu, and D. N. Spergel,  Phys. Rev. D94, 123526 (2016), 1605.02722.
  * [39] D. Munshi, I. T. Iliev, K. L. Dixon, and P. Coles,  MNRAS 463, 2425 (2016), 1511.03449.
  * [40] H. Park, M. A. Alvarez, and J. R. Bond,  ApJ 853, 121 (2018), 1710.02792.
  * [41] J. Kuruvilla, N. Aghanim, and I. G. McCarthy,  A&A 644, A170 (2020), 2010.05911.
  * [42] S. Amodeo et al.,  Phys. Rev. D103, 063514 (2021), 2009.05558.
  * [43] B. K. K. Lee, W. R. Coulton, L. Thiele, and S. Ho,  MNRAS 517, 420 (2022), 2205.01710.
  * [44] S. Zaroubi, Y. Hoffman, K. B. Fisher, and O. Lahav,  ApJ 449, 446 (1995), astro-ph/9410080.
  * [45] S. Zaroubi, A. Dekel, Y. Hoffman, and T. Kolatt,  arXiv e-prints , astro (1996), astro-ph/9603068.
  * [46] S. Zaroubi, Y. Hoffman, and A. Dekel,  ApJ 520, 413 (1999), astro-ph/9810279.
  * [47] E. Branchini et al.,  MNRAS 308, 1 (1999), astro-ph/9901366.
  * [48] H. Wang, H. J. Mo, X. Yang, and F. C. van den Bosch,  MNRAS 420, 1809 (2012), 1108.1008.
  * [49] F.-S. Kitaura, R. E. Angulo, Y. Hoffman, and S. Gottlöber,  MNRAS 425, 2422 (2012), 1111.6629.
  * [50] F.-S. Kitaura and R. E. Angulo,  MNRAS 425, 2443 (2012), 1111.6617.
  * [51] A. J. Benson, S. Cole, C. S. Frenk, C. M. Baugh, and C. G. Lacey,  MNRAS 311, 793 (2000), astro-ph/9903343.
  * [52] A. A. Berlind and D. H. Weinberg,  ApJ 575, 587 (2002), astro-ph/0109001.
  * [53] D. H. Weinberg, R. Davé, N. Katz, and L. Hernquist,  ApJ 601, 1 (2004), astro-ph/0212356.
  * [54] S. E. Nuza et al.,  MNRAS 432, 743 (2013), 1202.6057.
  * [55] F.-S. Kitaura et al.,  MNRAS 456, 4156 (2016), 1509.06400.
  * [56] S. A. Rodríguez-Torres et al.,  MNRAS 460, 1173 (2016), 1509.06404.
  * [57] J. C. Jackson,  MNRAS 156, 1P (1972), 0810.3908.
  * [58] N. Kaiser,  MNRAS 227, 1 (1987).
  * [59] S. E. Hong, D. Jeong, H. S. Hwang, and J. Kim,  ApJ 913, 76 (2021), 2008.01738.
  * [60] Z. Wu et al.,  ApJ 913, 2 (2021), 2105.09450.
  * [61] Z. Wu et al.,  MNRAS 522, 4748 (2023), 2301.04586.
  * [62] J. Zhou et al.,  arXiv e-prints , arXiv:1812.08434 (2018), 1812.08434.
  * [63] P. Villanueva-Domingo et al.,  ApJ 935, 30 (2022), 2111.08683.
  * [64] N. S. M. de Santi et al.,  ApJ 952, 69 (2023), 2302.14101.
  * [65] H. Shao et al.,  ApJ 956, 149 (2023), 2302.14591.
  * [66] DESI Collaboration,  arXiv e-prints , arXiv:1611.00036 (2016), 1611.00036.
  * [67] M. Hirschmann et al.,  MNRAS 442, 2304 (2014), 1308.0333.
  * [68] K. Dolag,  The Magneticum Simulations, from Galaxies to Galaxy Clusters,  in IAU General Assembly, p. 2250156, 2015.
  * [69] E. Komatsu et al.,  ApJS 192, 18 (2011), 1001.4538.
  * [70] P. Ade et al.,  JCAP 2019, 056 (2019), 1808.07445.
  * [71] K. Abazajian et al.,  arXiv e-prints , arXiv:1907.04473 (2019), 1907.04473.
  * [72] Y. Wang et al.,  arXiv e-prints , arXiv:1801.07829 (2018), 1801.07829.
  * [73] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama,  arXiv e-prints , arXiv:1907.10902 (2019), 1907.10902.
  * [74] O. Lahav, P. B. Lilje, J. R. Primack, and M. J. Rees,  MNRAS 251, 128 (1991).
  * [75] L. Wang and P. J. Steinhardt,  ApJ 508, 483 (1998), astro-ph/9804015.
  * [76] E. V. Linder,  Phys. Rev. D72, 043529 (2005), astro-ph/0507263.
  * [77] D. Huterer and E. V. Linder,  Phys. Rev. D75, 023519 (2007), astro-ph/0608681.
  * [78] P. G. Ferreira and C. Skordis,  Phys. Rev. D81, 104020 (2010), 1003.4231.
  * [79] M. Vargas-Magaña, S. Ho, S. Fromenteau, and A. J. Cuesta,  Mon. Not. Roy. Astron. Soc. 467, 2331 (2017), 1509.06384.
  * [80] B. Reid et al.,  MNRAS 455, 1553 (2016), 1509.06529.
  * [81] C. P. Ahn et al.,  ApJS 211, 17 (2014), 1307.7735.
  * [82] N. Padmanabhan et al.,  MNRAS 427, 2132 (2012), 1202.0090.
  * [83] B. Ried Guachalla, E. Schaan, B. Hadzhiyska, and S. Ferraro,  arXiv e-prints , arXiv:2312.12435 (2023), 2312.12435.
  * [84] N. A. Maksimova et al.,  MNRAS 508, 4017 (2021), 2110.11398.
  * [85] N. Tamura et al.,  Prime Focus Spectrograph (PFS) for the Subaru telescope: overview, recent progress, and future perspectives,  in Ground-based and Airborne Instrumentation for Astronomy VI, edited by C. J. Evans, L. Simard, and H. Takami, , Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9908, p. 99081M, 2016, 1608.01075.
  * [86] S. Appleby, R. Davé, D. Sorini, C. C. Lovell, and K. Lo,  MNRAS 525, 1167 (2023), 2301.02001.

##  Appendix A Bias correction

Figure 5: Distribution of the true and predicted LOS velocities of clusters
before (upper panel) and after (lower panel) the bias correction.

Figure 6: Relationship of the true and predicted LOS velocities of clusters
before (upper panel) and after (lower panel) the bias correction.

We found a positive correlation between the true and predicted values of the
LOS velocities of clusters, but the predicted values were slightly biased
lower than the true values. We corrected this bias by adding random Gaussian
noise to the predictions, as demonstrated by [86]. The scale of Gaussian
random noise was determined to match the standard deviations of the true and
predicted distributions. Fig. 6 shows the distributions of true and predicted
values before and after the bias corrections in the left and right panels,
respectively. In addition, Fig. 6 shows the scatter plots of the true and
predicted values before and after the bias correction in the left and right
panels, respectively. These figures show that the bias was statistically
corrected to match the true and predicted values.