# Reconstructing dark matter distribution with peculiar velocities: Bayesian
forward modelling with corrections for inhomogeneous Malmquist bias

Supranta S. Boruah1, Guilhem Lavaux2 Michael J. Hudson3,4,5
1 Department of Astronomy and Steward Observatory, University of Arizona, 933
N Cherry Ave, Tucson, AZ 85719, USA
2CNRS & Sorbonne Université, UMR7095, Institut d’Astrophysique de Paris,
F-75014, Paris, France
3Waterloo Centre for Astrophysics, University of Waterloo, 200, University Ave
W, Waterloo, ON N2L 3G1
4Department of Physics and Astronomy, University of Waterloo, Waterloo, ON,
N2L 3G1, Canada
5Perimeter Institute for Theoretical Physics, 31 Caroline St N, Waterloo, ON
N2L 2Y5 Contact e-mail<EMAIL_ADDRESS>e-mail:
<EMAIL_ADDRESS>e-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a forward-modelled velocity field reconstruction algorithm that
performs the reconstruction of the mass density field using only peculiar
velocity data. Our method consistently accounts for the inhomogeneous
Malmquist bias using analytic integration along the line-of-sight. By testing
our method on a simulation, we show that our method gives an unbiased
reconstruction of the velocity field. We show that not accounting for the
inhomogeneous Malmquist bias can lead to significant biases in the forward-
modelled reconstructions. We applied our method to a peculiar velocity data
set consisting of the SFI++ and 2MTF Tully-Fisher catalogues and the A2
supernovae compilation, thus obtaining a novel velocity reconstruction in the
local Universe. Our velocity reconstructions have a cosmological power
spectrum consistent with the theoretical expectation. Furthermore, we obtain a
full description of the uncertainties on reconstruction through samples of the
posterior distribution. We validate our velocity reconstruction of the local
Universe by comparing it to an independent reconstruction using the 2M++
galaxy catalogue, obtaining good agreement between the two reconstructions.
Using Bayesian model comparison, we find that our velocity model performs
better than the adaptive kernel smoothed velocity with the same peculiar
velocity data. However, our velocity model does not perform as well as the
velocity reconstruction from the 2M++ galaxy catalogue, due to the sparse and
noisy nature of the peculiar velocity tracer samples. The method presented
here provides a way to include peculiar velocity data in initial condition
reconstruction frameworks.

###### keywords:

Galaxy: kinematics and dynamics – cosmology: observations – large-scale
structure of Universe

††pubyear: 2021††pagerange: Reconstructing dark matter distribution with
peculiar velocities: Bayesian forward modelling with corrections for
inhomogeneous Malmquist bias–A

##  1 Introduction

The large-scale structure of the Universe sources the peculiar velocity of the
galaxies. The study of the peculiar velocity of galaxies in the local Universe
is important in cosmological applications for two reasons: i) peculiar
velocities are the only probe of the growth of large-scale structure in the
low-redshift ($z\lesssim 0.1$) Universe, ii) peculiar velocities are a
nuisance parameter in the measurement of the Hubble parameter $H_{0}$ and thus
need to be corrected for in such measurements. Correcting for the peculiar
velocity contributions to the redshifts in the measurement of the expansion
history relies on using a reconstruction of the peculiar velocity field in the
local Universe. Given the current tension in the measurement of the Hubble
constant, $H_{0}$ (Verde et al., 2019), we need to account for all sources of
uncertainties. In order to obtain unbiased estimates of the peculiar velocity
corrections for the local $H_{0}$ measurement, we need accurate models of the
local peculiar velocity field.

Therefore, reconstruction of the velocity field of the local Universe is
important for these cosmological applications. Many methods of reconstruction
of the peculiar velocity field have been proposed in the literature. We can
categorize the different methods according to two criteria: i) the data used
for the reconstruction, i.e, whether galaxy catalogues or peculiar velocity
tracers are used in the reconstruction and ii) whether the method involves
forward-modelling or a direct inversion from the data.

Most commonly used reconstructions of the local velocity field uses galaxy
catalogues to reconstruct the galaxy density field (e.g, Carrick et al., 2015;
Lilow & Nusser, 2021). Following the reconstruction of the galaxy density
field, the velocity field is calculated using linear perturbation theory.
However, since galaxies are biased tracers of the matter field, we need to fit
for the parameter $\beta=f/b$ to get the velocity field, where $f$ is the
logarithmic growth parameter and $b$ is the linear galaxy bias. This is done
by comparing the reconstructed velocity field to the peculiar velocity data.
Galaxy catalogues are much denser than peculiar velocity catalogues.
Furthermore, the velocity estimates obtained from peculiar velocity tracers
are very noisy. Therefore, reconstruction of the density field is usually an
easier task than the direct reconstruction of the velocity field. However,
connecting the galaxy density to the total matter density requires us to
assume a model of the galaxy bias, which is an important source of systematic
uncertainty in these reconstruction. On the other hand, using the peculiar
velocity tracer directly for reconstruction does not require that assumption,
thus providing a complementary method for reconstruction of the velocity
field. One well-known example is the POTENT method (Bertschinger & Dekel,
1989; Bertschinger et al., 1990) which uses peculiar velocity data directly to
reconstruct the velocity field assuming that the velocity field is the
gradient of the gravitational potential field. Another such widely used method
is the adaptive kernel-smoothed velocity reconstruction method (Springob et
al., 2014, 2016) that smooths the peculiar velocity data to obtain the
velocity reconstruction.

The other criteria that we can use to categorize the reconstruction methods is
whether the reconstruction is performed using forward-modelling or uses a
direct inversion from the data. Inverting non-linear problems from partial,
noisy, observations is an ill-posed inverse problem, which makes forward-
modelled methods particularly suitable for the task of reconstruction of high-
dimensional fields. Forward-modelled reconstruction methods have become
increasingly popular in cosmology and have been applied in a range of
different applications such as initial conditions reconstruction (Jasche &
Wandelt, 2013; Modi et al., 2018; Jasche & Lavaux, 2019), weak lensing
(Porqueres et al., 2021b, a; Fiedorowicz et al., 2021) and CMB lensing (Millea
et al., 2020a, b). Such methods have also been used for the local velocity
field reconstruction. The simplest of such methods uses a Wiener filtering
technique (Zaroubi et al., 1999). This approach assumes that the
density/velocity field is described as a Gaussian random field and the Wiener
filtered reconstruction is the maximum-a-posteriori (MAP) solution for the
problem. The Wiener filtering approach has been extended to account for
uncertainties using a constrained realization approach (Hoffman & Ribak, 1991;
Hoffman et al., 2015, 2018; Lilow & Nusser, 2021) and to account for biases in
Wiener filtering using the unbiased minimal variance approach (Zaroubi, 2002).
Another similar approach is the Bayesian hierarchical method, virbius (Lavaux,
2016) which is based on the constrained realization approach but accounts for
many different systematic effects in its analysis. This approach has been been
applied to the Cosmicflows-3 (Tully et al., 2016) dataset by Graziani et al.
(2019). However, this method fails to account for the inhomogeneous Malmquist
bias (IHM) which is an important source of systematic error in peculiar
velocity analysis. The IHM bias arises from an incorrect assumption on the
distribution of peculiar velocity tracers due to neglecting the line-of-sight
inhomogeneities. In this paper, we introduce a forward-modelled reconstruction
method that uses peculiar velocity data in its reconstruction while
consistently accounting for the inhomogeneous Malmquist bias in the analysis.
As we show in this paper, not accounting for the IHM bias in the
reconstruction may lead to substantial biases in the reconstruction. One
approach to correct for the IHM bias is by marginalizing over an accurate
model of galaxy distribution (e.g, Hudson, 1994). In this work, we use such an
approach to deal with the inhomogeneous Malmquist bias within our forward-
modelled peculiar velocity reconstruction method. Our reconstruction results
in samples of reconstruction that are compatible with the data from a given
peculiar velocity catalogue, and sampled from a posterior assuming a
$\Lambda$CDM Gaussian prior. Furthermore, the full correlated uncertainties in
the reconstruction can be estimated from the posterior samples. Our method
presented here can be extended in a straightforward manner to include more
complex gravity models such as initial condition reconstruction methods like
borg (Jasche & Wandelt, 2013; Jasche & Lavaux, 2019).

The paper is structured as follows: in Section 2, we describe the simulations
and data used in this paper. In Section 3, we describe our reconstruction
method in detail. Our method is validated with mock simulated catalogues in
Section 4 followed by an application of our method to reconstruct the velocity
field of the local Universe in Section 5. Finally, after a brief discussion in
Section 6, we summarize our results in Section 7. In Appendix A, we
investigate in further detail, the various sources contributing to the
reconstruction error.

##  2 Simulations and data

In this section, we describe the simulations and the peculiar velocity surveys
used in this work.

###  2.1 VELMASS simulation

In order to validate our method, we use an $N$-body simulation from the
VELMASS suite of simulation.111For more details on the simulation which we
used, see Kodi Ramanah et al. (2019) The simulation was performed in a cubic
box of size $(2~{}h^{-1}$ Gpc)3 with $2048^{3}$ particles with mass
$9.387\times 10^{10}~{}h^{-1}M_{\odot}$. The cosmological parameters used are:
$\Omega_{m}=0.315,\Omega_{b}=0.049,H_{0}=68$ km s-1 Mpc-1, $\sigma_{8}=0.81$,
$n_{s}=0.97$ and $Y_{\text{He}}=0.248$. The rockstar (Behroozi et al., 2013)
halo finding software was used to identify the dark matter halos in the
simulation. We only consider the halos with mass $M>3\times
10^{12}~{}h^{-1}\mathrm{M}_{\odot}$.

In order to create a mock peculiar velocity survey, we populate the halos with
mass $M>3\times 10^{12}~{}h^{-1}\mathrm{M}_{\odot}$ with a standard candle,
with different levels of intrinsic scatter. That is, we assign it an absolute
magnitude $\mathcal{M}=\mathcal{M}_{0}+\epsilon$, where $\mathcal{M}_{0}$ is
the standard candle absolute magnitude and $\epsilon$ is the intrinsic scatter
drawn from a normal distribution,
$\epsilon\sim\mathcal{N}(0,\sigma^{2}_{\text{int}})$. Here,
$\sigma_{\text{int}}$ is the intrinsic scatter. We can calculate the apparent
magnitude from this absolute magnitude and the distance to the halos. We make
mock catalogues with different levels of intrinsic scatter. Finally, we make
an apparent magnitude cut to create a magnitude limited catalogue. From this
magnitude limited sample, we randomly select $7\,000$ halos. While running the
forward likelihood method described in section 3.3, the method requires the
estimated galaxy density field as an input. For this purpose, we construct the
halo density field and smooth it using a Gaussian filter of smoothing length
$4~{}h^{-1}$ Mpc.

###  2.2 Peculiar velocity surveys

In this work, we use peculiar velocity data from various surveys in our
algorithm to reconstruct the velocity field of the local Universe. The three
peculiar velocity data sets that we use in this work are the A2 supernovae
compilation, SFI++ and 2MTF Tully-Fisher (TF) catalogues. Compared to our
treatment of these data-sets in Boruah et al. (2020), we do not fit the TF or
the SNe light curve fitting parameters during the reconstruction. Instead we
use the distances as provided in the catalogues. We however fit for a scaling
factor that may arise due to miscalibration in the fitting of the zero-points
of the TF relations as detailed in section 3.3. Furthermore, in order to limit
the edge effects in the reconstruction, we impose a distance cut of
$d<100~{}h^{-1}$ Mpc. Note that imposing such a distance cut does not
introduce a selection bias in the forward peculiar velocity analysis methods
(Strauss & Willick, 1995). The redshift distribution of peculiar velocity
tracers is shown in Figure 1. In the following, we provide a brief description
of the data used in this work. We refer readers to Boruah et al. (2020, 2021)
for more details on the selection and the treatment of outliers, choosing only
to highlight some main features of the data sets.

