# Joint Cosmic Density Reconstruction from Photometric and Spectroscopic
Samples

B. Horowitz,1 P. Melchior2,3
1Lawrence Berkeley National Lab, 1 Cyclotron Road, Berkeley, CA 94720, USA
2Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA
3Center for Statistics & Machine Learning, Princeton University, Princeton, NJ
08544, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We reconstruct the dark matter density field from spatially overlapping
spectroscopic and photometric redshift catalogs through a forward modelling
approach. Instead of directly inferring the underlying density field, we find
the best fitting initial Gaussian fluctuations that will evolve into the
observed cosmic volume. To account for the substantial uncertainty of
photometric redshifts we employ a differentiable continuous Poisson process.
In the context of the upcoming Prime Focus Spectrograph (PFS), we find
improvements in cosmic structure classification equivalent to 50-100% more
spectroscopic targets by combining relatively sparse spectroscopic with dense
photometric samples.

###### keywords:

methods: data analysis – large-scale structure of Universe – galaxies:
distances and redshifts

††pubyear: 2023††pagerange: Joint Cosmic Density Reconstruction from
Photometric and Spectroscopic Samples–Joint Cosmic Density Reconstruction from
Photometric and Spectroscopic Samples

##  1 Introduction

For the past half century, spectroscopic galaxy redshift surveys have been the
key method used to study the three-dimensional distribution of matter in our
universe (Gregory & Thompson, 1978; Kirshner et al., 1978; Davis et al.,
1982). This intricate structure, known as the cosmic web, consists of under-
dense void regions enclosed by sheet-like walls, which are embedded between
filamentary structures that connect dense dark matter halos of galaxy clusters
(Bond et al., 1996b). Signatures of these structures are well established
(e.g. Platen et al. (2011); Tempel et al. (2014)), and there is significant
interest in connecting galaxy properties, like inferred star formation rate,
with their cosmic environments (Poudel et al., 2017; Greene et al., 2022b).

Photometric surveys, on the other hand, provide significant cosmological
information through their clustering statistics (Nicola et al., 2020),
gravitational lensing effect (Kilbinger et al., 2013), and cross correlations
with spectroscopic catalogs (Coupon et al., 2015). Upcoming surveys, such as
the Legacy Survey of Space and Time (LSST, LSST Science Collaboration et al.,
2009), and ongoing surveys like the Hyper Suprime-Camera (HSC) on the Subaru
Telescope are greatly expanding the number of galaxies with multiwavelength
data suitable for photometric analysis (Wen & Han, 2021). However, the
substantial uncertainty in their inferred redshifts (“photo-z”) is a major
limitation for their utility in tracing cosmic structure.

The past few years has seen an explosion of work in foreword modelling of
galaxy density fields and applying those techniques to real survey data. These
include the Bayesian Origin Reconstruction from Galaxies method (Jasche &
Wandelt, 2013, BORG) applied to the SDSS/BOSS galaxy survey (Lavaux et al.,
2019), and the more recent COSMological Initial Conditions from Bayesian
Inference Reconstructions with THeoretical models (Kitaura et al., 2020,
COSMIC BIRTH) applied to the COSMOS field (Ata et al., 2021). Both these
techniques rely on computationally expensive Hamiltonian Monte Carlo (HMC)
sampling techniques, requiring on the order of 10,000 samples, each requiring
a number of forward model computations. An alternative method was proposed in
Seljak et al. (2017), where the authors pointed out that a Gaussian
approximation for the likelihood is fairly accurate for cosmological analysis,
and a simple quasi-Newtonian optimization scheme can achieve accurate results
in $\mathcal{O}(100)$ steps. An implementation of this method was first
applied to a realistic galaxy mock catalog in Horowitz et al. (2021b), finding
an accurate reconstruction of the cosmic web. These techniques have since been
combined with Lyman-$\alpha$ forest measurements and been used for map making
from data from the COSMOS Lyman Alpha Mapping And Tomographic
Observations(Horowitz et al., 2021a).

A critical aspect for both the quasi-Newtonian optimization and the HMC
analysis is the model to map from the evolved dark matter distribution to the
observed galaxies which act as biased tracers. These methods must be fully
differentiable; i.e. one must be able to calculate the derivative of the final
(gridded) tracer density field with respect to the initial matter density
field and/or particle positions. While their have been a number of fast
particle-mesh solvers developed over the past decade (Merz et al., 2005;
Tassev et al., 2013; Izard et al., 2016), only a few have been designed with
integrated derivative calculations (Feng et al., 2016; Modi et al., 2020). In
addition, one must still further map from the evolved density field (or the
found displacement field) to the tracer density in a differentiable way.

In this work, we present an approach to jointly modelling both spectroscopic
and photometric redshift catalogs with the goal of cosmic web reconstruction.
We extend the forward galaxy modelling framework in Horowitz et al. (2021b) to
account for significant uncertainty in galaxy positions through a Poissonian
likelihood. We present our likelihood model and mock catalog construction in
Section 2. We show our reconstructions of the cosmic web and compare with our
simulated truth in Section 3 and finally discuss the results and further
applications in Section 4. Throughout this work we assume a Planck 2015
cosmology as discussed in Planck Collaboration et al. (2016).

##  2 Methods

###  2.1 Dynamical Forward Model

We follow past works (Horowitz et al., 2019, 2021b) and use a dynamical model
to forward simulate our observed universe and iteratively find the maximum a
posteriori initial conditions. In short, we solve for initial density field,
$\hbox{\boldmath{$s$}}\in\mathbb{R}^{M}$, given some observed data
$\hbox{\boldmath{$d$}}\in\mathbb{R}^{N}$. We can map our initial density field
to the data with a forward operator, $f$. In general this forward operator
will depend on additional parameters (i.e. cosmological or astrophysical
parameters) which will need to be marginalized over depending on the final
quantity of interest. For this work we assume these are fixed (see Seljak et
al. (2017) for a discussion of this marginalization).

This is a highly dimensional optimization problem, approximately equal to the
number of initial particles in our reconstruction simulation. To use gradient
based optimization methods, like L-BFGS, we need to have a fully
differentiable forward model to map from our initial density field to the late
time observed quantities. The dynamical evolution of the field itself is done
with FlowPM (Modi et al., 2020), a tensorflow implementation of the FastPM
(Feng et al., 2016) package. This implementations has the advantage of GPU-
compatibility while also being completely differentiable. The FastPM framework
has been well tested to reproduce halo statistics, power spectrum, and
correlation functions as more traditional n-body solvers (TreePM). We use 8
timesteps to evolve the field from $z=10$ to $z=1.0$, linearly spaced in scale
factor, with a force resolution of 1. We choose $z=1.0$ as it is the mean
target redshift for the PFS-Galaxy Evolution low-$z$ survey, the main
motivation for this work.

For our differentiable galaxy forward model we use second order Lagrangian
bias using the displacement, $\psi$, calculated from our particle-mesh
simulation. In this framework, the galaxy field $\delta_{g}$ can be expressed
as a

$1+\delta_{g}(\hbox{\boldmath{$x$}})=\int
d^{3}\hbox{\boldmath{$q$}}w(\hbox{\boldmath{$q$}})\delta_{D}(\hbox{\boldmath{$x$}}-\hbox{\boldmath{$q$}}-\psi),$
(1)

where $w(\hbox{\boldmath{$q$}})$ is a weight function (i.e. response)
expressed as a sum of of the Lagrangian density as

$w(\hbox{\boldmath{$q$}})=1+b_{1}\delta(\hbox{\boldmath{$q$}})+b_{2}\left[\delta^{2}(\hbox{\boldmath{$q$}})-\langle\delta^{2}(\hbox{\boldmath{$q$}})\rangle\right].$
(2)

Note that, in practice, this integral expression can be done trivially at the
field level via painting our weighted field on the initial (i.e. Lagrangian)
density field and applying displacement calculated by the dynamical forward
model. Additional bias terms (such as tidal shear) could be added if they can
be expressed as differentiable operators of outputs from our particle mesh
simulation. We calibrate $b_{1}$ and $b_{2}$ by fitting the mock observed
power spectra to that from a separate mock catalog (see Horowitz et al.
(2021b) for additional discussion of this approach). Bias parameters could
also be fitted jointly with the optimization, with an increased computational
cost and possible decrease in reconstruction accuracy.

###  2.2 Galaxy and Photometric Redshift Likelihood Model

At each step in our optimization, we need to compare our reconstructed dark
matter field with the (true) observed photometric galaxy catalog. In Horowitz
et al. (2021b), the authors performed an density reconstruction of
spectroscopic observed galaxies. In that work the authors used a $\ell_{2}$
norm for comparing the observed galaxies with a forward modelled bias density
field. As studied in Nguyen et al. (2021), a Poisson likelihood term provides
a more accurate reconstruction of small-scale power and, in the case of coarse
galaxy resolution (i.e. $\lesssim 5$ galaxies per grid square) the Gaussian
approximation is arguably no longer mathematically valid.