Figure 1: Redshift distribution of the tracers in the different peculiar
velocity catalogues used in this work. We have imposed a distance cut of
$d<100~{}h^{-1}$ Mpc on all catalogues.

  * •

A2 Supernovae: In Boruah et al. (2020), we compiled the second amendment (A2)
dataset of nearby supernovae from publicly available supernovae from the CfA
supernovae sample (Hicken et al., 2009), Carnegie Supernovae Project - Data
Release 3 (CSP-DR3, Krisciunas et al., 2017), the Lick observatory Supernova
Survey (LOSS, Ganeshalingam et al., 2013) and the Foundation sample (Foley et
al., 2018; Jones et al., 2019) of supernovae. We only consider the supernovae
with estimated distance, $d<100~{}h^{-1}$ Mpc.

  * •

SFI++: We use peculiar velocity data from the SFI++ (Masters et al., 2006;
Springob et al., 2007) Tully-Fisher catalogue, which is an $I$-band TF survey.
We follow the treatment of Boruah et al. (2020) to remove outliers. For
galaxies in groups, we use the peculiar velocity data for the groups to
suppress non-linear contributions to the peculiar velocity data. We only
consider the galaxies/groups with estimated distance, $d<100~{}h^{-1}$ Mpc.

  * •

2MTF: The 2MTF survey (Masters et al., 2008; Hong et al., 2019) is an all sky
Tully-Fisher in the near-infrared $J,H$ and $K$ bands. The survey is limited
to $cz<10000$ km/s.

With the given distance cuts, we select $345$, $1682$, $556$, $1225$ objects
from the A2 catalogue, SFI++ field galaxy catalogue, SFI++ group catalogue and
2MTF catalogue respectively.

##  3 Methodology

We introduce a novel method for forward modelled velocity reconstruction which
can consistently account for inhomogeneous Malmquist (IHM) bias in the
reconstruction in this paper. In this section, we describe the theory and the
methodology behind our method for reconstructing the velocity field with the
peculiar velocity data.

###  3.1 Peculiar velocity theory

Figure 2: A probabilistic graph showing the dependency of the different
observables and the underlying variables. The observables (redshift,
$cz_{\text{obs}}$ and the distance modulus, $\mu$) are shown with a double
circle. In the forward methods of peculiar velocity analysis (see section
3.3), the redshift of the peculiar velocity tracers are predicted as a
function of the distance estimates. This is reflected in the relationship
between the observables. Note that, we use an external estimate of the galaxy
density field, $\hat{\delta}_{g}$ in order to correct for the inhomogeneous
Malmquist bias.

In the $\Lambda$CDM model, under linear perturbation theory, the present day
peculiar velocity field can be expressed in terms of the density field as
(Peebles, 1980),

${\bm{v}}_{{\bm{k}}}=\frac{ifH_{0}{\bm{k}}}{k^{2}}\delta_{{\bm{k}}},$ (1)

where $f\approx\Omega_{m}^{0.55}$ is the logarithmic growth rate in the
$\Lambda$CDM model, $i^{2}=-1$. Since a factor of $k$ appears in the
denominator, the peculiar velocity field is sensitive to the large-scale
density modes. Since these modes are well described with linear perturbation
theory, linear theory predictions of peculiar velocity provide a fairly
accurate result for the true velocity field. The peculiar velocity estimated
using linear theory has been calibrated with simulations (Carrick et al.,
2015; Hollinger & Hudson, 2021) and the uncertainty in the linear theory
velocity estimates due to the non-linearities was found to be
$\sigma_{\text{NL}}\approx 150$ km s-1. In this work, we use linear
perturbation theory to predict the velocity field. In the future, this can be
extended to include non-linear differentiable velocity field models such as
the Lagrangian perturbation theory (LPT, Bouchet et al., 1995) full particle
mesh or COLA (Tassev et al., 2013) simulations that have previously been used
for forward-modelled density reconstruction (e.g, Jasche & Lavaux, 2019).

###  3.2 Inhomogeneous Malmquist bias

Peculiar velocity analyses are impacted by a number of different statistical
biases, including the homogeneous and inhomogeneous Malmquist biases. The
homogeneous Malmquist bias arises due to the fact that the number of galaxies
grows along the line-of-sight according to the volume factor. In this paper,
we mainly discuss the impact of inhomogeneous Malmquist bias in forward
modelled reconstruction with peculiar velocity data. Therefore, we discuss the
basics of inhomogeneous Malmquist (IHM) bias in some detail in this section.

The IHM bias arises due to the lack of accounting for inhomogeneity in the
radial distribution of the peculiar velocity tracers. Neglecting the
overdensities along the line-of-sight leads to a bias in the inferred
distance, $d_{\text{est}}$ compared to the true distance, $d_{\text{true}}$.
The magnitude of the bias is given as (Lynden-Bell et al., 1988; Strauss &
Willick, 1995),

$\Delta_{\text{IHM}}=\langle d_{\text{true}}\rangle-\langle
d_{\text{est}}\rangle\approx d_{\text{est}}\gamma\Delta_{d}^{2},$ (2)

where, $\Delta_{d}$ is the fractional uncertainty in the distance estimates
and $\gamma$ is the logarithmic slope in the galaxy over-density along the
line-of-sight given by,

$\gamma=\frac{\mathrm{d}\ln n(r)}{\mathrm{d}\ln r}.$ (3)

Since the slope changes sign in front of and behind an overdensity, the
distance of the galaxies in front of the overdensity is underestimated and
those at the back are overestimated. This may be interpreted as a spurious
flow if the IHM bias correction is neglected and thus biases the magnitude of
the velocity field to a higher value compared to the true value. Therefore,
the inferred value of $\beta$ is biased high when the IHM bias is not
accounted for in velocity-velocity comparison methods (e.g. Boruah et al.,
2020). As we will show in section 4, in forward-modelled reconstruction from
peculiar velocity data, this effect shows up as an increment in the inferred
power spectrum of the reconstructed density field if the IHM bias is not
corrected for.

###  3.3 Probability model

The method we present in this paper is a forward-modelled Bayesian method to
reconstruct the velocity field using the data from a peculiar velocity survey.
One such reconstruction method is the virbius algorithm (Lavaux, 2016) which
has previously been applied to the Cosmicflows-3 (Tully et al., 2016) dataset
in Graziani et al. (2019). In their method, the true distance of the peculiar
velocity tracers are used as latent parameters and the value of the distance
is inferred in each Markov Chain Monte Carlo (MCMC) step. The homogeneous
Malmquist bias is dealt with by using a skewed prior on these distances.
However, the inhomogeneous Malmquist bias is not explicitly accounted for in
their analysis. In this work, we account for the inhomogeneous Malmquist bias
by using analytic integration over the line-of-sight uncertainty.
Marginalizing over the true distance also reduces the complexity since, we do
not have to infer their value in each MCMC steps. Furthermore, the posterior
distribution of the peculiar velocity tracer distances can potentially be
highly multi-modal due to the line of sight inhomogeneities, which makes the
analytical marginalization preferable to sampling based methods.

In the forward modelled reconstruction, we want to infer the velocity field,
$\\{{\bm{v}}\\}$ given the peculiar velocity data. The peculiar velocity data
set consists of the observed redshift, $\\{z_{\text{obs}}\\}$, and the
estimate of its distance modulus, $\\{\mu\\}$, for a set of peculiar velocity
tracers. The uncertainty associated with the distance modulus is denoted by
$\sigma_{\mu}$. In linear perturbation theory, the velocity field is
determined in terms of the density field, $\\{\delta\\}$ through equation (1).
Therefore, we can equivalently think about the process described here as a
density reconstruction algorithm. Furthermore, we need an estimate of the
galaxy density field to correct for the inhomogeneous Malmquist bias. For this
work, we provide the estimate of the galaxy density field,
$\\{\hat{\delta}_{g}\\}$, as an input to our method.

Figure 3: A visual comparison of the reconstructed radial velocity field and
the associated uncertainties with the true velocity field in the $Z=0$ plane
in our mock simulation. We show the reconstruction for $3$ different mock
surveys with a distance error of $\Delta_{d}=0.07$ (left), $\Delta_{d}=0.15$
(centre) and $\Delta_{d}=0.20$ (right). The top row is the true velocity field
in the simulation in the $Z=0$ plane. The black dashed circle is at a radius
of $150~{}h^{-1}$ Mpc, which is close to the $95\%$ completeness radius for
each survey scenario. The second and the third rows show the mean velocity for
the samples of our reconstructions with and without inhomogeneous Malmquist
bias corrections respectively. As can be seen, if we do not correct for the
IHM bias, the velocity field in the data region is boosted up compared to the
true velocity field. The bottom two rows shows the standard deviation in the
reconstruction samples with and without the IHM bias correction.

In this work, we use a variant of the forward VELMOD method (Willick et al.,
1997), which was introduced in Pike & Hudson (2005) to model the likelihood.
For simplicity, we call this likelihood the Forward Likelihood. The forward
likelihood method corrects for the inhomogeneous Malmquist bias by accounting
for the density dependent line-of-sight distribution of peculiar velocity
tracers. In the forward likelihood method, the observed redshift is predicted
as a function of the distance estimates. That is, we model the likelihood
using the conditional probability,
$\mathcal{P}(\\{cz_{\text{obs}}\\}|\\{\mu\\},\\{{\bm{v}}\\},\\{\hat{\delta}_{g}\\})$.
We show the dependency of the various variables and the observables in a
probabilistic graph in Figure 2. Using Bayes theorem, the posterior for the
density field $\delta$ can be written as

$\mathcal{P}(\\{\delta\\}|\\{cz_{\text{obs}}\\},\\{\mu\\},\\{\hat{\delta}_{g}\\})\propto\mathcal{P}(\\{cz_{\text{obs}}\\}|\\{\mu\\},\\{\delta\\},\\{\hat{\delta}_{g}\\})\mathcal{P}(\\{\delta\\}).$
(4)

In our runs, we use a $128^{3}$ cubic box where the density is inferred. That
leads to a total of $128^{3}$ density parameters. In order to sample from such
high-dimensional parameter space, we used the Hamiltonian Monte-Carlo (HMC,
Neal, 1993, 1996) algorithm. HMC can be used to sample from very high
dimensional parameter space by using the gradient information of the log-
posterior. In this paper, we implement our code by relying on the jax package
(Bradbury et al., 2018) which includes the gradients of the requisite function
by using automatic differentiation. Sampling with HMC is highly sensitive to
the choice of mass matrix. The optimal choice of mass matrix is to use inverse
of the width of the posterior as the mass matrix (Taylor et al., 2008). We use
automatic differentiation to calculate the second derivative of the posterior
for a number of Gaussian mocks and use it to calculate the mass matrix.

The RHS of equation (4) consists of two terms, i) the prior on the density
field, $\mathcal{P}(\\{\delta\\})$ and ii) The likelihood,
$\mathcal{P}(\\{cz_{\text{obs}}\\}|\\{\mu\\},\\{\delta\\},\\{\hat{\delta}_{g}\\})$.
We will first discuss the prior we use in this work. In this work, we assume a
Gaussian prior on the density field. In the Fourier space, we can write it as

$\mathcal{P}(\\{\delta\\})=\prod_{{\bm{k}}}\frac{1}{\sqrt{2\pi\sigma^{2}_{k}}}\exp\bigg{[}-\frac{|\delta_{{\bm{k}}}|^{2}}{2\sigma^{2}_{k}}\bigg{]},$
(5)

where $\sigma^{2}_{k}=V_{s}P(k)$ is the variance in $\delta_{{\bm{k}}}$. Here,
$P(k)$ is the $\Lambda$CDM power spectrum and $V_{s}$ is the volume of the
simulation box. In this work, we fix our cosmological parameters to the
fiducial cosmological parameters used in the VELMASS simulation. We note that
the present day density field is significantly non-Gaussian. Thus it cannot be
well-described using a Gaussian prior. In the future, we can incorporate the
method presented here with the borg framework which uses a non-linear,
differentiable gravity model to forward model the present day velocity field
from the density field at an initial time when the field is well-described as
a Gaussian field. Nevertheless, due to the sensitivity of the velocity field
to the large-scale density modes, the Gaussian prior still proves to be a
useful approximation.

Assuming that the random components of the radial velocity for the various
peculiar velocity tracers are independent, we can write,

$\mathcal{P}(\\{cz_{\text{obs}}\\}|\\{\mu\\},\\{\delta\\},\\{\hat{\delta}_{g}\\})=\prod_{i=1}^{N_{\text{PV}}}\mathcal{P}(cz_{\text{obs}}^{(i)}|\mu^{(i)},\\{\delta\\},\\{\hat{\delta}_{g}\\})$
(6)

Given the large uncertainties on the distance estimates, we marginalize over
the true line-of-sight distance, $r$, for the peculiar velocity tracer. The
likelihood for the individual peculiar velocity tracers can be expressed as

$\displaystyle\mathcal{P}(cz^{(i)}_{\text{obs}}|\mu^{(i)},$
$\displaystyle\\{\delta\\},\\{\hat{\delta}_{g}\\})$
$\displaystyle=\int_{0}^{\infty}\mathrm{d}r\mathcal{P}(cz^{(i)}_{\text{obs}}|r,\\{{\bm{v}}\\})\mathcal{P}(r|\mu^{(i)},\\{\hat{\delta}_{g}\\}).$
(7)

In the above, $\\{{\bm{v}}\\}$ is the linear velocity field calculated from
the density field $\\{\delta\\}$ using equation (1). The first term inside the
integral is modelled as a Gaussian with the uncertainty in the observed
redshift given by $\sigma_{\text{NL}}$, which is the uncertainty in the
velocity induced by the non-linearities in the density field. The likelihood
function in the integrand is given by

$\mathcal{P}(cz^{(i)}_{\text{obs}}|r,\\{{\bm{v}}\\})=\frac{1}{\sqrt{2\pi\sigma^{2}_{\text{NL}}}}\exp\bigg{[}-\frac{(cz^{(i)}_{\text{obs}}-cz^{(i)}_{\text{pred}}(r,{\bm{v}}))^{2}}{2\sigma^{2}_{\text{NL}}}\bigg{]},$
(8)

where

$1+z_{\text{pred}}(r,{\bm{v}})=[1+z_{\text{cos}}(r)]\bigg{(}1+\frac{v_{r}(r)}{c}\bigg{)}.$
(9)

In the above equation, $z_{\text{cos}}$ denotes the ‘cosmological’ recessional
redshift. In order to get the estimate of the velocity field, we smooth the
linear velocity field using a Gaussian smoothing filter with a smoothing scale
of $R_{\text{smooth}}=4h^{-1}$ Mpc. We use this smoothing scale for all
subsequently reported velocity fields. The velocity field smoothed at this
scale is found to be unbiased in $N$-body simulations and has a velocity
uncertainty of $\sigma_{\text{NL}}\approx 150$ km/s (Carrick et al., 2015;
Hollinger & Hudson, 2021). We fix the value of $\sigma_{\text{NL}}$ to this
value in our inference.

The ‘forward’ method for peculiar velocity analysis is susceptible to
homogeneous and inhomogeneous Malmquist biases (Hudson, 1994; Strauss &
Willick, 1995). As mentioned earlier, the homogeneous Malmquist bias arises if
we do not account for the fact that there are more galaxies at large radius.
Moreover, galaxies are not uniformly distributed in the Universe, instead are
clustered at regions of large overdensities. Neglecting these overdensities
along the line-of-sight leads to the inhomogeneous Malmquist bias as discussed
in section 3.2. We use an estimate of the line-of-sight galaxy density,
$\hat{\delta}_{g}(r)$, to account for the inhomogeneous Malmquist bias in the
radial distribution. Using the galaxy density field, the expected number of
peculiar velocity tracers along the line of sight can be estimated as

$n(r|\\{\hat{\delta}_{g}\\})\propto r^{2}[1+\hat{\delta}_{g}(r)].$ (10)

Accounting for this dependence on the line-of-sight density, we can write the
distribution,

$\mathcal{P}(r|\mu^{(i)},\\{\hat{\delta}_{g}\\})\propto
n(r|\\{\hat{\delta}_{g}\\})\exp\bigg{(}-\frac{[\mu^{(i)}-\mu(r)]^{2}}{2\sigma^{(i)2}_{\mu}}\bigg{)},$
(11)

where $\mu(r)=5\log_{10}(r/10~{}\text{pc})$ and $\sigma_{\mu}^{(i)}$ is the
uncertainty on the estimate of the distance modulus. Note that we assume a
Gaussian uncertainty on the distance modulus. One of the main sources of
systematic uncertainty in estimating the distances to the peculiar velocity
tracers is the calibration of the zero-point of the distance-luminosity
relation. In order to take care of this uncertainty, we introduce a scale
factor, $A$, so that all the distance estimates in a given catalog is scaled
by the same factor, $\mu\rightarrow\mu+5\log_{10}A$ (or equivalently
$d\rightarrow Ad$). We then fit and marginalize over this factor during the
reconstruction in a block sampling scheme. When using multiple catalogues, we
fit a separate scaling factor for different catalogues.

In this work, the estimated galaxy density field, $\\{\hat{\delta}_{g}\\}$ is
an additional input to our reconstruction code and is assumed to be
independent of the density field that is inferred. This is not fully self-
consistent and needs to be improved in the future. We note that we do not use
any cross-correlation information between the estimated galaxy density field
and the inferred velocity field. The estimated galaxy density field simply
provides an estimate for the radial distance of tracer. For the test with
simulation, we use the true underlying particle density to account for the IHM
correction. For our run with the real data, we use an iteratively
reconstructed density field from Carrick et al. (2015).

##  4 Validation with simulations

Figure 4: The power spectrum (top), ratio of the power spectrum to the
theoretical expectation (middle) and the cross-correlation of the
reconstructed density field with the true density field (bottom) for the three
mock surveys with distance uncertainties, $\Delta_{d}=0.07$ (left),
$\Delta_{d}=0.15$ (centre) and $\Delta_{d}=0.20$ (right). In the top panels,
we plot the power spectrum inferred from our reconstruction and compare it
with the true power spectrum of the simulation (shown with black line) and the
theoretical expectation (shown as a black dotted line). The red lines show the
power spectrum of the reconstruction without IHM bias correction and the blue
lines are the corresponding curves with IHM bias correction. The shaded region
show the $68\%$ confidence interval as calculated from the reconstruction
samples. In the middle panels, we plot the ratio of the mean of the power
spectrum of our reconstructions to the theoretical expectation. The error bars
are calculated assuming a Gaussian covariance for the power spectrum. As we
can see the power spectrum of the reconstruction without IHM bias correction
is biased high. The magnitude of the bias is higher for the reconstruction
with larger distance uncertainty. Finally, in the bottom panel, we plot the
cross-correlation of the reconstructed density field with the true density
field of the simulation. The cross-correlation is the largest at large scales
(small $k$) and it decreases to $0$ at small scales (high $k$). However, we
see that if the IHM bias is not corrected, the large scale cross-correlation
goes down, suggesting that these density modes are not reconstructed well.
Figure 5: (Top) Slope of regression of the reconstructed $V_{r}$ on the true
$V_{r}$ values from the simulation, binned in radial shells and (bottom) the
RMS error between the reconstructed velocity and the true velocity. We compare
the velocity reconstruction on grid points of the two reconstruction within a
given radial shell. The three horizontal panels show the results for mock
surveys with three distance uncertainties: $\Delta_{d}=0.07$ (left),
$\Delta_{d}=0.15$ (centre) and $\Delta_{d}=0.20$ (right). In the top panel, we
plot the slopes of regression in radial shells. The red lines show the results
without the IHM bias correction and the blue lines are the results after
correcting for the IHM bias. We see that the the slope of regression, $m\sim
1$ deep inside the data region and then drops to $m\sim 0$ at large radius.
Without the IHM bias correction, the slope can be however be higher than $1$
since IHM bias introduces spurious velocities. In the bottom panel, we plot
the reconstruction error in radial bins. The blue and the red lines show the
root-mean-square of the difference of reconstructed and true values of the
radial velocity, with and without the IHM bias correction as a function of
radius. The dotted line shows the standard deviation in our reconstruction
samples. The dash-dotted line show the velocity scatter of the halos in the
VELMASS simulation. At large radius, the velocity error in our reconstruction
approaches the velocity scatter of the halos.

Having introduced our method in the previous section, we now present the
results of running our code on a mock peculiar velocity survey. We validate
our code by running the code on mock surveys produced from $N$-body
simulations and comparing the density and velocity reconstructions with the
ground truth.

Table 1: Depth of the three mock peculiar velocity surveys used in this work. $\Delta_{d}$ is the fractional uncertainty of the peculiar velocity tracers. $r_{95}$ is the radius of $95\%$ completeness for the mock survey and $r_{\text{peak}}$ is the radius of peak number density of the peculiar velocity tracers. | $\Delta_{d}$ |  $r_{95}~{}(h^{-1}$ Mpc) |  $r_{\text{peak}}~{}(h^{-1}$ Mpc)
---|---|---|---
Mock1 | $0.07$ | $135$ | $120$
Mock2 | $0.15$ | $143$ | $117$
Mock3 | $0.20$ | $176$ | $113$

We create our mock peculiar velocity surveys from the VELMASS $N$-body
simulation. We select dark matter haloes from the simulation with masses,
$M>3\times 10^{12}~{}\mathrm{M}_{\odot}$ to avoid using haloes with small
number of particles. We then assign an absolute magnitude, $\mathcal{M}$, to
each halo with an scatter, $\sigma_{\text{int}}$. We can estimate the distance
to these haloes from their apparent magnitude. Finally, we select the halos
brighter than an apparent magnitude, $m_{\text{cut}}$, thus resulting in a
magnitude-limited peculiar velocity survey. From the final magnitude-limited
sample, we randomly select $7000$ halos. In this paper, we report the results
with three different intrinsic scatters, corresponding to distance
uncertainties of $7\%$, $15\%$ and $20\%$. As we showed in section 3.2, the
IHM bias increases for a larger distance uncertainty. Therefore, we anticipate
the impact of the IHM bias in our reconstruction to be the larger for the
larger distance uncertainty. In Table 1, we show the depth of the three mock
surveys we used in our reconstruction showing the radius of $95\%$
completeness of the mock surveys and the peak of tracer number density. As
mentioned previously, our method requires an estimate of the galaxy density
field to correct for the inhomogeneous Malmquist bias. In our runs with the
simulations, we use the halo density field of the simulation (smoothed with a
Gaussian filter of of scale $4~{}h^{-1}$ Mpc) as the estimate of the galaxy
density field.