With this in mind, we model our observed catalog as a draw from the
inhomogeneous Poisson process in $\hbox{\boldmath{$x$}}\in\mathbb{R}^{3}$,
defined by the observable rate density
$\Phi(\hbox{\boldmath{$x$}}\mid\hbox{\boldmath{$s$}})=f(\hbox{\boldmath{$x$}}\mid\hbox{\boldmath{$s$}})\mathcal{W}(\hbox{\boldmath{$x$}})$,
a product of the sample density predicted by the forward model and the survey
detection efficiency $\mathcal{W}$.

We can write the probability of the observed data for $K$ objects,
$\hbox{\boldmath{$d$}}=\\{\hbox{\boldmath{$x$}}_{1},\cdots,\hbox{\boldmath{$x$}}_{K}\\}$,
given some initial density field as

$\mathcal{P}(\hbox{\boldmath{$d$}}\mid\hbox{\boldmath{$s$}})=\exp\left(-\int
d^{3}x\,\Phi(\hbox{\boldmath{$x$}}\mid\hbox{\boldmath{$s$}})\right)\prod_{k=1}^{K}\Phi(\mathbf{x}_{k}\mid\hbox{\boldmath{$s$}}).$
(3)

In this formula, the exponential term acts as a normalization to express the
expected total number of observations. As our forward model is binned in
cells, we transform the product over galaxies to a product over grid cells:

$\displaystyle\log\mathcal{P}(\hbox{\boldmath{$d$}}\mid\hbox{\boldmath{$s$}})$
$\displaystyle=\sum_{i}n_{i}\log\phi_{i}-\phi_{i}\Delta_{i}^{3}\mathrm{,with}$
(4) $\displaystyle n_{i}$
$\displaystyle=\sum_{k}\iota(x_{k}\in\mathcal{C}_{i})\ \mathrm{and}$ (5)
$\displaystyle\phi_{i}$
$\displaystyle=\int_{\mathcal{C}_{i}}d^{3}x\,\Phi(\hbox{\boldmath{$x$}}\mid\hbox{\boldmath{$s$}})$
(6)

being the number of observed samples that fall in cell $\mathcal{C}_{i}$ and
the total expected rate in that cell, respectively. $\Delta^{3}_{i}$ is the
cell volume.

Equation 5 requires that we can directly count the number of times a galaxy
falls into a specific cell, which implies positional uncertainties much
smaller than the size of the cells. This requirement is concerning already for
spectroscopic samples, but it cannot be met with photometric ones with dozens
of Mpc uncertainties in the redshift direction. We can generalize the number
counts to account for the redshift uncertainty $p(z_{k})$ of each galaxy in
that cell:

$n^{\prime}_{i}=\sum_{k}\mathcal{P}\left(x_{k}\in\mathcal{C}_{i}\mid
p(z_{k})\right).$ (7)

This amounts to a redistribution of the counts according to the probability
mass every $p(z_{k})$ assigns to each cell. If all galaxies are assumed to
have the same redshift uncertainty, Equation 7 simplifies to a one-dimensional
convolution of the observed counts with the redshift kernel in the radial
direction, but our approach does not require this assumption. Because the
modified counts are no longer integer, we employ the generalization of the
Poisson distribution for continuous variables, which can be expressed in terms
of the incomplete Gamma function probability:

$\mathcal{P}(n^{\prime}_{i}\mid\hbox{\boldmath{$s$}})=\frac{1}{\Gamma(n^{\prime}_{i})}\int_{0}^{\phi(\hbox{\boldmath{$s$}})_{i}}t^{n^{\prime}_{i}-1}e^{-t}dt,$
(8)

where $\Gamma$ is the normal gamma function.

Figure 1: Projection along the $z$ axis for a 10 $h^{-1}$ Mpc slice of a simulated volume as reconstructed by a mock spectroscopic (PFS-Cosmo), photometric (HSC), and jointly both surveys. We find that the photometric reconstruction provides an accurate but diffuse map of the underlying structures, while the spectroscopic reconstruction is able to pick out the most massive clusters but lacks the number density to follow the low-density structures. The joint reconstruction is able to capture both features and shows good qualitative agreement to the true density field. | HSC | PFS-Cosmo | PFS-GE
---|---|---|---
Number Density (h-1 Mpc)-3 | 0.06 | 0.0006 | 0.025
Redshift Error (h-1Mpc) | 30.0 | 0.5 | 0.5
Table 1: Configurations chosen to mimic the upcoming Prime Focus Spectrograph
(PFS) survey and the Hyper Suprime-Cam (HSC) photometric survey.

###  2.3 Mock Catalog and Redshift Error Modeling