To compare the impact of the IHM bias correction, we also ran our method
without any IHM correction by setting $\hat{\delta}_{g}=0$. As we will see,
neglecting the IHM bias correction impacts our reconstruction significantly. A
visual comparison of the impact of the IHM correction in the velocity
reconstruction is shown in Figure 3, where we plot the reconstructed velocity
field and the associated uncertainty with our algorithm in the $Z=0$ plane. By
comparing the reconstructed mean field to the true velocity field in the
figure, we can see that our algorithm reconstructs the coherent large-scale
velocities well. However, we see that without the IHM bias correction, the
reconstructed velocity field has a higher amplitude compared to the true
velocity field. This is due to the fact that the IHM bias leads to a spurious
velocity in the reconstruction as discussed in section 3.2. As expected from
equation (2), this effect is most pronounced for a distance uncertainty of
$\Delta_{d}=0.20$. Finally, note that outside the data region, the mean
velocity field is suppressed to a smaller value of the velocity. The impact of
the IHM bias is most readily seen in the power spectrum of the reconstructed
density field. We show the power spectrum of our reconstruction (with and
without the IHM bias correction) in the top panel of Figure 4. We can see that
the reconstructed power spectrum without the IHM bias correction is biased
high. This can be seen more clearly in the middle panel of Figure 4 where we
plot the ratio of the mean of the inferred power spectrum to the theoretical
expectation. As we can see from the figure, the power spectrum of the
reconstruction with IHM bias correction is consistent with the true power
spectrum of the simulation. On the other hand, if we do not account for the
IHM bias, we see excess power in the reconstructions. This is due to the fact
that the IHM bias pushes the velocity to a larger value, which in turn
increases the variance in the density field. In order to quantify the bias in
the power spectrum, we calculate the quantity

$\Delta_{P}=\sum_{i=1}^{N_{k-\text{bins}}}\frac{\langle
P_{\textrm{sample}}(k_{i})\rangle-
P_{\text{theory}}(k_{i})}{\sigma[P(k_{i})]},$ (12)

where $\sigma[P(k_{i})]$ is the standard deviation in the power spectrum
assuming a Gaussian covariance. For an unbiased reconstruction, we expect this
quantity to have a value of zero. We compare the bias in the reconstruction by
computing the bias in the power spectrum measured in $15$ bins in $k$. For the
reconstruction with the IHM bias correction, we obtain values of
$\Delta_{P}=+0.46,-0.25,+0.55$ for the mock surveys with
$\Delta_{d}=0.07,0.15,0.20$ respectively. The same quantities when measured
for the reconstruction without the IHM bias correction are $\Delta_{P}=+8.78$,
$+15.24$ and $+16.99$ respectively showing that the power spectrum inferred
without IHM bias correction is biased high with high significance.

While the power spectrum gives a useful consistency check, the power spectrum
does not provide an assessment about the phases of the Fourier modes. In order
to test that the phases of the density modes are reconstructed correctly, we
compute the cross-correlation rate of the reconstructed density field with the
true underlying density field, which is defined as follows,

$\rho_{c}(k)=\frac{\langle\delta_{\text{rec}}({\bm{k}})\delta^{*}_{\text{sim}}({\bm{k}})\rangle_{{\bm{k}}\in
k~{}\text{bin}}}{V_{s}\sqrt{P_{\text{rec}}(k)P_{\text{sim}}(k)}},$ (13)

where $P_{\text{rec}}$ is the power spectrum of the reconstruction and
$P_{\text{sim}}$ is the true power spectrum of the simulation. $V_{s}$ is the
volume of the survey. We show the cross-correlation of our reconstructed
density samples with the true underlying density field of the simulation in
the bottom panels of Figure 4. As can be seen, the reconstructed large scale
density modes shows a cross-correlation of $\sim 0.5$ with the true density
field, showing that the large scale density modes are reconstructed well with
our method. The cross-correlation approaches zero for the small-scale
($k\gtrsim 0.04~{}h$ Mpc-1) modes. We can see the impact of the IHM bias on
the cross-correlation from this plot. For the mock survey with larger distance
uncertainty, for the reconstruction without IHM bias correction, cross-
correlations of the large-scale density modes are reduced compared to the
reconstruction with IHM bias correction. This effect is most readily seen in
the reconstruction with $\Delta_{d}=0.2$

We however note that, while the cross-correlation in Fourier space is an
interpretable metric, it may not be the best metric to compare the two fields.
This is because we only have the peculiar velocity data in a small volume of
our simulation box. Therefore, we also compare the radial velocity in our
reconstruction to the true radial velocity in shells of radial distances. In
each shell, we compute the slope, $m$ and the scatter, $\sigma_{v}$ between
the velocity estimates by regressing the reconstructed velocity on the true
velocity field. We use this slope as a proxy for the cross-correlation of the
two fields in real space. We show the result from this comparison in Figure 5.
As it can be seen from the figure, IHM bias corrected reconstruction has a
slope of $m\sim 1$ in the data region ($r\lesssim 120~{}h^{-1}$ Mpc) and then
decreases to a value of $m\sim 0$ at large $R$. On the other hand, the IHM
bias leads the reconstructed velocity to be biased high, which is reflected by
the fact that the slope, $m$ can be higher than $1$ for the reconstructions
with large IHM bias. Furthermore, we compare the reconstruction error for the
different runs of our method. The reconstruction error is defined as the Root
mean squared (RMS) error between the reconstruction and the true velocity,
$(\Delta V_{r})^{2}=|V^{\text{samples}}_{r}-V^{\text{true}}_{r}|^{2}$. We see
that the IHM bias corrected reconstruction leads to a smaller error in the
reconstructed velocity compared to the reconstruction without IHM bias
correction. At large radius, the velocity scatter between the reconstructed
velocity and the true velocity approaches the scatter of the halos in the
VELMASS simulation. We discuss in more detail the different sources
contributing to this reconstruction error in Appendix A.

The above tests show that we get an unbiased reconstruction of the velocity
field with our forward-modelled reconstruction method after correcting for the
IHM bias. On the other hand, not correcting for the inhomogeneous Malmquist
bias can substantially bias the velocity field reconstruction. We now move to
use our reconstruction method on real peculiar velocity surveys and
reconstruct the velocity field of the local Universe.

Figure 6: Power spectrum of the samples of reconstruction of the local
Universe. The reconstruction was performed using the 2MTF and the SFI++ Tully-
Fisher catalogues and the A2 supernovae compilation. The density
reconstruction from Carrick et al. (2015) was used for the IHM bias
correction. Figure 7: A comparison of our reconstruction with the
reconstruction of Carrick et al. (2015) in the galactic $X$ (left), $Y$
(centre) and $Z$ (right) planes. The top panel show the reconstruction of
Carrick et al. (2015), the middle and the bottom panels show the mean and the
standard deviation of the velocity reconstruction from our method. We can see
similar large scale features in the velocity field in both our reconstruction
and the C15 reconstruction. We show the radial distance of $100~{}h^{-1}$ Mpc
with a black dotted circle. This is the distance cut we use for our peculiar
velocity tracer sample. Note that the uncertainty in the velocity field
reconstruction increases drastically beyond this boundary.

##  5 Velocity reconstruction in the local Universe

Figure 8: The density and the velocity field reconstruction in the
supergalactic plane. The top plot shows the density field smoothed using a
Gaussian filter of smoothing scale $10~{}h^{-1}$ Mpc and the bottom panel
shows the velocity field reconstruction. On the left panel we show the density
and the velocity field of a randomly selected sample in our reconstruction.
The middle panels show the mean density and velocity field. The rightmost
panels show the standard deviation in our density and velocity field
reconstruction samples. We show the location of several prominent clusters in
the supergalactic plane, namely, Virgo, Norma, Perseus-Pisces (PP), Coma and
Shapley, with a black star. The black dotted line shows the radius at
$100~{}h^{-1}$ Mpc, which is the selection limit imposed on the sample. Figure
9: Comparison of (left) the slope and (right) the RMS difference between the
velocity estimates of our reconstruction and the C15 reconstruction for voxels
binned in radial shells. The slope can be interpreted as the cross-correlation
between the two velocity fields. We see a large degree of cross-correlation in
the region with the peculiar velocity data ($R\lesssim 100~{}h^{-1}$ Mpc),
while it rapidly goes to zero beyond this region. The right panel shows the
root-mean-squared error between the two velocity reconstructions. At very
small radius, $R\lesssim 50h^{-1}$ Mpc, the scatter approaches the value,
$\sigma_{\text{NL}}\approx 150$ km/s. The black dashed line show the standard
deviation in our reconstruction sample with $\sigma_{\text{NL}}$ added in
quadrature. Figure 10: Comparison of the bulk flow measured in our
reconstruction and compared to the $\Lambda$CDM expectation as well as the
other results in the literature. The left panel shows the magnitude of the
bulk flow, while the middle and the right panel shows the direction of the
bulk flow in terms of the galactic longitude and latitude respectively. As we
can see, the magnitude of the bulk flow in our reconstruction (shown in blue)
is consistent with both the $\Lambda$CDM expectation (shown in green hatch) as
well as other results in the literature. The direction of the bulk flow as
inferred from our reconstruction is also consistent with other results in the
literature. The results in the literature compared in this results are:
Scrimgeour et al. (2016) (S16), Qin et al. (2019) (Q17) Boruah et al. (2020)
(B20), Lilow & Nusser (2021) (L21)

In this section, we apply our method on real peculiar velocity data to
reconstruct the velocity field of the local Universe.

We use the SFI++ and 2MTF Tully-Fisher catalogues and the A2 supernovae data
set in this work. We use the galaxy density reconstruction of Carrick et al.
(2015) to correct for the inhomogeneous Malmquist bias. Since this
reconstruction is limited to $125~{}h^{-1}$ Mpc in certain sky directions, we
only consider peculiar velocity tracers within estimated distances,
$d<100h^{-1}$ Mpc in this work. Note that, this does not lead to selection
effects in the forward methods of peculiar velocity analysis (Strauss &
Willick, 1995). We marginalize over a scaling factor to account for
uncertainties in the zero point calibration. We use a different scaling factor
for each peculiar velocity survey. By sampling density fields from the
posterior (4), we get samples of reconstruction that are consistent with the
observed data. One useful side product of sampling from the posterior is that
we get an estimate of the correlated uncertainties in the reconstructed
velocity field. For the reconstruction of the local Universe, we use a cubic
sampling box of side length, $L_{\text{box}}=500~{}h^{-1}$ Mpc with $128^{3}$
grids. In Figure 6, we plot the power spectrum of the samples of our
reconstruction. Since we assume a Gaussian prior on the density field, we
expect the reconstruction to have a cosmological power spectrum. However, as
we saw in section 4, IHM bias may artificially enhance the inferred power
spectrum. Therefore, consistency between the theoretical power spectrum and
the inferred power spectrum shows that the IHM bias is corrected in our
reconstruction. As can be seen from the figure, this is indeed true and the
power spectrum computed from the samples of our posteriors agree very well
with the theoretical prediction from the $\Lambda$CDM model.

###  5.1 Comparison with 2M++ reconstruction

Unlike the simulations, for the real data, we do not have access to the true
velocity field. Nonetheless, we can compare our reconstruction to velocity
fields reconstructed using other methods. Therefore, we compare our
reconstructed velocities to the reconstruction of Carrick et al. (2015) who
used an iterative reconstruction procedure, where the galaxy density is
reconstructed from the galaxies in the 2M++ galaxy compilation (Lavaux &
Hudson, 2011). The velocity is predicted from the galaxy density field using
linear perturbation theory. However, since galaxies are a biased tracer of the
underlying density field, the velocity estimate from the galaxy field needs to
be scaled by a factor of $\beta=f/b$. The value of $\beta$ is then fitted by
comparing the predicted velocities to the the peculiar velocity data from a
peculiar velocity survey. In this work, we use the best fit $\beta$ and
${\bm{V}}_{\text{ext}}$ value fitted in Carrick et al. (2015). We show a
visual comparison of our velocity reconstruction with the reconstruction of
Carrick et al. (2015) (hereafter C15) in the galactic $X,Y,Z$ planes in Figure
7. The top panels of Figure 7 shows the radial velocity field in the
reconstruction of C15. The middle and the bottom panels show the mean and the
standard deviation of the velocity field samples reconstructed using our
algorithm. By comparing our velocity reconstruction to the 2M++
reconstruction, we see the similarity between some of the same large scale
features in both reconstruction. Note that in the region without data, the
mean field approaches very small values. This is similar to the Weiner
filtered fields, where in the region without data, the field is strongly
suppressed such that in the limit of no data, the velocity field approaches
zero. In the bottom panel, we plot the uncertainty in the velocity estimates
in these planes. Since our method produces samples of reconstruction from the
posterior with the Gaussian prior, we get the correlated uncertainties in
velocity estimates of our reconstruction. From the figure, we see that the
uncertainties in the velocity estimates is substantially lower in the inner
region which has the peculiar velocity data than the outer region which does
not have peculiar velocity data. We also plot the velocity field and the
density field in the supergalactic plane comparing also the position of some
prominent galaxy clusters in Figure 8. From the mean field, we see that while
the over-densities around the Perseus-Pisces, Virgo and Coma clusters are
well-reconstructed, we do not see prominent overdensities around the Norma and
the Shapley clusters. Since the Shapley cluster lies outside the
$100~{}h^{-1}$ Mpc, the threshold of our reconstruction and the Norma cluster
lies along the galactic plane, there is lack of peculiar velocity data around
these clusters. The lack of peculiar velocity data around the cluster leads to
poor reconstruction quality around these regions.

Finally, we compare our reconstruction to the C15 reconstruction by comparing
the velocity fields in radial shells. For this comparison, we compared the
radial velocity from C15 and our reconstruction in radial shells of width
$20~{}h^{-1}$ Mpc. In each shell, we compare the radial velocities from our
reconstruction samples and the reconstruction from C15. We fit the slope and
the standard deviation between the two velocity estimates by regressing the
velocity estimates of our reconstruction on the C15 velocity estimates. The
fitted slope can be thought of as a proxy for the cross-correlation between
the two velocity estimates. The slope and the scatter between the two velocity
fields are shown in Figure 9. As we can see from the figure, the two velocity
fields have a large degree of correlation in the nearby Universe ($R\lesssim
120~{}h^{-1}$ Mpc). Since the peculiar velocity data is restricted to
$d<100h^{-1}$ Mpc, we do not expect the reconstructed velocity field to be a
good estimate beyond this boundary. This is reflected in the fact that the
slope rapidly approaches zero beyond the peculiar velocity data boundary. We
also see that the scatter between the two increases with increasing radius. At
very low radius ($R\lesssim 50~{}h^{-1}$ Mpc), the scatter between the two
velocity fields is less than or comparable to the uncertainty due to the non-
linearities, $\sigma_{\text{NL}}\approx 150$ km/s. At large radius, the
scatter approaches the sample standard deviation ($\sim 250$ km/s).

###  5.2 Bulk flow

Table 2: Comparison of the bulk flow results with other studies in the literature. The results are also shown in Figure 10. We see that our bulk flow result is in excellent agreement with the other results in literature. Reference | Abbreviation | Effective radius |  $|{\bm{V}}_{\textrm{bulk}}|$ (km/s) |  $l$ (deg) |  $b$ (deg)
---|---|---|---|---|---
Scrimgeour et al. (2016) | S16 |  $50\ h^{-1}$ Mpc | $248\pm 58$ | $318\pm 20$ | $40\pm 13$
Qin et al. (2019) | Q17 |  $37\ h^{-1}$ Mpc | $259\pm 15$ | $300\pm 4$ | $23\pm 3$
Boruah et al. (2020) | B20 |  $40\ h^{-1}$ Mpc | $252\pm 11$ | $293\pm 5$ | $14\pm 5$
Lilow & Nusser (2021) | L21 |  $50\ h^{-1}$ Mpc | $274\pm 50$ | $287\pm 9$ | $11\pm 10$
This work | — |  $\bm{40\ h^{-1}}$ Mpc | $\bm{220\pm 21}$ | $\bm{295\pm 6}$ | $\bm{21\pm 5}$

Due to the sensitivity to large scale density modes, the velocity field is
correlated on very large scales. The large scale velocity field is often
expanded in terms of its kinematic moments, such as the dipole and the
quadrupole of the velocity field (Jaffe & Kaiser, 1995),

${\bm{V}}(d)={\bm{U}}+d{\mathbf{Q}}\cdot\hat{{\bm{r}}}+...,$ (14)

where, ${\bm{U}}$ is the dipolar bulk flow and ${\mathbf{Q}}$ is the trace-
free, symmetric quadrupolar shear moment and $d$ is the radial distance. The
bulk flow is sensitive to the large scale power of the density field and
therefore can also be used to constrain cosmological parameters due to its
sensitivity to the power spectrum shape parameter (Feldman & Watkins, 2008).

We can measure the bulk flow and the shear moments directly from our
reconstruction. In order to measure the bulk flow in our reconstructions at a
given scale, we smooth the velocity field with a window function of
characteristic scale, $R_{\text{bulk}}$. We use a Gaussian filter in order to
measure the bulk flow in our reconstruction. The bulk flow at a scale of
$R_{\text{bulk}}=40h^{-1}$ Mpc for our reconstruction is
$V_{\text{bulk}}=220^{+20}_{-21}$ km/s in the direction $l=295^{\circ}\pm
6^{\circ}$, $b=21^{\circ}\pm 5^{\circ}$ degree. The bulk flow as measured in
our reconstruction is compared to a number of different results in the
literature Table 2. We also show the total magnitude and the direction of the
bulk flow as measured from our reconstruction samples in Figure 10. As we can
see from the figure, both the magnitude and the direction of the bulk flow is
consistent with other results in the literature. Of the compared results,
Scrimgeour et al. (2016) and Qin et al. (2019) use an estimator to estimate
the bulk flow directly from the peculiar velocity data. On the other hand,
Boruah et al. (2020) and Lilow & Nusser (2021) explicitly fit an external
dipole in their flow model to account for velocity contribution from outside
the survey volume. Note that the approach presented here is different from
these approaches. We do not explicitly model an additional dipole. The density
modes we fit for in our reconstruction not only accounts for the bulk flow,
but also other higher kinematic moments of the large scale velocity field.

###  5.3 Velocity field comparison

Table 3: Comparison of the logarithm of Bayes factor for various redshift selection. A positive log of Bayes factor implies that our model is favored over the compared model, while a negative log for the Bayes factor implies the compared model is preferred over our model. From the table, we can see that our model performs better than the adaptive kernel smoothed velocity field. On the other hand C15 velocity field performs better than our forward modelled velocity field. Test set | Selection | $\ln\bigg{(}\frac{\mathcal{P}_{\text{\text{fwd-PV}}}}{\mathcal{P}^{\text{TF}}_{\text{adaptive}}}\bigg{)}$ | $\ln\bigg{(}\frac{\mathcal{P}_{\text{fwd-PV}}}{\mathcal{P}_{\text{2M++}}}\bigg{)}$ | $N_{\text{tracers}}$
---|---|---|---|---
A2 | $cz<0.01$ | $7.44$ | $-2.12$ | $49$
| $z<0.015$ | $10.88$ | $-5.92$ | $92$
| $z<0.02$ | $21.73$ | $-20.69$ | $168$
| $z<0.03$ | $41.16$ | $-27.65$ | $310$

We introduced a Bayesian model comparison framework in Boruah et al. (2021) to
compare the performance of different velocity reconstruction models. In this
model comparison framework, we look at the Bayes factor between two models,
$\mathcal{M}_{1}$ and $\mathcal{M}_{2}$,

$\text{Bayes
factor}=\frac{\mathcal{P}(\mathcal{D}|\mathcal{M}_{1})}{\mathcal{P}(\mathcal{D}|\mathcal{M}_{2})}.$
(15)

If the Bayes factor is greater than 1, the model $\mathcal{M}_{1}$ is
preferred over the model $\mathcal{M}_{2}$ and vice-versa. We use this model
comparison framework to assess the quality of the velocity reconstruction
method introduced in this paper. We compare our velocity model to two
different velocity reconstruction of the local Universe - i) The
reconstruction of C15, ii) an adaptive kernel-smoothed velocity
reconstruction, where the peculiar velocity data is directly smoothed
(Springob et al., 2014, 2016). We use an adaptive-kernel smoothed velocity
field obtained from a combined Tully-Fisher catalog consisting of the SFI++
and the 2MTF velocity field. The Bayesian model comparison method relies on
our capacity to predict new data points from the model. This is called a
posterior predictive test, and thus requires what is sometimes called a “test
set” of peculiar velocities. Therefore for this comparison, we perform a
reconstruction using only the Tully-Fisher catalogs for the reconstruction and
use the A2 sample as a test set for the posterior predictive analysis. We note
that this is the same peculiar velocity data set used for the adaptive kernel-
smoothed reconstruction.

The results of the model comparison is shown in Table 3. Using our model
comparison comparison, we find that the forward-modelled velocity
reconstruction performs better than the kernel smoothed velocity field, even
though both reconstructions use the same data sets. Our velocity field however
does not perform as well at predicting new velocity data than the
reconstruction of Carrick et al. (2015), which uses the 2M++ galaxies to map
the density field, and linear perturbation theory to predict peculiar
velocities. This is not a failure of our method, but rather a limitation of
the peculiar velocity data used for the reconstruction. Since the peculiar
velocity data used in this work are very noisy and sparse, the quality of
reconstruction is not comparable to the 2M++ reconstruction. In the future,
with denser peculiar velocity samples, we expect the quality of reconstruction
to be significantly better. There is also the possibility that the forward
linear modeling approximation may fail which would require more complex models
based, e.g., on borg.

##  6 Discussion

Our work provides a way towards including peculiar velocity data into the
initial condition reconstruction framework, borg. A density reconstruction by
combining both galaxy and peculiar velocity data in borg potentially has
multiple advantages. First, using borg, we can access the fully non-linear
velocity field, which is accurate to much smaller scales as compared to the
velocity field from linear perturbation theory. For example, borg particle
mesh runs has been shown to perform better than linear velocity fields in
$N$-body simulations (Mukherjee et al., 2021). Furthermore, the vorticity of
the velocity field, which is a uniquely non-linear phenomenon was already
reconstructed in Jasche & Lavaux (2019). Second, adding the velocity data to
the forward-modelled reconstruction can potentially make the reconstruction
robust to the details of the galaxy bias model. As was shown in Nguyen et al.
(2021), the amplitude of the inferred power spectrum in borg reconstruction
with galaxy/halo catalogues alone may be substantially biased. This arises
because the higher order bias terms are not sufficiently informative to break
the degeneracy between linear galaxy bias and the amplitude of matter
clustering, $\sigma_{8}$. Since peculiar velocity data is sensitive to the
total matter density, the addition of peculiar velocity data to galaxy data
would make the borg reconstruction robust to the details of the galaxy bias.
Adding velocity information to the forward-modelled reconstruction can
potentially lead to improved reconstruction of the full phase space structure
(density + velocity) of dark matter (Leclercq et al., 2017). Such phase space
reconstruction can then potentially guide observation efforts for new
discoveries (Kostić et al., 2021).

In the future, we can also extend our method to simultaneously constrain
cosmological parameters from peculiar velocity data. As we highlight in this
paper, not accounting for the IHM bias can lead to the amplitude of the
peculiar velocities and so the inferred power spectrum of the density/velocity
field to be biased high. This is crucial for the inference of the cosmological
parameters, since this bias will likely translate to the biased inference of
cosmological parameters. Thus, this paper provides the foundation of inferring
cosmological parameters with peculiar velocity data in a forward modelled
framework. Also note that our method is sensitive to the large-scale flows,
which in turn are a sensitive probe of the break in the matter power spectrum
(e.g, Feldman & Watkins, 2008). Thus, the large-scale power spectrum can
constrain the shape parameter, $\Gamma\approx\Omega_{m}h$. Combining these
constraints from large-scale flows with traditional velocity-velocity
comparison method can potentially provide additional cosmological information.
The use of the non-linear velocity field using borg can also potentially break
the $f\sigma_{8}$ degeneracy through the non-linearities. This feature was
already noted in the context of weak lensing by using a forward model based on
Lagrangian Perturbation theory within the borg framework for mock weak lensing
data (Porqueres et al., 2021a).