We construct a mock simulated volume using FastPM, which has been well tested
to reconstruct proper halo statistics and profiles for the range of interest
in this work (Feng et al., 2016). We use a simulated volume of $128^{3}$
$\,h^{-3}\,\mathrm{Mpc}^{3}$with particle resolution of $1024^{3}$, with time
resolution of 10 steps. We use a friends of friends halo finder from Nbodykit
(Hand et al., 2018), finding 2,097,149 halos in our volume. We use the Zheng
et al. (Zheng et al., 2005) HOD model from halotools (Hearin et al., 2016)
with associated HOD parameters $\log{M_{\textrm{min}}}=12.02$,
$\sigma_{\log{M}}=0.26$, $M_{0}=11.38$, $\log{M_{1}}=13.31$ and $\alpha=1.06$,
matching those found in Zheng et al. (2007). This results in 91588 central and
satellite galaxies, from which we down-sample to form our mock catalogs for
each survey.

For our mock catalog we model the Prime Focus Spectrograph’s cosmology survey
or PFS Galaxy Evolution survey combined with the overlapping HSC targetting
photometric catalog. The PFS Cosmology survey is designed to reconstruct the
growth rate of structure and geometry via a variety of probes on large scales
(Bayronic acoustic oscillations, redshift space distortions, clustering)
(Takada et al., 2014) while the Galaxy Evolution survey is currently designed
to reach a significantly larger number density in order to resolve cosmic web
structures and identify galaxy groups and clusters (Greene et al., 2022a). At
$z\sim 1$, each tasks require different number densities of spectroscopic
targets, and we show the corresponding number densities for each sample in
Table 1.

We model our photometric redshift errors inspired by the average photo-z
calibration error described in HSC PDR2 (Nishizawa et al., 2020). We assume
that all galaxies have a similar photometric error distribution, with Gaussian
error of 30 Mpc/h corresponding to a $\Delta z\sim 0.015$ at $z=1.0$. More
complex photo-z distributions can be included in a straightforward fashion via
re-expression as a Gaussian mixture model and adjusting the rate density in
Eq. LABEL:eq:prob2. When combining HSC with a spectroscopic data set, we
assume all spectroscopically identified galaxies are also in the HSC catalog
and remove them in order to not double count a given galaxy.

We assume all galaxies have the same bias properties, but it would be simple
within our formalism to use multiple different types of biased galaxy subsets,
each with their own log likelihood that is then added together. We smooth the
observed three-dimensional galaxy field by 0.05 $h^{-1}$ Mpc because sharp
discontinuities can lead to a “stuck" optimization at a non-optimal solution.

###  2.4 Optimization Method and Prior

Our method relies on an high-dimensional optimization for the underlying
initial density field, $\hbox{\boldmath{$s$}}\in\mathbb{R}^{M}$. The
probability of some given observed data $d$ can be expressed as

$\mathcal{P}(\hbox{\boldmath{$d$}})=\int
d\hbox{\boldmath{$s$}}\mathcal{P}(\hbox{\boldmath{$s$}})\mathcal{P}(\hbox{\boldmath{$d$}}\mid\hbox{\boldmath{$s$}}),$
(9)

where $\mathcal{P}(\hbox{\boldmath{$s$}})$ is a prior on the initial density
field and $\mathcal{P}(\hbox{\boldmath{$d$}}\mid\hbox{\boldmath{$s$}})$ is the
likelihood with respect to the fields. This expression can further be
constrained by imposing a constraint on the signal field based on cosmological
model. Assuming the initial density is a Gaussian random field described by
covariance matrix, $S$, we can write down the prior as

$\mathcal{P}(\hbox{\boldmath{$s$}}\mid\hbox{\boldmath{$S$}})=\frac{\exp{\left(-\frac{1}{2}\hbox{\boldmath{$s$}}^{\dagger}\hbox{\boldmath{$S$}}^{-1}\hbox{\boldmath{$s$}}\right)}}{\det{2\pi\hbox{\boldmath{$S$}}}},$
(10)