Another possible application of our forward-modelled velocity reconstructions
is in using the kinetic Sunyaev-Zeldovich (kSZ) data to reconstruct the
largest scale velocity modes. kSZ-based velocity reconstruction has been shown
to improve the constraints on primordial non-Gaussianities through multi-
tracer analyses (Giri & Smith, 2020).

##  7 Summary

In this paper, we introduced a forward-modelled velocity reconstruction method
which uses peculiar velocity data to reconstruct the velocity field. Using an
external estimate of the line of sight density, we consistently correct for
the inhomogeneous Malmquist (IHM) bias. Using mocks created from an $N$-body
simulation, we validated that our method leads to unbiased velocity field
reconstruction after accounting for the inhomogeneous Malmquist bias. However,
as we show using mock simulations, not accounting for the inhomogeneous
Malmquist bias in forward-modelled methods can lead to significantly biased
reconstruction. Since the IHM bias may induce spurious flows, without the
correction for the IHM bias, the inferred power spectrum is biased high. We
then applied our method to the 2MTF, SFI++ Tully-Fisher catalogues and the A2
supernovae compilation, resulting in a novel forward-modelled velocity field
reconstruction of the local Universe. Since we sample our velocity field
reconstruction samples from the field-level posterior, we get an estimate of
the full correlated uncertainties in the peculiar velocity estimates.
Furthermore, we found that the bulk flow calculated from our reconstruction is
consistent with other results in the literature. Using a Bayesian model
comparison framework, we showed that the reconstructed forward modelled
velocity fields perform better than the widely used adaptive kernel-smoothed
velocity fields constructed from the same data. However, the reconstruction
presented here does not perform as well as the velocity field reconstructed
from the 2M++ galaxy catalogue. This is likely because of peculiar velocity
data sets are sparser and noisier compared to the galaxy catalogues. Our
method provides a way to extend forward-modelled initial condition
reconstruction algorithms such as borg by including peculiar velocity data
sets.

## Acknowledgement

We thank Eleni Tsaprazi, Florent Leclercq and James Prideaux-Ghee for useful
feedback on the manuscript. This work has been done as part of the activities
of the Domaine d’Intérêt Majeur (DIM) “Astrophysique et Conditions
d’Apparition de la Vie” (ACAV), and received financial support from Région
Ile-de-France. GL acknowledges financial support from the ANR BIG4, under
reference ANR-16-CE23-0002. MH acknowledges the support of an NSERC Discovery
Grant. This work is done within the Aquila Consortium222https://www.aquila-
consortium.org/. This work has made use of the Infinity Cluster hosted by
Institut d’Astrophysique de Paris. We thank Stephane Rouberol for running this
cluster smoothly for us.

## Data availability statement

The Second Amendment (A2) supernovae compilation is available with the
supplementary data of https://doi.org/10.1093/mnras/staa2485. The 2MTF and the
SFI++ catalogues are publicly available with their respective publications
that are cited in section 2. The 2M++ reconstruction used in this work is
publicly available at https://cosmicflows.iap.fr/. The data products generated
in this work will be shared on reasonable request to the authors.

## References

  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109
  * Bertschinger & Dekel (1989) Bertschinger E., Dekel A., 1989, ApJ, 336, L5
  * Bertschinger et al. (1990) Bertschinger E., Dekel A., Faber S. M., Dressler A., Burstein D., 1990, ApJ, 364, 370
  * Boruah et al. (2020) Boruah S. S., Hudson M. J., Lavaux G., 2020, MNRAS, 498, 2703
  * Boruah et al. (2021) Boruah S. S., Hudson M. J., Lavaux G., 2021, MNRAS, 507, 2697
  * Bouchet et al. (1995) Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995, A&A, 296, 575
  * Bradbury et al. (2018) Bradbury J., Frostig R., Hawkins P., Johnson M. J., Leary C., Maclaurin D., Wanderman-Milne S., 2018, JAX: composable transformations of Python+NumPy programs, http://github.com/google/jax
  * Carrick et al. (2015) Carrick J., Turnbull S. J., Lavaux G., Hudson M. J., 2015, MNRAS, 450, 317
  * Feldman & Watkins (2008) Feldman H. A., Watkins R., 2008, MNRAS, 387, 825
  * Fiedorowicz et al. (2021) Fiedorowicz P., Rozo E., Boruah S. S., Chang C., Gatti M., 2021, arXiv e-prints, p. arXiv:2105.14699
  * Foley et al. (2018) Foley R. J., et al., 2018, MNRAS, 475, 193
  * Ganeshalingam et al. (2013) Ganeshalingam M., Li W., Filippenko A. V., 2013, MNRAS, 433, 2240
  * Giri & Smith (2020) Giri U., Smith K. M., 2020, arXiv e-prints, p. arXiv:2010.07193
  * Graziani et al. (2019) Graziani R., Courtois H. M., Lavaux G., Hoffman Y., Tully R. B., Copin Y., Pomarède D., 2019, MNRAS, 488, 5438
  * Hicken et al. (2009) Hicken M., Wood-Vasey W. M., Blondin S., Challis P., Jha S., Kelly P. L., Rest A., Kirshner R. P., 2009, ApJ, 700, 1097
  * Hoffman & Ribak (1991) Hoffman Y., Ribak E., 1991, ApJ, 380, L5
  * Hoffman et al. (2015) Hoffman Y., Courtois H. M., Tully R. B., 2015, MNRAS, 449, 4494
  * Hoffman et al. (2018) Hoffman Y., et al., 2018, Nature Astronomy, 2, 680
  * Hollinger & Hudson (2021) Hollinger A. M., Hudson M. J., 2021, MNRAS, 502, 3723
  * Hong et al. (2019) Hong T., et al., 2019, MNRAS, 487, 2061
  * Hudson (1994) Hudson M. J., 1994, MNRAS, 266, 468
  * Jaffe & Kaiser (1995) Jaffe A. H., Kaiser N., 1995, ApJ, 455, 26
  * Jasche & Lavaux (2019) Jasche J., Lavaux G., 2019, A&A, 625, A64
  * Jasche & Wandelt (2013) Jasche J., Wandelt B. D., 2013, MNRAS, 432, 894
  * Jones et al. (2019) Jones D. O., et al., 2019, ApJ, 881, 19
  * Kodi Ramanah et al. (2019) Kodi Ramanah D., Charnock T., Lavaux G., 2019, Phys. Rev. D, 100, 043515
  * Kostić et al. (2021) Kostić A., Jasche J., Kodi Ramanah D., Lavaux G., 2021, arXiv e-prints, p. arXiv:2107.00657
  * Krisciunas et al. (2017) Krisciunas K., et al., 2017, AJ, 154, 211
  * Lavaux (2016) Lavaux G., 2016, MNRAS, 457, 172
  * Lavaux & Hudson (2011) Lavaux G., Hudson M. J., 2011, MNRAS, 416, 2840
  * Leclercq et al. (2017) Leclercq F., Jasche J., Lavaux G., Wandelt B., Percival W., 2017, J. Cosmology Astropart. Phys., 2017, 049
  * Lilow & Nusser (2021) Lilow R., Nusser A., 2021, MNRAS, 507, 1557
  * Lynden-Bell et al. (1988) Lynden-Bell D., Faber S. M., Burstein D., Davies R. L., Dressler A., Terlevich R. J., Wegner G., 1988, ApJ, 326, 19
  * Masters et al. (2006) Masters K. L., Springob C. M., Haynes M. P., Giovanelli R., 2006, ApJ, 653, 861
  * Masters et al. (2008) Masters K. L., Springob C. M., Huchra J. P., 2008, AJ, 135, 1738
  * Millea et al. (2020a) Millea M., et al., 2020a, arXiv e-prints, p. arXiv:2012.01709
  * Millea et al. (2020b) Millea M., Anderes E., Wandelt B. D., 2020b, Phys. Rev. D, 102, 123542
  * Modi et al. (2018) Modi C., Feng Y., Seljak U., 2018, J. Cosmology Astropart. Phys., 2018, 028
  * Mukherjee et al. (2021) Mukherjee S., Lavaux G., Bouchet F. R., Jasche J., Wandelt B. D., Nissanke S., Leclercq F., Hotokezaka K., 2021, A&A, 646, A65
  * Neal (1993) Neal R., 1993, Technical Report CRG-TR-93-1. Dept. of Computer Science, University of Toronto.
  * Neal (1996) Neal R., 1996, Lecture Notes in Statistics, No. 118. Springer-Verlag
  * Nguyen et al. (2021) Nguyen N.-M., Schmidt F., Lavaux G., Jasche J., 2021, J. Cosmology Astropart. Phys., 2021, 058
  * Peebles (1980) Peebles P. J. E., 1980, The large-scale structure of the universe. Princeton University Press
  * Pike & Hudson (2005) Pike R. W., Hudson M. J., 2005, ApJ, 635, 11
  * Porqueres et al. (2021a) Porqueres N., Heavens A., Mortlock D., Lavaux G., 2021a, MNRAS,
  * Porqueres et al. (2021b) Porqueres N., Heavens A., Mortlock D., Lavaux G., 2021b, MNRAS, 502, 3035
  * Qin et al. (2019) Qin F., Howlett C., Staveley-Smith L., Hong T., 2019, MNRAS, 482, 1920
  * Scrimgeour et al. (2016) Scrimgeour M. I., et al., 2016, MNRAS, 455, 386
  * Springob et al. (2007) Springob C. M., Masters K. L., Haynes M. P., Giovanelli R., Marinoni C., 2007, ApJS, 172, 599
  * Springob et al. (2014) Springob C. M., et al., 2014, MNRAS, 445, 2677
  * Springob et al. (2016) Springob C. M., et al., 2016, MNRAS, 456, 1886
  * Strauss & Willick (1995) Strauss M. A., Willick J. A., 1995, Phys. Rep., 261, 271
  * Tassev et al. (2013) Tassev S., Zaldarriaga M., Eisenstein D. J., 2013, J. Cosmology Astropart. Phys., 2013, 036
  * Taylor et al. (2008) Taylor J. F., Ashdown M. A. J., Hobson M. P., 2008, MNRAS, 389, 1284
  * Tully et al. (2016) Tully R. B., Courtois H. M., Sorce J. G., 2016, AJ, 152, 50
  * Verde et al. (2019) Verde L., Treu T., Riess A. G., 2019, Nature Astronomy, 3, 891
  * Willick et al. (1997) Willick J. A., Strauss M. A., Dekel A., Kolatt T., 1997, ApJ, 486, 629
  * Zaroubi (2002) Zaroubi S., 2002, MNRAS, 331, 901
  * Zaroubi et al. (1999) Zaroubi S., Hoffman Y., Dekel A., 1999, ApJ, 520, 413

##  Appendix A Physical origin of the reconstruction error

In this appendix, we try to isolate different physical origins of the
uncertainties in our velocity reconstruction. As shown in section 4, without
the use of the correct line-of-sight (LOS) distribution for the peculiar
velocity tracers, the reconstructed velocity field is biased high. This bias
also leads to increased uncertainty in the reconstructed velocity. However,
another factor contributes to the increased error in the velocity
reconstruction. In equation (11), there are two contributions to the
distribution of peculiar velocity tracers along the LOS - i) The uncertainty
in the distance, ii) The line of sight density field. Due to the latter
contribution, peculiar velocity tracers are clustered around regions of high
density. Therefore, when we take the LOS density into account, the radial
distribution may be more peaked than naively expected from the measured
distance uncertainty. This is shown in Figure 11, where we plot an example
from our mocks used in section 4. In this figure, the radial distribution of
the tracer including the LOS density is much sharply peaked compared to the
distribution is we ignored the LOS.

Therefore, when we account for the LOS density, we not only correct for the
IHM bias, but we also sharpen the distance estimate and therefore decrease the
uncertainty in the reconstruction. In order to investigate the impact of this
factor, we use our method on two different mock surveys. In mock survey A, we
populate the tracers according to equation (11), accounting for the LOS
density field. In mock B, we populate tracers homogeneously, i.e, we do not
account for the LOS inhomogeneity. Therefore, by design, inference with mock B
does not suffer from IHM bias. Since the velocity tracers are not clustered
around density peaks, we are not able to determine the position of the tracer
better than the distance uncertainty. We then make three runs of our code on
the two mocks:

  * •

Run 1: Mock A, corrected for IHM bias.

  * •

Run 2: Mock A, not corrected for IHM bias.

  * •

Run 3: Mock B, not corrected for IHM bias (by design, Mock B is unaffected by
IHM bias).

The results of the three runs are shown in Figure 12, where we plot the Root
Mean Squared (RMS) error in the velocity reconstruction binned in radial bins.
We see that the RMS error for the Run 1 is substantially lower than Run 2 and
3. Furthermore the RMS error on Run2 is higher than Run3. Note that only Run2
is impacted by the IHM bias. We see that the RMS error in Run3 is closer to
Run2 than to Run1, thus suggesting that extra error of the IHM uncorrected run
over the IHM corrected run is driven by the fact that the distance estimates
are better localized when we account for the LOS density field, rather than
the IHM bias. Thus, by accounting for the LOS density, we not only correct for
the IHM bias, but also reduce the reconstruction error by better localization
of the peculiar velocity tracers.

Figure 11: The impact of the line-of-sight density field on the radial
distribution of the peculiar velocity tracer. Here we show the radial
distribution of a peculiar velocity tracer in our mock survey. The red curve
shows the expected radial distribution ignoring the line-of-sight
inhomogeneities and the blue curve shows the distribution after accounting for
the LOS inhomogeneities. As we can see, radial distribution after accounting
for the line-of-sight inhomogeneities is more peaked than the naive
distribution expected from the measured distance error. Figure 12: The RMS
error in the velocity reconstructions in our three runs. The peculiar velocity
tracers in Mock A are sampled by accounting for the LOS inhomogeneities and
therefore the reconstruction is susceptible to IHM bias. On the other hand, in
Mock B, the peculiar velocity tracers are sampled homogeneously and therefore
does not suffer from IHM bias. We notice that the RMS error in Run 1 is much
lower than Run 2 and 3. On the other hand, the RMS error in Run 3 is closer to
the error in Run 2 (compared to Run1). This suggests that the error in the IHM
corrected reconstruction has the additional advantage of reducing the
reconstruction error by better localization of the peculiar velocity tracers.