Usually the dimension of the initial field, $M$ will be quite large, resulting
in a very high-dimensional inverse problem. For example, in Horowitz et al.
(2019) this procedure was used to reconstruct the initial density field from
Lyman Alpha forest measurements and required $M=128^{3}$ for an accurate
reconstruction over a $64$ $h^{-1}$ Mpc box. Optimization over such a large
number of parameters would be impractical without using derivative information
to inform each update in the optimization. We use an Limited Memory BFGS
optimization method, implemented in scipy, which uses derivative information
at each step to update a Hessian approximation at each step (a “quasi-
Newtonian" method).

In general there is no guarantee the likelihood surface is well behaved and
convex. However, for the examples shown in this work we find little
quantitative difference in our optimization by starting our optimization at
different initial density fields. Alternative methods, like Hamiltonian Monte
Carlo, could be used if this is a major concern (Kitaura et al., 2020; Jasche
& Wandelt, 2013). This is particularly true if one wants to use the
reconstructed density field to update the priors of a redshift of a given
galaxy (similar to clustering redshift calculations).

##  3 Results

We have applied the reconstruction algorithm to our mock catalogs representing
HSC and PFS-Cosmo and show a projected slice of the reconstructions at $z=1.0$
in Figure 1. Qualitatively, we find strong agreement between the fields, with
structures in the reconstructions appearing more diffuse than the
corresponding true structures. The spectroscopic reconstructions suffers from
sparse sampling and can thus not recover smaller structures. In contrast, the
higher-density photometric catalog allows to probe smaller and lower-density
structures but lacks the spatial sharpness to accurately recover compact
massive structures resulting in washed out features.

For a quantitative comparison between the reconstructions at the field level
we can compare the fields pixel by pixel via their Pearson correlation
coefficient, defined as

$r_{xy}=\frac{\sum^{n}_{i=1}(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\sum^{n}_{i=1}(x_{i}-\bar{x})^{2}}\sqrt{\sum^{n}_{i=1}(y_{i}-\bar{y})^{2}}}.$
(11)

We plot this coefficient as a function of smoothing scale in Figure 2 for both
spectroscopic and photometric samples, as well as the joint analysis. We also
show the result for using the higher density PFS-Galaxy Evolution survey as
well as using half the density of that survey. We find, at a scale of 1
$\,h^{-1}\,\mathrm{Mpc}\,$, the joint reconstruction performs comparably to
the half of the PFS-Galaxy Evolution survey alone despite requiring
significantly fewer spectroscopic targets (1260 vs. 26140).

Figure 2: Pearson correlation coefficient as a function of smoothing scales
between the reconstructed volume and the simulated truth. As we broaden the
Gaussian smoothing kernel, we are able to match the simulated truth better, at
the cost of small scale detail. We include reconstructions using half of the
PFS-GE sample density and the full PFS-GE sample. See Table 1.  Figure 3:
Qualitative cosmic structure plot for a sample x-z slice in the reconstructed
volume classified using the deformation tensor approach for the spectroscopic
redshift (PFS-Cosmo), photometric redshift (HSC) and joint surveys compared to
the simulated truth. Dark blue corresponds to nodes (i.e. clusters), light
blue to filaments, green to sheets, and yellow to voids. Figure 4: Cosmic web
classification confusion matrix for our reconstruction for the spectroscopic
redshift (PFS-Cosmo), photometric redshift (HSC) and joint surveys Perfect
reconstruction would correspond to no off-diagonal terms. The joint
reconstruction uniformly improves all cosmic structure classifications. Note
that there is significant confusion between notes and filaments (and a lesser
extent between sheets and voids) due to difficulty in exactly locating
structures along the line of sight. Figure 5: Histogram showing the
eigenvector orientation recovery, corresponding to the principle inflow
(outflow) directions of the cosmic structures, by computing the dot product of
the reconstructed eigenvector with the true eigenvector. Perfect alignment
corresponds to all weight in the far right side bin, and random alignment
corresponds to the black dotted line. We also show (in red) the structured
inferred by taking the spectroscopic galaxy field and smoothing it at the same
resolution as the other fields, 2 $h^{-1}$Mpc.

###  3.1 Structure Classification: T-Web

Following Horowitz et al. (2019, 2021b), we use an the eigenvalues and vectors
of the pseudo-deformation tensor, sometimes called the “T-web" classification,
as described in Lee & White (2016); Krolewski et al. (2017) and based on work
in Bond et al. (1996a); Hahn et al. (2007); Forero-Romero et al. (2009). The
pseudo-deformation tensor is the Hessian of the gravitation potential and its
eigenvectors describe the principle inflow (outflow) directions assuming the
corresponding eigenvalue is positive (negative) assuming the Zel’dovich
approximation (Zel’dovich, 1970). This can be efficiently computed in Fourier
space with matter density, $\delta_{k}$, as

$\tilde{D}_{ij}=\frac{k_{i}k_{j}}{k^{2}}\delta_{k}.$ (12)

We use this definition, as opposed to using the reconstructed nonlinear
velocity field, in order to make a direct comparison to existing works in the
literature which use this definition.

To classify cosmic structure into discrete classes we compare the eigenvalues
to some thresh-hold value, $\lambda_{th}$. While one could use only the
signature of the eigenvalues (i.e. $\lambda_{th}=0$), this would result in a
very low node fraction as only the pixel at local maxima would be considered
true nodes. Following the approach in Horowitz et al. (2019, 2021b), we
instead use set our thresh-hold for each map such that all our reconstructed
fields have the same void volume fraction of 25%. We use voids, rather than
nodes, as they are less sensitive to the effective smoothing scale of the
tracers. We show the classified fields in Figure 3. We can compare pixel by
pixel to compute a confusion matrix to see how often given cosmic structures
are mis-classified, which is shown in Figure 4. We find the joint
reconstruction performs significantly better than the photometric catalog
alone and moderately better than the spectroscopic catalog. In particular we
find substantially fewer significant misclassifications (i.e. node
misidentified as void or visa versa).

In addition to discrete classification based on eigenvalues, the eigenvectors
of the deformation tensor can be used to categorize the orientations of the
cosmic structures. Each eigenvector indicates the principle inflow or outflow
direction from the structure based on their signature. We can then determine
the angle $\theta$ between these eigenvectors and the simulated true
eigenvector, with the ideal result being a distribution that peaks at
$\cos(\theta)=1$. We show this alignment in Figure 5. We find the HSC
photometric catalogs contribute moderately to the reconstructed alignment in
the joint catalog with PFS-Cosmo but are poor tracers alone of the alignments.

###  3.2 Halo Recovery

Figure 6: Reconstructed halo mass from the joint PFS-Cosmo + HSC catalog
defined by a spherical overdensity vs. the true mass in the corresponding
overdensity in our simulated density field. We find reasonable reconstruction
of halo mass, with $\frac{\delta M}{M}=0.42$ for halos with $M>2\times
10^{13}M_{\odot}$.

Instead of the continuous eigenvalue description used in Sec 3.1, we can
perform a discrete classification based on a halo-finding algorithm. This
allows us to identify and easily compare individual cosmic structures, for
example to identify galaxies to the same massive halo.

There are a number of possible ways of determining halos, including using a
friends-of-friends algorithm on the particle data (as done in Horowitz et al.
(2021b)) or using spherical over-density (for example in Stark et al. (2015)).
For this work, we use halos identified by a friends-of-friends algorithm (from
nbodykit (Hand et al., 2018)) from our reconstructed particle catalog. As we
are interested in accuracy of the overall environment of galaxy clusters as
opposed to construction of a specific halo catalog for cosmological analysis,
we perform our comparison not on the friends-of-friends mass but on a
spherical volume surrounding the center of each identified halo. This allows
direct comparison of the halos identified in the reconstruction volume to
those in the simulated truth without needing to account for specific halo
identification between the two catalogs. This also minimizes the effects of
minor variations in halo detection parameters causing a true halo to
“fragment” into two smaller halos. We compare the recovered mass of our
catalog to the true mass in Figure 6. We find $\frac{\delta M}{M}=0.42$ for
halos with $M>2\times 10^{13}M_{\odot}$.

##  4 Conclusion and Discussion

In this work, we have presented the first joint analysis of a spectroscopic
and photometric galaxy samples for reconstruction of the initial density field
of an observed cosmic volume. We have shown that even a fairly sparse
spectroscopic galaxy sample, when combined with a deep photometric galaxy
catalog, can provide an accurate reconstruction of the cosmic web and recover
major over-densities. We have demonstrated the competitive performance with
the cosmic web classification by a deformation tensor approach as well as via
the statistics of spherical over/under-densities. We have shown this method
with emphasis on the upcoming Prime Focus Spectroscopic (PFS) Survey, combined
with the Subaru Hyper Suprime-Cam surveys photometry. While cosmic web
reconstruction is a the focus for the PFS - Galaxy Evolution survey, the
method presented here is entirely general and could also be applied to
combined the upcoming Legacy Survey of Space and Time (LSST) (LSST Science
Collaboration et al., 2009) survey combined with a variety of spectroscopic
surveys.

In Horowitz et al. (2021b), a similar approach was used to perform a joint
reconstruction using Ly$\alpha$ forest data and spectroscopic redshift
catalogs. In that case the data was highly complementary and the
reconstruction synergistic as each probe provided information on different
spatial scales. In this work, we find the effects are more additive since the
probes have comparable bias properties and probe comparable physical scales.
In Modi et al. (2021), the authors combined photometric surveys and intensity
mapping in a comparable forward modelling framework. However, the focus of
that work was on reconstructing large scales for cosmological analysis and
didn’t consider late-time cosmic web statistics, and used significantly
different likelihoods due to that application.

Recently, several works employed modern machine learning techniques to perform
a more accurate photometric redshift prediction (e.g. Eriksen et al. (2020);
Mu et al. (2020); Abul Hayat et al. (2021))111For an updated more complete
list, see https://github.com/georgestein/ml-in-cosmology.. Approaches that use
individual galaxies photometry to infer redshifts can be trivially
incorporated into these reconstruction techniques. More nuanced approaches
that rely on correlations between galaxies, such as clustering redshift
approaches, require more care as they likely will include correlations in the
redshift errors which would need to be propagated through the likelihood
function.

In this work we use a biasing scheme to map from density to galaxy field.
While these approaches have been demonstrated to work well at describing
structure on large scales, on small scales they are known to break down
(Desjacques et al., 2018). An alternative approach could be to include in our
forward model an explicit halo finding and halo population step.
Differentiable versions of halo finding has been explored using neural network
approaches in (Modi et al., 2018, 2020; Kodi Ramanah et al., 2019) and a
differential HOD model has been developed in Horowitz et al. (2022). We plan
to further explore applying these methods in upcoming work.

While not directly addressed in this work, the use of photometric redshift
catalogs could help mitigate the effects of fiber collisions in cosmic web
reconstructions (Dawson et al., 2013; Makarem et al., 2014). In spectroscopic
surveys covering dense cluster environments there is difficulty in achieving
uniform sampling due to mechanical limitations of the fiber positioners. Since
these environments are still well sampled by photometric surveys (assuming
proper deblending), our reconstruction technique could substantially reduce
the bias from this localized incompleteness.

## Acknowledgements

This work is supported by the AI Accelerator program of the Schmidt Futures
Foundation.

This research used resources of the National Energy Research Scientific
Computing Center, a DOE Office of Science User Facility supported by the
Office of Science of the U.S. Department of Energy under Contract No.
DEC02-05CH11231.

This work made extensive use of scipy (Virtanen et al., 2020) and Tensorflow
(Abadi et al., 2015).

## References

  * Abadi et al. (2015) Abadi M., et al., 2015, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, https://www.tensorflow.org/
  * Abul Hayat et al. (2021) Abul Hayat M., Harrington P., Stein G., Lukić Z., Mustafa M., 2021, arXiv e-prints, p. arXiv:2101.04293
  * Ata et al. (2021) Ata M., Kitaura F.-S., Lee K.-G., Lemaux B. C., Kashino D., Cucciati O., Hernández-Sánchez M., Le Fèvre O., 2021, MNRAS, 500, 3194
  * Bond et al. (1996a) Bond J. R., Kofman L., Pogosyan D., 1996a, Nature, 380, 603
  * Bond et al. (1996b) Bond J. R., Kofman L., Pogosyan D., 1996b, Nature, 380, 603
  * Coupon et al. (2015) Coupon J., et al., 2015, MNRAS, 449, 1352
  * Davis et al. (1982) Davis M., Huchra J., Latham D. W., Tonry J., 1982, The Astrophysical Journal, 253, 423
  * Dawson et al. (2013) Dawson K. S., et al., 2013, AJ, 145, 10
  * Desjacques et al. (2018) Desjacques V., Jeong D., Schmidt F., 2018, Physics reports, 733, 1
  * Eriksen et al. (2020) Eriksen M., et al., 2020, MNRAS, 497, 4565
  * Feng et al. (2016) Feng Y., Chu M.-Y., Seljak U., McDonald P., 2016, MNRAS, 463, 2273
  * Forero-Romero et al. (2009) Forero-Romero J. E., Hoffman Y., Gottlöber S., Klypin A., Yepes G., 2009, MNRAS, 396, 1815
  * Greene et al. (2022a) Greene J., Bezanson R., Ouchi M., Silverman J., et al., 2022a, arXiv preprint arXiv:2206.14908
  * Greene et al. (2022b) Greene J., Bezanson R., Ouchi M., Silverman J., the PFS Galaxy Evolution Working Group 2022b, arXiv e-prints, p. arXiv:2206.14908
  * Gregory & Thompson (1978) Gregory S. A., Thompson L. A., 1978, ApJ, 222, 784
  * Hahn et al. (2007) Hahn O., Carollo C. M., Porciani C., Dekel A., 2007, MNRAS, 381, 41
  * Hand et al. (2018) Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, The Astronomical Journal, 156, 160
  * Hearin et al. (2016) Hearin A. P., Zentner A. R., van den Bosch F. C., Campbell D., Tollerud E., 2016, MNRAS, 460, 2552
  * Horowitz et al. (2019) Horowitz B., Lee K.-G., White M., Krolewski A., Ata M., 2019, ApJ, 887, 61
  * Horowitz et al. (2021a) Horowitz B., et al., 2021a, arXiv e-prints, p. arXiv:2109.09660
  * Horowitz et al. (2021b) Horowitz B., Zhang B., Lee K.-G., Kooistra R., 2021b, ApJ, 906, 110
  * Horowitz et al. (2022) Horowitz B., Hahn C., Lanusse F., Modi C., Ferraro S., 2022, arXiv e-prints, p. arXiv:2211.03852
  * Izard et al. (2016) Izard A., Crocce M., Fosalba P., 2016, MNRAS, 459, 2327
  * Jasche & Wandelt (2013) Jasche J., Wandelt B. D., 2013, MNRAS, 432, 894
  * Kilbinger et al. (2013) Kilbinger M., et al., 2013, MNRAS, 430, 2200
  * Kirshner et al. (1978) Kirshner R., Oemler Jr A., Schechter P., 1978, The Astronomical Journal, 83, 1549
  * Kitaura et al. (2020) Kitaura F.-S., Ata M., Rodríguez-Torres S. A., Hernández-Sánchez M., Balaguera-Antolínez A., Yepes G., 2020, MNRAS,
  * Kodi Ramanah et al. (2019) Kodi Ramanah D., Charnock T., Lavaux G., 2019, Phys. Rev. D, 100, 043515
  * Krolewski et al. (2017) Krolewski A., Lee K.-G., Lukić Z., White M., 2017, ApJ, 837, 31
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201
  * Lavaux et al. (2019) Lavaux G., Jasche J., Leclercq F., 2019, arXiv e-prints, p. arXiv:1909.06396
  * Lee & White (2016) Lee K.-G., White M., 2016, ApJ, 831, 181
  * Makarem et al. (2014) Makarem L., Kneib J.-P., Gillet D., Bleuler H., Bouri M., Jenni L., Prada F., Sanchez J., 2014, A&A, 566, A84
  * Merz et al. (2005) Merz H., Pen U.-L., Trac H., 2005, New Astron., 10, 393
  * Modi et al. (2018) Modi C., Feng Y., Seljak U., 2018, J. Cosmology Astropart. Phys., 2018, 028
  * Modi et al. (2020) Modi C., Lanusse F., Seljak U., 2020, arXiv e-prints, p. arXiv:2010.11847
  * Modi et al. (2021) Modi C., White M., Castorina E., Slosar A., 2021, J. Cosmology Astropart. Phys., 2021, 056
  * Mu et al. (2020) Mu Y.-H., Qiu B., Zhang J.-N., Ma J.-C., Fan X.-D., 2020, Research in Astronomy and Astrophysics, 20, 089
  * Nguyen et al. (2021) Nguyen N.-M., Schmidt F., Lavaux G., Jasche J., 2021, J. Cosmology Astropart. Phys., 2021, 058
  * Nicola et al. (2020) Nicola A., et al., 2020, J. Cosmology Astropart. Phys., 2020, 044
  * Nishizawa et al. (2020) Nishizawa A. J., Hsieh B.-C., Tanaka M., Takata T., 2020, arXiv e-prints, p. arXiv:2003.01511
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Platen et al. (2011) Platen E., van de Weygaert R., Jones B. J. T., Vegter G., Calvo M. A. A., 2011, MNRAS, 416, 2494
  * Poudel et al. (2017) Poudel A., Heinämäki P., Tempel E., Einasto M., Lietzen H., Nurmi P., 2017, A&A, 597, A86
  * Seljak et al. (2017) Seljak U., Aslanyan G., Feng Y., Modi C., 2017, preprint,  (arXiv:1706.06645)
  * Stark et al. (2015) Stark C. W., White M., Lee K.-G., Hennawi J. F., 2015, MNRAS, 453, 311
  * Takada et al. (2014) Takada M., et al., 2014, PASJ, 66, 1
  * Tassev et al. (2013) Tassev S., Zaldarriaga M., Eisenstein D. J., 2013, J. Cosmology Astropart. Phys., 2013, 036
  * Tempel et al. (2014) Tempel E., Stoica R. S., Martínez V. J., Liivamägi L. J., Castellan G., Saar E., 2014, MNRAS, 438, 3465
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wen & Han (2021) Wen Z. L., Han J. L., 2021, MNRAS, 500, 1003
  * Zel’dovich (1970) Zel’dovich Y. B., 1970, A&A, 5, 84
  * Zheng et al. (2005) Zheng Z., et al., 2005, ApJ, 633, 791
  * Zheng et al. (2007) Zheng Z., Coil A. L., Zehavi I., 2007, ApJ, 667, 760