* Sung  _et al._ (2021) Y. Sung, L. Ding, J. Braumüller, A. Vepsäläinen, B. Kannan, M. Kjaergaard, A. Greene, G. O. Samach, C. McNally, D. Kim, A. Melville, B. M. Niedzielski, M. E. Schwartz, J. L. Yoder, T. P. Orlando, S. Gustavsson, and W. D. Oliver, “Realization of High-Fidelity CZ and ZZ-Free iSWAP Gates with a Tunable Coupler,” Physical Review X 11, 021058 (2021), arXiv:2011.01261 .
  * Foxen  _et al._ (2020) B. Foxen _et al._ , “Demonstrating a Continuous Set of Two-qubit Gates for Near-term Quantum Algorithms,” Physical Review Letters 125, 120504 (2020), arXiv:2001.08343 .
  * Xu  _et al._ (2021) H. Xu, W. Liu, Z. Li, J. Han, J. Zhang, K. Linghu, Y. Li, M. Chen, Z. Yang, J. Wang, T. Ma, G. Xue, Y. Jin, and H. Yu, “Realization of adiabatic and diabatic CZ gates in superconducting qubits coupled with a tunable coupler,” Chinese Physics B 30, 044212 (2021), arXiv:2010.14053 .
  * McKay  _et al._ (2016) D. C. McKay, S. Filipp, A. Mezzacapo, E. Magesan, J. M. Chow, and J. M. Gambetta, “Universal Gate for Fixed-Frequency Qubits via a Tunable Bus,” Physical Review Applied 6, 064007 (2016), arXiv:1604.03076 .
  * Roth  _et al._ (2017) M. Roth, M. Ganzhorn, N. Moll, S. Filipp, G. Salis, and S. Schmidt, “Analysis of a parametrically driven exchange-type gate and a two-photon excitation gate between superconducting qubits,” Physical Review A 96, 062323 (2017), arXiv:1708.02090 .
  * Bengtsson  _et al._ (2020) A. Bengtsson, P. Vikstål, C. Warren, M. Svensson, X. Gu, A. F. Kockum, P. Krantz, C. Križan, D. Shiri, I.-M. Svensson, G. Tancredi, G. Johansson, P. Delsing, G. Ferrini, and J. Bylander, “Improved Success Probability with Greater Circuit Depth for the Quantum Approximate Optimization Algorithm,” Physical Review Applied 14, 034010 (2020), arXiv:1912.10495 .
  * Ganzhorn  _et al._ (2020) M. Ganzhorn, G. Salis, D. J. Egger, A. Fuhrer, M. Mergenthaler, C. Müller, P. Müller, S. Paredes, M. Pechal, M. Werninghaus, and S. Filipp, “Benchmarking the noise sensitivity of different parametric two-qubit gates in a single superconducting quantum computing platform,” Physical Review Research 2, 033447 (2020), arXiv:2005.05696 .
  * Sjöqvist (2016) E. Sjöqvist, “Nonadiabatic holonomic single-qubit gates in off-resonant $\Lambda$ systems,” Physics Letters A 380, 65 (2016), arXiv:1511.00911 .
  * Sjöqvist  _et al._ (2012) E. Sjöqvist, D. M. Tong, L. M. Andersson, B. Hessmo, M. Johansson, and K. Singh, “Non-adiabatic holonomic quantum computation,” New Journal of Physics 14, 103035 (2012), arXiv:1107.5127 .
  * Fredkin and Toffoli (1982) E. Fredkin and T. Toffoli, “Conservative logic,” International Journal of Theoretical Physics 21, 219 (1982).
  * Milburn (1989) G. J. Milburn, “Quantum optical Fredkin gate,” Physical Review Letters 62, 2124 (1989).
  * Abdumalikov Jr  _et al._ (2013) A. A. Abdumalikov Jr, J. M. Fink, K. Juliusson, M. Pechal, S. Berger, A. Wallraff, S. Filipp, A. A. Abdumalikov, J. M. Fink, K. Juliusson, M. Pechal, S. Berger, A. Wallraff, and S. Filipp, “Experimental realization of non-Abelian non-adiabatic geometric gates,” Nature 496, 482 (2013), arXiv:1304.5186 .
  * Danilin  _et al._ (2018) S. Danilin, A. Vepsäläinen, and G. S. Paraoanu, “Experimental state control by fast non-Abelian holonomic gates with a superconducting qutrit,” Physica Scripta 93, 055101 (2018), arXiv:1804.01759 .
  * Egger  _et al._ (2019) D. J. Egger, M. Ganzhorn, G. Salis, A. Fuhrer, P. Müller, P. K. Barkoutsos, N. Moll, I. Tavernelli, and S. Filipp, “Entanglement Generation in Superconducting Qubits Using Holonomic Operations,” Physical Review Applied 11, 014017 (2019), arXiv:1804.04900 .
  * Ganzhorn  _et al._ (2019) M. Ganzhorn, D. J. Egger, P. Barkoutsos, P. Ollitrault, G. Salis, N. Moll, M. Roth, A. Fuhrer, P. Mueller, S. Woerner, I. Tavernelli, and S. Filipp, “Gate-Efficient Simulation of Molecular Eigenstates on a Quantum Computer,” Physical Review Applied 11, 044092 (2019), arXiv:1809.05057 .
  * Abrams  _et al._ (2020) D. M. Abrams, N. Didier, B. R. Johnson, M. P. da Silva, and C. A. Ryan, “Implementation of XY entangling gates with a single calibrated pulse,” Nature Electronics 3, 744 (2020), arXiv:1912.04424 .
  * Toffoli (1980) T. Toffoli, “Reversible computing,” in  _Automata, Languages and Programming_, Lecture Notes in Computer Science, Vol. 85, edited by J. de Bakker and J. van Leeuwen (Springer Berlin Heidelberg, Berlin, Heidelberg, 1980).
  * Fedorov  _et al._ (2012) A. Fedorov, L. Steffen, M. Baur, M. P. da Silva, and A. Wallraff, “Implementation of a Toffoli gate with superconducting circuits,” Nature 481, 170 (2012), arXiv:1108.3966 .
  * Horodecki  _et al._ (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, “Quantum entanglement,” Reviews of Modern Physics 81, 865 (2009), arXiv:0702225 [quant-ph] .
  * Briegel  _et al._ (2009) H. J. Briegel, D. E. Browne, W. Dür, R. Raussendorf, and M. Van den Nest, “Measurement-based quantum computation,” Nature Physics 5, 19 (2009), arXiv:0910.1116 .
  * Preskill (2012) J. Preskill, “Quantum computing and the entanglement frontier,”  (2012), arXiv:1203.5813 .
  * Gisin and Thew (2007) N. Gisin and R. Thew, “Quantum communication,” Nature Photonics 1, 165 (2007), arXiv:0703255 [quant-ph] .
  * Wang  _et al._ (2018) Y. Wang, Y. Li, Z.-Q. Yin, and B. Zeng, “16-qubit IBM universal quantum computer can be fully entangled,” npj Quantum Information 4, 46 (2018), arXiv:1801.03782 .
  * Mooney  _et al._ (2019) G. J. Mooney, C. D. Hill, and L. C. L. Hollenberg, “Entanglement in a 20-Qubit Superconducting Quantum Computer,” Scientific Reports 9, 13465 (2019), arXiv:1903.11747 .
  * Wei  _et al._ (2020) K. X. Wei, I. Lauer, S. Srinivasan, N. Sundaresan, D. T. McClure, D. Toyli, D. C. McKay, J. M. Gambetta, and S. Sheldon, “Verifying multipartite entangled Greenberger-Horne-Zeilinger states via multiple quantum coherences,” Physical Review A 101, 032343 (2020), arXiv:1905.05720 .
  * Mooney  _et al._ (2021b) G. J. Mooney, G. A. L. White, C. D. Hill, and L. C. L. Hollenberg, “Whole-device entanglement in a 65-qubit superconducting quantum computer,”  (2021b), arXiv:2102.11521 .
  * Yang  _et al._ (2021) B. Yang, R. Raymond, H. Imai, H. Chang, and H. Hiraishi, “Testing Scalable Bell Inequalities for Quantum Graph States on IBM Quantum Devices,”  (2021), arXiv:2101.10307 .
  * Friis  _et al._ (2018) N. Friis, O. Marty, C. Maier, C. Hempel, M. Holzäpfel, P. Jurcevic, M. B. Plenio, M. Huber, C. Roos, R. Blatt, and B. Lanyon, “Observation of Entangled States of a Fully Controlled 20-Qubit System,” Physical Review X 8, 021012 (2018), arXiv:1711.11092 .
  * Pogorelov  _et al._ (2021) I. Pogorelov, T. Feldker, C. D. Marciniak, L. Postler, G. Jacob, O. Kriegelsteiner, V. Podlesnic, M. Meth, V. Negnevitsky, M. Stadler, B. Höfer, C. Wächter, K. Lakhmanskiy, R. Blatt, P. Schindler, and T. Monz, “A compact ion-trap quantum computing demonstrator,”  (2021), arXiv:2101.11390 .
  * Omran  _et al._ (2019) A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Montangero, T. Calarco, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, “Generation and manipulation of Schrödinger cat states in Rydberg atom arrays,” Science 365, 570 (2019), arXiv:1905.05721 .
  * Neeley  _et al._ (2010) M. Neeley, R. C. Bialczak, M. Lenander, E. Lucero, M. Mariantoni, A. D. O’Connell, D. Sank, H. Wang, M. Weides, J. Wenner, Y. Yin, T. Yamamoto, A. N. Cleland, and J. M. Martinis, “Generation of three-qubit entangled states using superconducting phase qubits,” Nature 467, 570 (2010), arXiv:1004.4246 .
  * Dür  _et al._ (2000) W. Dür, G. Vidal, and J. I. Cirac, “Three qubits can be entangled in two inequivalent ways,” Physical Review A 62, 062314 (2000), arXiv:0005115 [quant-ph] .
  * Pezzè  _et al._ (2018) L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, “Quantum metrology with nonclassical states of atomic ensembles,” Reviews of Modern Physics 90, 035005 (2018), arXiv:1609.01609 .
  * Paulisch  _et al._ (2019) V. Paulisch, M. Perarnau-Llobet, A. González-Tudela, and J. I. Cirac, “Quantum metrology with one-dimensional superradiant photonic states,” Physical Review A 99, 043807 (2019), arXiv:1805.00712 .
  * Prevedel  _et al._ (2009) R. Prevedel, G. Cronenberg, M. S. Tame, M. Paternostro, P. Walther, M. S. Kim, and A. Zeilinger, “Experimental Realization of Dicke States of up to Six Qubits for Multiparty Quantum Networking,” Physical Review Letters 103, 020503 (2009), arXiv:0903.2212 .
  * Wieczorek  _et al._ (2009) W. Wieczorek, R. Krischek, N. Kiesel, P. Michelberger, G. Tóth, and H. Weinfurter, “Experimental Entanglement of a Six-Photon Symmetric Dicke State,” Physical Review Letters 103, 020504 (2009).
  * Chiuri  _et al._ (2012) A. Chiuri, C. Greganti, M. Paternostro, G. Vallone, and P. Mataloni, “Experimental Quantum Networking Protocols via Four-Qubit Hyperentangled Dicke States,” Physical Review Letters 109, 173604 (2012).
  * Miguel-Ramiro and Dür (2020) J. Miguel-Ramiro and W. Dür, “Delocalized information in quantum networks,” New Journal of Physics 22, 043011 (2020), arXiv:1912.12935 .
  * Cook  _et al._ (2020) J. Cook, S. Eidenbenz, and A. Bartschi, “The Quantum Alternating Operator Ansatz on Maximum k-Vertex Cover,” in  _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (IEEE, 2020) pp. 83–92, arXiv:1910.13483 .
  * Bartschi and Eidenbenz (2020) A. Bartschi and S. Eidenbenz, “Grover Mixers for QAOA: Shifting Complexity from Mixer Design to State Preparation,” in  _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (IEEE, 2020) pp. 72–82, arXiv:2006.00354 .
  * Hume  _et al._ (2009) D. B. Hume, C. W. Chou, T. Rosenband, and D. J. Wineland, “Preparation of Dicke states in an ion chain,” Physical Review A 80, 052302 (2009), arXiv:0909.0046 .
  * Stockton  _et al._ (2004) J. K. Stockton, R. van Handel, and H. Mabuchi, “Deterministic Dicke-state preparation with continuous measurement and control,” Physical Review A 70, 022106 (2004), arXiv:0402137 [quant-ph] .
  * Xiao  _et al._ (2007) Y.-F. Xiao, X.-B. Zou, and G.-C. Guo, “Generation of atomic entangled states with selective resonant interaction in cavity quantum electrodynamics,” Physical Review A 75, 012310 (2007).
  * Lamata  _et al._ (2013) L. Lamata, C. E. López, B. P. Lanyon, T. Bastin, J. C. Retamal, and E. Solano, “Deterministic generation of arbitrary symmetric states and entanglement classes,” Physical Review A 87, 032325 (2013), arXiv:1211.0404 .
  * Ivanov  _et al._ (2013) S. S. Ivanov, N. V. Vitanov, and N. V. Korolkova, “Creation of arbitrary Dicke and NOON states of trapped-ion qubits by global addressing with composite pulses,” New Journal of Physics 15, 023039 (2013).
  * Wu  _et al._ (2017) C. Wu, C. Guo, Y. Wang, G. Wang, X.-L. Feng, and J.-L. Chen, “Generation of Dicke states in the ultrastrong-coupling regime of circuit QED systems,” Physical Review A 95, 013845 (2017).
  * Chakraborty  _et al._ (2014) K. Chakraborty, B.-S. Choi, A. Maitra, and S. Maitra, “Efficient quantum algorithms to construct arbitrary Dicke states,” Quantum Information Processing 13, 2049 (2014), arXiv:1209.5932 .
  * Bärtschi and Eidenbenz (2019) A. Bärtschi and S. Eidenbenz, “Deterministic Preparation of Dicke States,” in  _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_, Vol. 11651 LNCS (2019) pp. 126–139, arXiv:1904.07358 .
  * Mukherjee  _et al._ (2020) C. S. Mukherjee, S. Maitra, V. Gaurav, and D. Roy, “On Actual Preparation of Dicke State on a Quantum Computer,”  (2020), arXiv:2007.01681 .
  * Tavis and Cummings (1968) M. Tavis and F. W. Cummings, “Exact Solution for an N-Molecule-Radiation-Field Hamiltonian,” Physical Review 170, 379 (1968).
  * Klimov and Chumakov (2009) A. B. Klimov and S. M. Chumakov,  _A Group-Theoretical Approach to Quantum Optics_ (Wiley, 2009).
  * Noguchi  _et al._ (2012) A. Noguchi, K. Toyoda, and S. Urabe, “Generation of Dicke States with Phonon-Mediated Multilevel Stimulated Raman Adiabatic Passage,” Physical Review Letters 109, 260502 (2012).
  * Unanyan  _et al._ (1999) R. G. Unanyan, B. W. Shore, and K. Bergmann, “Laser-driven population transfer in four-level atoms: Consequences of non-Abelian geometrical adiabatic phase factors,” Physical Review A 59, 2910 (1999).
  * Mousolou and Sjöqvist (2018) V. A. Mousolou and E. Sjöqvist, “Entangling power of holonomic gates in atom-based systems,” Journal of Physics A: Mathematical and Theoretical 51, 475303 (2018), arXiv:1803.00308 .
  * Kobayashi  _et al._ (2014) T. Kobayashi, R. Ikuta, S. K. Özdemir, M. Tame, T. Yamamoto, M. Koashi, and N. Imoto, “Universal gates for transforming multipartite entangled Dicke states,” New Journal of Physics 16, 023005 (2014).
  * Reagor  _et al._ (2018) M. Reagor _et al._ , “Demonstration of universal parametric entangling gates on a multi-qubit lattice,” Science Advances 4, eaao3603 (2018), arXiv:1706.06570 .
  * Caldwell  _et al._ (2018) S. A. Caldwell _et al._ , “Parametrically Activated Entangling Gates Using Transmon Qubits,” Physical Review Applied 10, 034050 (2018), arXiv:1706.06562 .
  * Nielsen (2002) M. A. Nielsen, “A simple formula for the average gate fidelity of a quantum dynamical operation,” Physics Letters A 303, 249 (2002), arXiv:0205035 [quant-ph] .
  * Willsch (2020) D. Willsch,  _Supercomputer simulations of transmon quantum computers_ , Ph.D. thesis, RWTH Aachen University (2020), arXiv:2008.13490 .
  * Koch  _et al._ (2007) J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, “Charge-insensitive qubit design derived from the Cooper pair box,” Physical Review A 76, 042319 (2007), arXiv:0703002 [cond-mat] .
  * Liu  _et al._ (2019) B.-J. Liu, X.-K. Song, Z.-Y. Xue, X. Wang, and M.-H. Yung, “Plug-and-Play Approach to Nonadiabatic Geometric Quantum Gates,” Physical Review Letters 123, 100501 (2019), arXiv:1806.07904 .
  * Abraham  _et al._ (2019) H. Abraham _et al._ , “Qiskit: An open-source framework for quantum computing,”  (2019).
  * Neeley (2010) M. G. Neeley,  _Generation of Three-Qubit Entanglement Using Josephson Phase Qubits_ , Ph.D. thesis, University of California Santa Barbara (2010).
  * Galiautdinov and Martinis (2008) A. Galiautdinov and J. M. Martinis, “Maximally entangling tripartite protocols for Josephson phase qubits,” Physical Review A 78, 010305 (2008), arXiv:0804.3159 .
  * Rigetti and Devoret (2010) C. Rigetti and M. Devoret, “Fully microwave-tunable universal gates in superconducting qubits with linear couplings and fixed transition frequencies,” Physical Review B 81, 134507 (2010).
  * Chow  _et al._ (2011) J. M. Chow, A. D. Córcoles, J. M. Gambetta, C. Rigetti, B. R. Johnson, J. A. Smolin, J. R. Rozen, G. A. Keefe, M. B. Rothwell, M. B. Ketchen, and M. Steffen, “Simple All-Microwave Entangling Gate for Fixed-Frequency Superconducting Qubits,” Physical Review Letters 107, 080502 (2011), arXiv:1106.0553 .
  * Farhi  _et al._ (2014) E. Farhi, J. Goldstone, and S. Gutmann, “A Quantum Approximate Optimization Algorithm,”  (2014), arXiv:1411.4028 .
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature Communications 5, 4213 (2014).
  * Johansson  _et al._ (2012) J. R. Johansson, P. D. Nation, and F. Nori, “QuTiP: An open-source Python framework for the dynamics of open quantum systems,” Computer Physics Communications 183, 1760 (2012), arXiv:1110.0573 .
  * Johansson  _et al._ (2013) J. R. Johansson, P. D. Nation, and F. Nori, “QuTiP 2: A Python framework for the dynamics of open quantum systems,” Computer Physics Communications 184, 1234 (2013), arXiv:1211.6518 .
  * Kay (2018) A. Kay, “Tutorial on the Quantikz package,”  (2018), arXiv:1809.03842 .
  * Niu  _et al._ (2020) Y. Niu, V. Smelyanskiy, and S. Boixo Castrillo, “Three Qubit Entangling Gate Through Two-Local Hamiltonian Control,”  (2020), US Patent Application PCT/US2019/016047.
  * Kim  _et al._ (2021) Y. Kim, A. Morvan, L. B. Nguyen, R. K. Naik, C. Jünger, L. Chen, J. M. Kreikebaum, D. I. Santiago, and I. Siddiqi, “High-fidelity iToffoli gate for fixed-frequency superconducting qubits,”  (2021), arXiv:2108.10288 .
  * Roth  _et al._ (2019) M. Roth, N. Moll, G. Salis, M. Ganzhorn, D. J. Egger, S. Filipp, and S. Schmidt, “Adiabatic quantum simulations with driven superconducting qubits,” Physical Review A 99, 022323 (2019), arXiv:1808.04666 .
  * Buckle  _et al._ (1986) S. J. Buckle, S. M. Barnett, P. L. Knight, M. A. Lauder, and D. T. Pegg, “Atomic Interferometers,” Optica Acta: International Journal of Optics 33, 1129 (1986).