# Cosmic variance suppression in radiation-hydrodynamic modeling of the
reionization-era 21-cm signal

Anshuman Acharya1, Enrico Garaldi1, Benedetta Ciardi1, Qing-bo Ma2,3
1Max-Planck-Institut für Astrophysik, Garching 85748, Germany
2School of Physics and Electronic Science, Guizhou Normal University, Guiyang
550001, PR China
3Guizhou Provincial Key Laboratory of Radio Astronomy and Data Processing,
Guizhou Normal University, Guiyang 550001, PR China E-mail: anshuman@mpa-
garching.mpg.de

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The 21-cm line emitted by neutral hydrogen is the most promising probe of the
Epoch of Reionization (EoR). Multiple radio interferometric instruments are on
the cusp of detecting its power spectrum. It is therefore essential to deliver
robust theoretical predictions, enabling sound inference of the coeval
Universe properties. The nature of this signal traditionally required the
modeling of $\mathcal{O}(10^{7-8}\,{\rm Mpc}^{3})$ volumes in order to
suppress the impact of cosmic variance. However, the recently-proposed Fixed &
Paired (F&P) approach (Pontzen et al., 2016) uses carefully-crafted simulation
pairs to achieve equal results in smaller volumes. In this work, we thoroughly
test the applicability of and improvement granted by this technique to
different observables of the 21-cm signal from the EoR. We employ radiation-
magneto-hydrodynamics simulations to ensure the most realistic physical
description of this epoch, greatly improving over previous studies using a
semi-numerical approach without accurate galaxy formation physics and
radiative transfer. We estimate the statistical improvement granted by the F&P
technique on predictions of the skewness, the power spectrum, the bispectrum
and ionized regions size distribution of the 21-cm signal at redshift $7\leq
z\leq 10$. We find that the effective volume of F&P simulations is at least
3.5 times larger than in traditional simulations. This directly translates
into an equal improvement in the computational cost (both in terms of time and
memory). Finally, we confirm that a combination of different observables like
skewness, power spectrum and bispectrum across different redshifts can be
utilised to maximise the improvement brought by the F&P technique.

###### keywords:

cosmology: dark ages, reionization, first stars; cosmology: theory; galaxies:
formation; methods: data analysis

††pubyear: 2023††pagerange: Cosmic variance suppression in radiation-
hydrodynamic modeling of the reionization-era 21-cm signal–Cosmic variance
suppression in radiation-hydrodynamic modeling of the reionization-era 21-cm
signal

##  1 Introduction

One of the most important probes for studying the Epoch of Reionization (EoR;
the period of the Universe when the first astrophysical objects formed and
emitted UV photons that ionized the neutral hydrogen of the intergalactic
medium) is the brightness fluctuation of the 21-cm line of neutral hydrogen as
observed in emission or absorption against the cosmic microwave background
radiation (CMB; Hogan & Rees 1979; Madau et al. 1997; Shaver et al. 1999;
Tozzi et al. 2000; Ciardi & Madau 2003; Zaroubi 2013). The 21-cm signal is a
highly-forbidden hyperfine transition produced by the spin-flip transition of
the ground state of neutral hydrogen. Thanks to the long lifetime of the
excited state combined with the abundance of neutral hydrogen in the
intergalactic medium (IGM) at the inception of the EoR, it can be used to
trace the gas progressive ionization during the EoR. A statistical detection
of the strength of the fluctuations of its brightness temperature can allow us
to constrain our models of the early Universe and the formation of the first
stars and galaxies.

Multiple interferometric low-frequency radio telescopes have been designed to
search for this signal, like LOFAR111Low-Frequency Array,
http://www.lofar.org, HERA222Hydrogen Epoch of Reionization Array,
https://reionization.org/, MWA333Murchison Widefield Array,
http://www.mwatelescope.org, PAPER444Precision Array to Probe EoR,
http://eor.berkeley.edu and the upcoming SKA 555Square Kilometre Array,
https://www.skao.int/en. While the signal is yet to be detected, upper limits
of the 21-cm power spectrum have been obtained (for e.g., Mertens et al.,
2020; The HERA Collaboration et al., 2022; Trott et al., 2020). These have
allowed to rule out some extreme astrophysical models by comparison with
simulations (e.g., Ghara et al., 2020; Ghara et al., 2021; Greig et al.,
2021b; Greig et al., 2021a; Abdurashidova et al., 2022). Further efforts are
being made to minimise observational errors by improving signal extraction
techniques (e.g. Mertens et al., 2023; Acharya et al., 2023). It is therefore
imperative to concurrently improve numerical modeling, readying it for the
detection of the 21cm signal. As of now, the main limitation is the trade-off
between simulated volume and resolution achieved, forced by limits in
computational power.

In fact, as the thermal noise is proportional to the wave-modes (Koopmans et
al., 2015), the best results from upcoming surveys with LOFAR, HERA, MWA and,
subsequently, SKA are expected at the lowest $k$ values, i.e., 0.15 $h~{}{\rm
cMpc}^{-1}\lesssim k\lesssim 1.5~{}h~{}\rm cMpc^{-1}$, after sufficient
mitigation of foreground contamination. This corresponds to large physical
scales, and indeed e.g. Iliev et al. (2014) found that simulations with box
sizes greater than 100 $h^{-1}$ cMpc are necessary to accurately model the
evolution and distribution of ionized regions. However, a simulation with a
box size corresponding to $k\sim 0.15~{}h~{}\rm cMpc^{-1}$ is still limited by
sample variance. The straightforward solution to this issue is to employ
larger boxes in order to increase the sampling of the long wave modes
observationally relevant. For example, in Ghara et al. (2020) a side length of
$500~{}h^{-1}~{}\rm cMpc$ has been simulated, while the semi-numerical
approach of Mesinger et al. (2016) and Greig et al. (2022) has modeled Gpc
scales. Such large boxes also allow to study the topology of large ionized
regions ($\gtrsim 10$ cMpc), which is especially important to understand the
role of quasars. However, even these approaches have issues, as extremely
large computational resources for running the simulations and storing the
generated data are required, and running large number of such simulations
would only be computationally affordable with the semi-numerical approaches.
While they can be appropriate for some analyses, they also have their own
limitations (see for e.g., Choudhury & Paranjape, 2018, for details).

Additionally, there are other limitations on the extent of physics that can be
explored with such simulations. For example, mini halos ($\lesssim
10^{8}~{}\rm M_{\rm\odot}$) are expected to play a significant role in
reionization (see Iliev et al., 2005; Haiman et al., 2001), but their
inclusion would either require a high enough mass resolution to be able to
capture them in large boxes, or running significantly smaller boxes. The
former exacerbates the computational requirements, while the latter limits the
wave modes that can be sampled. A possible approach to tackle these problems
is to use simulations with different box sizes, and combine their results. A
recent example of this is an analysis of galaxy populations at $z\geq 8$ in
Kannan et al. (2023), which, however, does not provide a truly combined
picture, as the different boxes are manually re-scaled to account for
resolution differences without actual convergence between simulations. There
have also been suggestions to use techniques like deep learning to increase
the resolution of large box simulations (see Kodi Ramanah et al. 2020 and Li
et al. 2021 for applications to dark matter only $N$-body simulations),
however these haven’t been carefully tested for EoR redshifts.

Another approach is to run computationally-cheaper boxes but include the
properties of the more-expensive, full-physics ones. Recently, Sinigaglia et
al. (2021) developed a method to map baryonic properties of the IGM onto DM-
only simulations, but its applicability to the reionization epoch is unclear.

An alternative approach is to simulate the full physics of reionization in
small boxes that are augmented with techniques that compensate the limited
volume. One way to do this is the Fixed and Paired (F&P hereafter) approach
described in Angulo & Pontzen (2016) and Pontzen et al. (2016). It has been
shown to substantially reduce the cosmic variance in the matter power spectrum
(see Chartier et al., 2021; Klypin et al., 2020; Villaescusa-Navarro et al.,
2018, for examples) with respect to the traditional approach for the same box
size, bypassing the need to run a large number of smaller boxes. However, this
reduction is substantial only for statistical properties (e.g., the stellar
mass function), and not for local ones like the gas distribution of individual
galaxies (check Villaescusa-Navarro et al., 2018, for more details). The 21-cm
signal is influenced by both the large-scale distribution of neutral hydrogen,
and the local properties of the ionized regions around galaxies. Exploring the
extent of improvement (if any) on the effective volume for estimating the
summary statistics of the 21-cm signal using the F&P approach can be helpful
in minimising the required computational needs. This idea was explored in Giri
et al. (2023) by running a large number of realizations of the F&P approach
using 21cmFAST (Mesinger & Furlanetto, 2007; Greig & Mesinger, 2015). By
comparing them with randomly generated simulations, they found that F&P boxes
could obtain the same precision in the 21-cm signal power spectrum as
traditional boxes twice their size, allowing at least a factor of 4 reduction
in computing costs. However, 21cmFAST does not take into account baryonic
hydrodynamics, thus implicitly assuming that baryons track dark matter. It
also does not include a proper implementation of radiative transfer (RT) or
galaxy properties, and thus becomes unreliable at scales below 1 cMpc. Hence
it cannot provide an accurate picture of galaxy-scale effects on the 21-cm
signal.

To increase confidence in the applicability of the F&P method for improving
EoR 21-cm signal studies, it is necessary to use a more realistic framework,
i.e. one which includes baryonic hydrodynamics, RT and which models galaxy
properties more accurately. For this, we employ a setup similar to the THESAN
simulations (Kannan et al., 2022; Garaldi et al., 2022; Smith et al., 2022;
Garaldi et al., 2023), detailed in Section 2. We discuss the results of
various summary statistics in Section 3, and the improvement on using the F&P
approach in Section 4, while we give our conclusions in Section 5.

##  2 Methodology

###  2.1 The 21-cm signal

The brightness temperature fluctuations of the 21-cm signal are given relative
to the CMB temperature for any patch of the IGM as (see Furlanetto et al.,
2006):

$\delta T_{\rm b}=27x_{\rm HI}(1+\delta_{\rm B})\left(1-\frac{T_{\rm
CMB}}{T_{S}}\right)\\\ \times\biggl{[}\left(\frac{\Omega_{\rm
B}h^{2}}{0.023}\right)\left(\frac{0.15}{\Omega_{\rm
m}h^{2}}\frac{1+z}{10}\right)^{1/2}\biggr{]}~{}\rm mK$ (1)

where $x_{\rm HI}$ is the fraction of neutral hydrogen, $\delta_{\rm B}$ is
the fractional overdensity of baryons, $T_{\rm S}$ is the hydrogen spin
temperature, $T_{\rm CMB}$ is the temperature of the CMB photons at redshift
$z$, $\Omega_{\rm m}$ is the total matter density, $\Omega_{\rm B}$ is the
baryon density, and $h$ is the Hubble constant in units of 100 kms-1cMpc-1.

In the next section, we introduce the simulation used to generate mock
differential brightness temperature ($\delta T_{\rm b}$) maps.

###  2.2 Simulations

Our setup is inspired by the THESAN simulations. We run a suite of radiation-
magneto-hydrodynamic simulations that utilize the moving-mesh hydrodynamics
code AREPO (Springel, 2010; Weinberger et al., 2020), which includes a gravity
solver based on the hybrid Tree-PM method (Barnes & Hut, 1986), a quasi-
Lagrangian Godunov method (Godunov & Bohachevsky, 1959) based hydrodynamics
solver implemented on an unstructured Voronoi mesh grid (Vogelsberger et al.,
2020) and the radiative transfer extension AREPO-RT (Kannan et al., 2019) for
a self-consistent treatment of ionizing radiation. We include the production
and propagation of ionizing photons in three energy bins relevant for hydrogen
and helium photoionization ([13.6, 24.6, 54.4, $\infty$] eV). Further, we
emply a non-equilibrium thermochemistry solver to model the coupling of
radiation fields to gas. For the luminosity and spectral energy density of
stars, we use a complex function of age and metallicity calculated using the
Binary Population and Spectral Synthesis models (BPASS v2.2.1; Eldridge et
al., 2017), modeling the unresolved birth cloud with a uniform escape fraction
of $f_{\rm esc}$ = 0.37. We note that we do not perform a recalibration of
this parameter with respect to THESAN since (i) it mostly impacts the final
phases of reionization, while we are interested in the initial ones, and (ii)
our goal is to compare methods for initial condition generation, so a
slightly-inaccurate reionization history is not expected to affect at all our
results. For further details on the THESAN simulations, see Kannan et al.
(2019); Kannan et al. (2022).

All our simulations have a box-size of $L=95.5~{}\rm cMpc$, and $N=2\times
525^{3}$ particles, giving a dark matter and baryonic particle mass of $m_{\rm
DM}=2.0\times 10^{8}~{}\rm M_{\odot}$ and $m_{\rm gas}=4.7\times 10^{6}~{}\rm
M_{\odot}$, respectively. The gravitational softening length for the star and
dark matter particles is set to 6.0 ckpc, which is also the minimum value for
the adaptively softened gas cells according to cell radius. The cosmological
parameters are taken from Planck Collaboration et al. (2016) as $h=0.6774$,
$\Omega_{\rm m}=0.3089$, $\Omega_{\Lambda}=0.6911$, $\Omega_{\rm b}=0.0486$,
$\sigma_{\rm 8}=0.8159$ and $n_{s}=0.9667$. In the initial conditions the gas
is assumed to follow the DM distribution perfectly, with a primordial hydrogen
and helium fractions of $X=0.76$ and $Y=0.24$, respectively. We start the
simulations from $z_{\rm ini}=49$, generating 54 snapshots between $z=20~{}\rm
to~{}7$. As no further improvement when using the F&P approach is expected
once a significant ionization of neutral hydrogen has been reached, we only
run the simulations down to $z_{\rm fin}=7$, which corresponds to $\langle
x_{\rm HI}\rangle\approx 0.8$.

Lastly, as done in Kannan et al. (2022), we also save Cartesian data output at
a higher cadence (243 outputs between $z=16$ to $z=7$) by gridding the
simulation data onto a regular Cartesian grid employing a first order
Particle-In-Cell approach. We use a $256^{3}$ grid, i.e. each cell is $\sim
372$ ckpc in size.

Figure 1: Evolution of the volume averaged neutral hydrogen fraction versus
redshift for GIC simulations (grey solid), their ensemble average (black
solid) and five F&P averages (orange, magenta, red, blue and green dashed). We
note that at $z$=7 $\langle x_{\rm HI}\rangle$ ranges between 0.72 and 0.82
for the GIC simulations.

####  2.2.1 The Fixed & Paired approach

Figure 2: Middle slices of the $\delta T_{\rm b}$ maps of (from left to
right): one of the fixed simulations, its corresponding pair, their average,
and the average of 2 GIC simulations, at $z=10,8.3,7.6,7$ (with $\langle
x_{\rm HI}\rangle=0.99,0.95,0.90,0.80$ respectively). The average of 2 GIC
simulations shows clear regions of high (and low) $\delta T_{\rm b}$, unlike
the F&P average. The differences between the two averages reduce with
decreasing redshift.

The F&P approach was first proposed by Angulo & Pontzen (2016) and is based on
the creation of a special pair of initial conditions (ICs) that, when employed
together, significantly reduce the impact of cosmic variance. In the
traditional approach to ICs creation, an initially uniform and isotropic
distribution of tracers particles is perturbed (following e.g. Zel’dovich,
1970) to induce density perturbations:

$\delta(\mathbf{k},z_{\rm ini})=\sqrt{P(\mathbf{k},z_{\rm
ini})}e^{i\theta_{\mathbf{k}}}$ (2)

where the phase $\theta_{\mathbf{k}}$ is sampled from a flat distribution in
the range $[0,2\pi]$, and the power spectrum amplitude $P(\mathbf{k},z_{\rm
ini})$ is sampled from a Gaussian distribution centered on its expectation
value $E[P(\mathbf{\hat{k}},z_{\rm ini})]$ (necessary to produce a Gaussian
random field).

In the F&P approach, the power spectrum modes are fixed to their expectation
values. This produces the fixed initial conditions, while the paired one is
obtained by reversing the phase associated to each particle displacement. By
combining simulations run with these two sets of ICs, the variance-induced
fluctuations (which in traditional ICs mainly affect the large-scale modes,
where the sampling of the power spectrum amplitude is scarce and therefore
more susceptible to deviations from its expectation value) are suppressed up
to the fourth perturbative order (Angulo & Pontzen, 2016). Additionally,
despite breaking the Gaussianity of the generated field, this approach does
not induce unwanted features, as shown by e.g., Angulo & Pontzen (2016) and
Chartier et al. (2021) for the matter power spectrum, bispectrum and the halo
mass function, by Anderson et al. (2019) for Lyman-$\alpha$ forest power
spectra, and for a variety of other quantities by Villaescusa-Navarro et al.
(2018) and Klypin et al. (2020).

In this work, we generate 5 such pairs of F&P simulations, in order to explore
the effects of the random sampling of density perturbations. We compare the
averages of these 5 pairs against 35 traditional Gaussian-sampled initial
conditions based simulations (hereafter referred to as GIC). As an example, in
Figure 1 we show the reionization histories of the 5 pairs (magenta, orange,
red, blue, green dashed), as well as that of the GIC simulations (grey solid).
From the figure, we note that the reionization histories of the F&P averages
cluster close to the average of those of the GIC simulations. This is a
consequence of the fact that the F&P method ensures that the F&P averages
closely match the halo mass function, while GIC simulations can spuriously
have an excess/dearth of very bright sources.

In Figure 2, we show maps of $\delta T_{\rm b}$ at different redshifts
($z=10,8.3,7.6,7$; these correspond to $\langle x_{\rm
HI}\rangle=0.99,0.95,0.9,0.8$) for one of the fixed simulations, its
corresponding pair, their average, and the average of two randomly chosen GIC
simulations. Note that due to the phase inversion used for generating the
initial conditions of the pair, the regions of high $\delta T_{\rm b}$ in the
fixed simulation roughly overlaps with regions of low $\delta T_{\rm b}$ in
the pair. This is evident in their average, which does not have specific
regions of high (or low) $\delta T_{\rm b}$, unlike the average of the two GIC
simulations. This difference reduces with decreasing redshift.

In the next section, we compare the ensemble average of the GIC simulations
(taken to be the “true" value) against the F&P average simulations. We show
qualitative comparisons between two of the F&P averages and the true value in
Section 3, and a more quantitative analysis using all five of the averaged
simulations in Section 4.

##  3 Analysis and Results

In this section, we compare the various summary statistics of the 21-cm signal
for the GIC simulations and their ensemble average versus the F&P averages. In
particular, we focus on the skewness, the power spectrum, and the bispectrum.
In Sections 3.1, 3.2 and 3.3 we analyse their behaviour at various redshifts,
focusing on $z=10,9,8,~{}\rm and~{}7$. We also compare the ionized bubble size
distribution across these redshifts in Section 3.4.

###  3.1 Skewness

The skewness ($\tilde{\mu}_{\rm 3,b}$) is a statistical measure of asymmetry
in a distribution, i.e., as the name suggests, it quantifies how skewed a
distribution is around its mean value. $\tilde{\mu}$ = 0 indicates a symmetric
distribution. Here, we use the definition of skewness from Ma et al. (2021),
given as:

$\tilde{\mu}_{\rm 3,b}=\frac{\mu_{3}(\delta T_{\rm b})}{\mu_{2}(\delta T_{\rm
b})^{3/2}}=\rm E\Bigg{[}\frac{(\delta T_{\rm b}-\langle\delta T_{\rm
b}\rangle)^{3}}{\sigma_{\delta T_{\rm b}}^{3}}\Bigg{]}$ (3)

where $\mu_{i}$ is the $i$-th central moment, and $\rm E[]$ is the expectation
value at each snapshot where we calculate the skewness. In our case, this
would be the volume average.

In Figure 3, we compare the evolution of skewness versus redshift from
$z=10~{}\rm to~{}7$ for the 35 GIC simulations, their ensemble average, and
the 5 F&P average simulations. We note that the F&P averages are a good
estimate for the true value of the skewness in most cases. For the 2nd F&P
average, the deviation is larger at $z\lesssim 8$, which is indicative of more
inhomogeneity in the $\delta T_{\rm b}$. We find that this anomalous skewness
is caused by the fixed simulation in the second pair having a chance
association of galaxies with strong black hole feedback and high output of
ionizing photons. While such cases are rare, they are physically possible and
highlight how the relevance of galactic processes in the production of the
21-cm signal hinders the improvements granted by the F&P approach. We also
note that such cases are not completely captured by e.g. the 21cmFAST
simulation used in Giri et al. (2023).

This, however, is not expected to cause a major effect on summary statistics
like the power spectrum and bispectrum, which we explore in the next sections.

Figure 3: Redshift evolution of skewness for 35 GIC simulations (grey solid),
their ensemble average (black solid), and the five F&P averages (orange,
magenta, red, blue and green; dashed).

###  3.2 Power Spectrum

Figure 4:  Top: Power spectra of GIC simulations (grey solid), their ensemble
average (black solid), two F&P pairs (purple and green, dashed and dotted) and
their averages (magenta and orange, solid) for $z=10,8.3,7.6,~{}\rm and~{}7$
from left to right. Bottom: The normalised deviation of each random simulation
and P&F average from the ensemble average.

The power spectrum is expected to be the first detectable statistic of the
21-cm signal during the EoR, and thus is of particular interest. So analysing
the improvement in modelling, if any, with the F&P method is important. The
power spectrum of $\delta T_{\rm b}$ as defined in Equation 1 is given by:

$P_{\rm
21cm}(\mathbf{k})=\delta_{D}(\mathbf{k}+\mathbf{k}^{\prime})\langle\delta
T_{\rm b}(\mathbf{k})\delta T_{\rm b}(\mathbf{k}^{\prime})\rangle$ (4)

where $\delta_{D}$ is the Dirac function, $\delta T_{\rm b}(\mathbf{k})$ is
the differential brightness temperature in Fourier space, and
$\langle...\rangle$ is the ensemble average. In this work, we use the
normalized form of the power spectrum given by:

$\Delta^{2}_{\rm 21cm}=\frac{k^{3}}{2\pi^{2}}\times P_{\rm 21cm}.$ (5)

Figure 5:  Top panels: Equilateral triangle bispectra ($B_{\rm equi}$) for
$0.15h/{\rm cMpc}\leq k\leq 2.0h/\rm cMpc$ of GIC simulations (grey solid),
their ensemble average (black solid), two F&P pairs (purple and green, dashed
and dotted) and their averages (magenta and orange, solid) for
$z=10,8.3,7.6,~{}\rm and~{}7$, clockwise from upper left. Bottom panels: The
normalised deviation of each random simulation and P&F average from the
ensemble average.

In Figure 4 we show the power spectra of the GIC simulations (in grey), their
ensemble average (in black), and the F&P averages (magenta and orange) for 2
pairs (purple and green, with the fixed simulation in dashed, and their
corresponding pairs in dotted lines) out of the 5 generated in Section 2.2.1.
While we have analysed results for all of the five generated F&P pairs (see
Section 4.1 for a quantitative comparison), for the sake of visual clarity we
show only two out of the five pairs.

We note that the F&P averages match the ensemble average more closely for
$z=10~{}\rm and~{}8.3$ than any random GIC simulation, however they become
more comparable to them at $z=7.6~{}\rm and~{}7$. This is understandable as a
combination of two effects. First, the F&P approach is expected to yield
improvements on the predictions of features that are primarily governed by the
large-scale structure, while it does not provide a statistical improvement on
the prediction of galaxy properties (Villaescusa-Navarro et al., 2018), such
as ionizing photons output, as they are entirely dominated by local physics.
Secondly, the period of emergence of ionized regions is one where the 21-cm
signal becomes increasingly non-Gaussian. Thus, the information content in the
power spectrum is reduced in this period as compared to earlier and later
redshifts, which have a homogenized distribution of neutral and ionized
hdyrogen, respectively. So, any possible improvement in computational expense
obtained by using the F&P method rather than a traditional IC simulation due
to the contribution of the large-scale structure to the 21-cm signal is not
going to be visible in the Gaussian components of the signal. For this reason,
we observe an increase in deviation of the F&P averages from the ensemble
average with decreasing redshift (see lower panels of Figure 4). While we
visually note that the F&P averages match more closely the ensemble average as
compared to individual GIC simulations, we quantify the improvement in Section
4.1, where we consider all our F&P pair averages. In contrast, Giri et al.
(2023) reported an improvement by a factor of at least 4 at $z=9$, when
$\langle x_{\rm HI}\rangle$=0.8.

While discrepancies between our results and Giri et al. (2023) could possibly
be due to the difference in methodologies adopted for comparing the F&P method
versus traditional generation of initial conditions, contributions from better
handling of galaxy-driven physics in our simulations are also possible. It is
however difficult to discern the extent of the effect galaxy-driven physics
has on our results.

###  3.3 Bispectrum

Figure 6:  Top panels: Isosceles triangle bispectra ($B_{\rm isoc}$) for
$k_{1}=k_{2}=0.2~{}h~{}\rm cMpc^{-1}$ for $z=10,8.3,7.6,~{}\rm and~{}7$,
clockwise from upper left. Colours and linestyles follow Figure 5. Bottom
panels: The normalised deviation of each random simulation and P&F average
from the ensemble average.

As discussed in the previous section, the non-Gaussianity of the $\delta
T_{\rm b}$ increases with the growth of ionized regions as redshift decreases.
This means that a statistic that focuses only on the Gaussian parts of the
signal, like the power spectrum, encapsulates less and less information as we
move to lower redshift, into the regime of $\langle x_{\rm HI}\rangle\lesssim
0.9$. Therefore, we turn to higher-order statistics to assess whether the F&P
average is still a good approximation of the ensemble average of simulations
in this regime, since they are able to capture the non-Gaussian features of
the $\delta T_{\rm b}$. While skewness as discussed in Section 3.1 would show
some broad non-Gaussian features of the signal, it is still a one point
statistic, i.e. it will not quantify the correlation of the signal between
different Fourier modes. Thus, we now focus on the bispectrum as defined in
Majumdar et al. (2018):

$b_{\rm
21cm}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})=\delta_{D}(\mathbf{k_{1}}+\mathbf{k_{2}}+\mathbf{k_{3}})\langle\delta
T_{\rm b}(\mathbf{k_{1}})\delta T_{\rm b}(\mathbf{k_{2}})\delta T_{\rm
b}(\mathbf{k_{3}})\rangle$ (6)

where $\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}}$ are the Fourier space
wave numbers, $\delta_{D}$ is the Dirac delta function and $\langle\delta
T_{\rm b}(\mathbf{k_{1}})\delta T_{\rm b}(\mathbf{k_{2}})\delta T_{\rm
b}(\mathbf{k_{3}})\rangle$ is a measure of the number of triangles (weighted
by the $\delta T_{\rm b}$ values at their vertices) of different
configurations formed by wave numbers $\mathbf{k_{1}},\mathbf{k_{2}}~{}\rm
and~{}\mathbf{k_{3}}$. The different triangles can be formed by varying the
magnitude of the wave numbers.

To evaluate the bispectra from our simulations, we use the publicly available
BiFFT package (Watkinson et al., 2017), which employs a Fourier transform
based technique (as described in Scoccimarro, 2015; Sefusatti et al., 2016)
much faster rather than the more traditional approach of counting individual
triangles, while still providing consistent results. Following Majumdar et al.
(2020), we normalize $b_{\rm
21cm}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})$ as $B_{\rm
21cm}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})=k_{2}^{3}k_{3}^{3}/(2\pi^{2})^{2}b_{\rm
21cm}(\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}})$.

As done in Figure 4, we focus on wave numbers between $0.15\leq k\leq 2h~{}\rm
cMpc^{-1}$ and show only two of the five F&P averages for the sake of visual
clarity. Further, we consider only two reference cases:

  * •

Equilateral triangles ($B_{\rm equi}$): Here we set $k_{1}=k_{2}=k_{3}=k$,
where $k$ goes from 0.15 to $2~{}h~{}\rm cMpc^{-1}$. This allows us to explore
the non-Gaussian features of the signal across various physical scales. The
results are shown in Figure 5, with the same colours and linestyles as Figure
4. We note that the F&P averages are a close match to the ensemble average
across all redshifts. The apparent large deviations seen at some wave modes of
the F&P average in comparison to the ensemble average (lower panels of Figure
5) arise because at those scales the bispectra approach 0. The normalisation
by the ensemble average thus exaggerates the small differences significantly.
Giri et al. (2023) carried out a similar analysis for $z=9$ (corresponding to
$\langle x_{\rm HI}\rangle=0.8$ in their simulations), and found the F&P
averages to be a close match to the ensemble average for $k>0.1~{}h~{}\rm
cMpc^{-1}$. Thus our results are consistent with their conclusions.

  * •

Isosceles triangles ($B_{\rm isoc}$): We set $k_{1}=k_{2}=k=0.2~{}h~{}\rm
cMpc^{-1}\neq k_{3}$ to explore large physical scales more thoroughly. We plot
$B_{\rm isoc}$ versus the opening angle between the vectors $\mathbf{k}_{1}$
and $\mathbf{k}_{2}$ given as
$\theta=\mathrm{cos}^{-1}(\mathbf{k}_{1}\cdot\mathbf{k}_{2}/(k_{1}k_{2}))$ in
Figure 6, with the same colours and linestyles as Figure 4. As expected, the
GIC simulations show large sample variance, which reduces with decreasing
redshift. However, the F&P averages are an even closer match to the ensemble
average as compared to the GIC simulations, and thus continue to provide an
improvement across all four of the redshift bins used. Similarly to what
observed for $B_{\rm equi}$, the large deviations seen at some scales at
$z=10$ are due to $B_{\rm isoc}$ approaching 0.

We also repeat the process for other $B_{\rm isoc}$, by varying
$k_{1}=k_{2}=k$ from 0.15 to 1.5 $h~{}\rm cMpc^{-1}$, and find that the trends
hold across all physical scales. This is an interesting result, as it provides
a useful statistic for modelling the 21-cm signal using F&P averages at
redshifts with $\langle x_{\rm HI}\rangle\leq 0.9$. It thus works well for
comparing with observations too, as at these redshifts the bispectrum is a
more useful statistic than the power spectrum.

However, it is necessary to verify that the improvement in estimating the
“true" bispectrum when using an F&P average indeed correlates with statistical
estimates of the physical properties (i.e., the distribution of neutral
hydrogen) at the observed redshifts. A good way to check this, is to analyse
the properties and distribution of the source of the non-Gaussian features of
the $\delta T_{\rm b}$, i.e., the ionized regions.

###  3.4 Bubble size distribution

Figure 7: Number of ionized bubbles in different radii $r_{\rm eff}$ bins at
$z=10,8.3,7.6,~{}\rm and~{}7$. We plot their distribution for the GIC
simulations (cyan violin), the ensemble average (black circles), and the two
F&P averages used in Sections 3.2 and 3.3 (magenta and orange triangles).

There is no universal consensus on how to identify ionized regions, with
several methods having been used in literature (see Giri et al., 2018, for a
detailed comparison of different methods). Here, firstly, we choose to define
cells with $x_{\rm HI}\lesssim 0.5$ as ionized. Next, we use a Friend-of-
Friends algorithm based on the ndimage package of SciPy (Virtanen et al.,
2020) to identify regions with clusters of such ionized cells, which are then
considered as “bubbles”. From their volume, we derive the radius of an
equivalent sphere and use it as an effective bubble radius ($r_{\rm eff}$).
The threshold value of $x_{\rm HI}$=0.5 chosen to identify a cell as neutral
or ionized is arbitrary. However, we find that varying this value does not
affect our qualitative results.

The maximum value of $r_{\rm eff}$ is determined by the box size, while we
ignore bubbles equivalent to an individual cell as they are resolution
limited. We then find radii in the range $0.16\lesssim r_{\rm
eff}/(h^{-1}~{}\rm cMpc)\lesssim 40$. As our simulations are run only down to
$z=7$ and quasars are not included as sources of ionizing photons, we do not
expect many bubbles with $r_{\rm eff}>5~{}h^{-1}~{}\rm cMpc$, and thus do not
need to worry about the other extreme of the resolution range. We thus limit
our analysis to the range $0.3\leq r_{\rm eff}/(h^{-1}~{}\rm cMpc)\leq 5$, and
build histograms of number of bubbles ($n_{\rm bubbles}$) binned according to
$r_{\rm eff}$. In Figure 7, we construct violin plots (cyan) for the different
values of $n_{\rm bubbles}$ for the GIC simulations, and also show the
ensemble average (black circles). To compare this with the two F&P averages
used in Section 3.3, we plot the averages for the number of bubbles for the
individual fixed and paired simulations using the same colour scheme (magenta
and orange triangles). Again, we only show two out of the generated five
averages for the sake of easy visual comparison, as all five give similar
results. We note that both F&P averages are a good match for the ensemble
average, and even when they deviate, they remain well within the violins,
showcasing the ranges of the GIC simulations. This confirms that the
improvement in the bispectrum noted when using the F&P average is obtained
because they closely match the number and sizes of ionized regions of the
ensemble average.

Further, as expected, we see that the violin plots get narrower as we approach
lower redshifts, as small scale variability between the GIC simulations grows
with decreasing redshift. Lastly, we note that the number of smaller bubbles
grows when going from $z=10$ to $z=7.6$, but by $z=7$, many of them would have
begun merging, leading to a fall in the number of bubbles with $1.0\leq r_{\rm
eff}/(h^{-1}~{}\rm cMpc)\leq 2.0$.

##  4 Discussion

###  4.1 Advantage of the F&P method

Figure 8:  The standard deviation curve (grey) generated by interpolating the
standard deviations of $\langle P_{\rm ws}\rangle_{\mathrm{m~{}GIC}}$ (grey
squares) versus number of sampled GIC simulations $m$ going from 1 to 35 for
$z=10,8.3,7.6~{}\rm and~{}7$, clockwise from top left. The five F&P averages
are plotted as triangles (circles) with the assumption of them being
1-$\sigma$ (2-$\sigma$) away from $\langle P_{\rm
ws}\rangle_{\mathrm{35~{}GIC}}$ using the same colours as used in Figures 1
and 3. Their closest matches in the standard deviation curve are shown with
black dotted vertical lines gives $m_{\rm eq}$.

Qualitatively, we note that the F&P average provides a closer estimate of the
ensemble average, as compared to any individual GIC simulation, for statistics
like the skewness and the power spectrum, at least for $\langle x_{\rm
HI}\rangle\geq 0.9$. For $\langle x_{\rm HI}\rangle\lesssim 0.9$, the
improvements for the power spectrum is reduced and the bispectrum becomes the
better statistic to use.

However, it is necessary to quantify this improvement as compared to the
average of multiple GIC simulations. For this reason, here we lay out the
methodology to find the number of GIC simulations needed to match the
performance of one F&P average with respect to the power spectrum. For this,
we define “performance” as the extent of deviation of the average of multiple
GIC simulations or of an F&P pair from the power spectrum of the ensemble
average. If the number of GIC simulations required to match the performance of
an F&P average is greater than 2, this means more of them need to be run to
achieve the same performance. In this case, using the F&P average which just
needs 2 simulations to be run would reduce computational costs.

Ideally, for this comparison one should run a large number of GIC simulations
as well as F&P averages, and compare the extent of their deviation from the
ensemble average at specific wave-modes. However, as radiation-hydrodynamical
simulations are computationally expensive, we utilise the 35 GIC simulations
discussed in previous sections, and compare them to all 5 F&P averages
generated in Section 2.2.1. To mimic large number statistics, we proceeds as
follows:

  1. 1.

Wave-mode window: The availability of a large number of simulations allows the
investigation of specific wave-modes, while the same cannot be done when we
are limited to just a small number of them as in this case individual wave-
modes are affected as well by random noise. We thus choose a window of wave-
modes and sum all the power spectra at $ks$ contained within it for each GIC
simulation ($P_{\rm ws}$), their ensemble average ($\langle P_{\rm
ws}\rangle_{\mathrm{35~{}GIC}}$), as well as each F&P average ($\langle P_{\rm
ws}\rangle_{\mathrm{F\&P~{}i}}$, where $i$ goes from 1 to 5). We choose this
window to be $0.15\leq k/(h~{}\rm cMpc^{-1})\leq 0.4$ to ensure that the
largest physical scales covered correspond to those at which the current and
next-generation radio telescopes are most sensitive.

  2. 2.

Choosing GIC simulations: To evaluate the ‘effective volume’ of the P&F
simulations, we need to determine the number $m_{\rm eq}$ of GIC ones that
match the statistical power of a single pair of P&F simulations. However,
different subsets of $m<35$ GIC simulations will produce different results,
especially for small values of $m$. Therefore, in the following we consider
all the ${}^{35}C_{m}$ possible subsets, which is the number of combinations
of $m$ objects out of 35. But for computational reasons, we cap the number of
combinations to $10^{4}$.

  3. 3.

Standard deviation curve: For every $m$, we compute the average $P_{\rm ws}$
for all ${}^{35}C_{m}$ combinations. These $\langle P_{\rm
ws}\rangle_{\mathrm{m~{}GIC}}$ values form a Gaussian distribution centred
around the ensemble average value $\langle P_{\rm
ws}\rangle_{\mathrm{35~{}GIC}}$. We measure the standard deviation of this
Gaussian distribution ($\sigma_{\langle P_{\rm
ws}\rangle_{\mathrm{m~{}GIC}}}$) for each $m$. We linearly interpolate these
values to generate the grey curves in Figure 8, with the grey squares
referring to $\sigma_{\langle P_{\rm ws}\rangle_{\mathrm{m~{}GIC}}}$.

  4. 4.

Standard deviation of the F&P runs: Although the above procedure should be
repeated for the F&P runs, this is prohibitively expensive from a
computational point of view. Therefore, we explicitly parameterize our
ignorance of the true width of the distribution of $\langle P_{\rm
ws}\rangle_{\mathrm{F\&P}}$ by assuming that each F&P average lies exactly
1-$\sigma$ away from its center. In other words, we assume that $\langle
P_{\rm ws}\rangle_{\mathrm{F\&P~{}i}}-\langle P_{\rm
ws}\rangle_{\mathrm{35~{}GIC}}$ is a measure of the width of the Gaussian
distribution of $\langle P_{\rm ws}\rangle_{\mathrm{F\&P}}$. This is done
independently for each F&P run, and is shown in Figure 8 with triangles. We
then repeat this procedure but assuming that each F&P average lies exactly
2-$\sigma$ away from the center of the distribution, and show the results with
circles.

  5. 5.

Closest $m$ matching: For the F&P runs, we compare their standard deviation
with the curve generated from the GIC simulations in step (iii) and determine
the closest $m$ value (black dotted vertical lines in Figure 8). This is our
estimate of $m_{\rm eq}$ for each F&P simulation pair.

  6. 6.

Improvement factor ($f_{\rm imp}$): Finally, we define the improvement factor
$f_{\rm imp}=m_{\rm eq}/2$ for all 5 F&P runs. This corresponds to the ratio
between the number of simulation runs (i.e. simulated volumes, since all our
runs have the same box size) necessary with GIC and F&P ICs, for which we
report the minimum, maximum and average values in Table 1.

From the Figure, we note that the extent of improvement provided by the F&P
method reduces with decreasing redshift, as was expected from the qualitative
results of Section 3.2. However, interestingly, this trend seems to stop at
$z=7.6$, with higher values of $f_{\rm imp}$ noted at $z=7$. This indicates
that while the F&P average may be performing worse at lower redshifts, it
still does better than an average of a few GIC simulations. In fact, we note
that $\langle f_{\rm imp}\rangle$ is $\sim 6$ at $z=10$ (thus one F&P average
is better than running 12 GIC simulations), but $\sim 8$ at $z=7$ (equivalent
to 16 GIC simulations) for the case of F&P averages being 1-$\sigma$ away. The
lowest average improvement is 3.5 at $z=7.6$, indicating that at least 7 GIC
simulations are needed to match an F&P average.

As we sum the power spectra across the aforementioned wave-mode window, a
direct comparison between our results and those of Giri et al. (2023) is
difficult. Nevertheless, we note that in the range $10\geq z\geq 7$, we obtain
$f_{\rm imp}\geq 3.5$, which agrees with their result of an improvement of at
least a factor of 4 at $k=0.1~{}h~{}\rm cMpc^{-1}$ at $z=9,~{}\langle x_{\rm
HI}\rangle=0.8$. The worst possible value of $f_{\rm imp}$ is 0.5, which
corresponds to the case when we have a single GIC simulation perform better.
This is expected, as pure randomness does allow such chance events. However,
as the average improvement is above 2, we believe this still supports the use
of F&P averages for modelling the 21-cm signal power spectrum rather than
running GIC simulations.

A similar analysis can be trivially run for the equilateral triangles
bispectrum. We do not show this calculation, as Figure 5 already demonstrates
clearly that the F&P average bispectra are excellent estimates of the ensemble
average across all redshifts as compared to the individual GIC simulations,
with the scaled deviation being $<1$ across most wave modes.

Table 1:  $f_{\rm imp}$ is the factor of improvement on running an F&P average over running multiple GIC simulations. We report the minimum, maximum and average value of this quantity for the two cases discussed in (iv), for the five F&P averages.  $z$ |  $1-\sigma$ away |  $2-\sigma$ away
---|---|---
| $f_{\rm imp,min}$ | $f_{\rm imp,max}$ | $\langle f_{\rm imp}\rangle$ | $f_{\rm imp,min}$ | $f_{\rm imp,max}$ | $\langle f_{\rm imp}\rangle$
10.0 | 2.0 | 15.5 | 6.4 | 6.5 | 17.0 | 10.9
8.3 | 1.0 | 10.5 | 4.7 | 3.0 | 15.0 | 9.0
7.6 | 0.5 | 8.5 | 3.5 | 1.0 | 14.0 | 7.1
7.0 | 1.0 | 17.0 | 8.3 | 2.5 | 17.5 | 11.2
Figure 9: Redshift evolution of the improvement factor in computational
expense, $f_{\rm imp}$, when using the F&P approach. The errorbars show the
range from $f_{\rm imp,min}$ to $f_{\rm imp,max}$ for the F&P averages being
1-$\sigma$ away (yellow) or 2-$\sigma$ away (blue) from an assumed
distribution of F&P averages, and the points show the average, $\langle f_{\rm
imp}\rangle$. The top axis shows the corresponding $\langle x_{\rm HI}\rangle$
and the right axis shows the number $m_{\rm eq}$ of GIC simulations
corresponding to $f_{\rm imp}$. The black dashed line indicates no
improvement.

###  4.2 Limitations and future applications

As we have shown in the previous Section, the F&P method can be a powerful
tool to extend the statistical accuracy of limited-volume simulations of the
EoR 21-cm signal. In combination with the advancement of computational
techniques and hardware, this method can deliver accurate predictions over
volumes of interest for the study of the 21-cm signal, especially in the early
phases of cosmic reionization. However, we caution that the F&P method does
not improve our predictive ability concerning individual galaxy properties, as
they are dominated by local processes and environment. Therefore, the
improvement granted on statistical quantities (e.g. the power spectrum and
bispectrum) does not necessarily mirror exactly on different observables that
might be more affected by individual objects (see e.g. the discussion in Sec.
3.1).

We foresee numerous application of this technique. For instance, the lower
computational cost required (with respect to traditional approaches) entails
e.g. that larger volumes, as well as a broader range of cosmologies and/or
physical models can be explored with accurate simulations, greatly improving
the reliability of predictions and – eventually – inference from 21-cm signal
data. While basing the entire inference process on RMHD simulations remains
prohibitive, they will be needed in order to confirm and improve constraints
obtained through computationally-cheaper less-accurate methods. Additionally,
they are necessary to explore the coupling between small and large scales,
e.g. spatial correlations between galaxies and the 21-cm signal, that the
advent of SKA will enable. Our results suggest that such accurate predictions
can be obtained from RMHD simulations even in statistically-significant
volumes of the Universe.

Finally, the F&P approach operates orthogonally to super-resolution
techniques; therefore these two approaches should be considered complementary
rather than in opposition.

##  5 Summary

Running simulations for EoR is computationally very expensive. This is
exacerbated when aiming at resolving low-mass galaxies (whose importance has
been recently shown observationally in Atek et al., 2023) in volumes large
enough to be statistically significant for reionization studies. In this work,
we explored the Fixed & Paired (F&P) approach (Angulo & Pontzen, 2016; Pontzen
et al., 2016) to investigate the possibility of reducing the number of
simulations needed (and thus the overall computational expense) to produce
unbiased models of the 21-cm signal. While past efforts have used semi-
numerical approaches to implement this, we have shown more rigorous results by
using radiation hydrodynamic simulations that model more accurately galaxy-
scale effects on the 21-cm signal. We focus on the wave modes in the range
$0.15\leq k/(h~{}\rm cMpc^{-1})\leq 2$, as the best measurements from present
and upcoming radio telescopes like LOFAR, HERA, MWA and SKA are expected in
this regime. Further, we focus on redshifts $10\geq z\geq 7$, which in our
case correspond to $1.0>\langle x_{\rm HI}\rangle\geq 0.8$.

To explore the improvement with respect to various 21-cm signal statistics
obtained by adopting the F&P approach rather than running Gaussian-sampled
initial conditions based (GIC) simulations, we use a setup similar to that of
the THESAN project (Kannan et al., 2022; Garaldi et al., 2022; Smith et al.,
2022; Garaldi et al., 2023). In particular, we investigate the impact on the
skewness, power spectrum, bispectrum and the ionized region size distribution,
and introduce a novel method to quantify the improvement.

We find that the skewness and power spectrum are well-estimated by the F&P
averages for $\langle x_{\rm HI}\rangle\geq 0.9$, and their performance is
good also in the range $0.9>\langle x_{\rm HI}\rangle\geq 0.8$, with an
improvement in computational cost better than 3.5 for the generation of the
power spectrum. We find that the bispectrum is well estimated for $\langle
x_{\rm HI}\rangle\geq 0.8$, with the scaled deviation between the F&P averages
and the ensemble average being below $1$. The only exception are those modes
where the ensemble average bispectrum approaches zero, artificially increasing
the small differences between the F&P average and the ensemble average. This
confirms that the F&P average bispectrum is a great complement to the power
spectrum for studies of the 21-cm signal when $\langle x_{\rm HI}\rangle\geq
0.8$. Finally, we show that in this regime the F&P averages also provide a
good estimate of the HII regions size distribution, with the F&P averages
being within $2-\sigma$ deviation of the ensemble average for bubbles with
radius $\leq 5h^{-1}$ cMpc.

Thus, the F&P method can be used to model the 21-cm signal summary statistics
with significantly reduced effective volume and computational expense.

## Acknowledgements

The authors would like to thank Sambit K. Giri and Raul Angulo for useful
comments and discussions.

## Data Availability

The data used in this work can be shared upon reasonable request.

## References

  * Abdurashidova et al. (2022) Abdurashidova Z., et al., 2022, ApJ, 924, 51
  * Acharya et al. (2023) Acharya A., et al., 2023, submitted.
  * Anderson et al. (2019) Anderson L., Pontzen A., Font-Ribera A., Villaescusa-Navarro F., Rogers K. K., Genel S., 2019, ApJ, 871, 144
  * Angulo & Pontzen (2016) Angulo R. E., Pontzen A., 2016, MNRAS, 462, L1
  * Atek et al. (2023) Atek H., et al., 2023, arXiv e-prints, p. arXiv:2308.08540
  * Barnes & Hut (1986) Barnes J., Hut P., 1986, Nature, 324, 446
  * Chartier et al. (2021) Chartier N., Wandelt B., Akrami Y., Villaescusa-Navarro F., 2021, MNRAS, 503, 1897
  * Choudhury & Paranjape (2018) Choudhury T. R., Paranjape A., 2018, MNRAS, 481, 3821
  * Ciardi & Madau (2003) Ciardi B., Madau P., 2003, ApJ, 596, 1
  * Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publ. Astron. Soc. Australia, 34, e058
  * Furlanetto et al. (2006) Furlanetto S. R., Oh S. P., Pierpaoli E., 2006, Phys. Rev. D, 74, 103502
  * Garaldi et al. (2022) Garaldi E., Kannan R., Smith A., Springel V., Pakmor R., Vogelsberger M., Hernquist L., 2022, MNRAS, 512, 4909
  * Garaldi et al. (2023) Garaldi E., et al., 2023, arXiv e-prints, p. arXiv:2309.06475
  * Ghara et al. (2020) Ghara R., et al., 2020, MNRAS, 493, 4728
  * Ghara et al. (2021) Ghara R., Giri S. K., Ciardi B., Mellema G., Zaroubi S., 2021, MNRAS, 503, 4551
  * Giri et al. (2018) Giri S. K., Mellema G., Dixon K. L., Iliev I. T., 2018, MNRAS, 473, 2949
  * Giri et al. (2023) Giri S. K., Schneider A., Maion F., Angulo R. E., 2023, A&A, 669, A6
  * Godunov & Bohachevsky (1959) Godunov S. K., Bohachevsky I., 1959, Matematičeskij sbornik, 47(89), 271
  * Greig & Mesinger (2015) Greig B., Mesinger A., 2015, MNRAS, 449, 4246
  * Greig et al. (2021a) Greig B., Trott C. M., Barry N., Mutch S. J., Pindor B., Webster R. L., Wyithe J. S. B., 2021a, MNRAS, 500, 5322
  * Greig et al. (2021b) Greig B., et al., 2021b, MNRAS, 501, 1
  * Greig et al. (2022) Greig B., Wyithe J. S. B., Murray S. G., Mutch S. J., Trott C. M., 2022, MNRAS, 516, 5588
  * Haiman et al. (2001) Haiman Z., Abel T., Madau P., 2001, ApJ, 551, 599
  * Hogan & Rees (1979) Hogan C. J., Rees M. J., 1979, MNRAS, 188, 791
  * Iliev et al. (2005) Iliev I. T., Shapiro P. R., Raga A. C., 2005, MNRAS, 361, 405
  * Iliev et al. (2014) Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, MNRAS, 439, 725
  * Kannan et al. (2019) Kannan R., Vogelsberger M., Marinacci F., McKinnon R., Pakmor R., Springel V., 2019, MNRAS, 485, 117
  * Kannan et al. (2022) Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M., Hernquist L., 2022, MNRAS, 511, 4005
  * Kannan et al. (2023) Kannan R., et al., 2023, MNRAS, 524, 2594
  * Klypin et al. (2020) Klypin A., Prada F., Byun J., 2020, MNRAS, 496, 3862
  * Kodi Ramanah et al. (2020) Kodi Ramanah D., Charnock T., Villaescusa-Navarro F., Wandelt B. D., 2020, MNRAS, 495, 4227
  * Koopmans et al. (2015) Koopmans L., et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 1 (arXiv:1505.07568), doi:10.22323/1.215.0001
  * Li et al. (2021) Li Y., Ni Y., Croft R. A. C., Di Matteo T., Bird S., Feng Y., 2021, Proceedings of the National Academy of Science, 118, e2022038118
  * Ma et al. (2021) Ma Q.-B., Ciardi B., Eide M. B., Busch P., Mao Y., Zhi Q.-J., 2021, ApJ, 912, 143
  * Madau et al. (1997) Madau P., Meiksin A., Rees M. J., 1997, ApJ, 475, 429
  * Majumdar et al. (2018) Majumdar S., Pritchard J. R., Mondal R., Watkinson C. A., Bharadwaj S., Mellema G., 2018, MNRAS, 476, 4007
  * Majumdar et al. (2020) Majumdar S., Kamran M., Pritchard J. R., Mondal R., Mazumdar A., Bharadwaj S., Mellema G., 2020, MNRAS, 499, 5090
  * Mertens et al. (2020) Mertens F. G., et al., 2020, MNRAS, 493, 1662
  * Mertens et al. (2023) Mertens F. G., Bobin J., Carucci I. P., 2023, arXiv e-prints, p. arXiv:2307.13545
  * Mesinger & Furlanetto (2007) Mesinger A., Furlanetto S., 2007, ApJ, 669, 663
  * Mesinger et al. (2016) Mesinger A., Greig B., Sobacchi E., 2016, MNRAS, 459, 2342
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 596, A108
  * Pontzen et al. (2016) Pontzen A., Slosar A., Roth N., Peiris H. V., 2016, Phys. Rev. D, 93, 103519
  * Scoccimarro (2015) Scoccimarro R., 2015, Phys. Rev. D, 92, 083532
  * Sefusatti et al. (2016) Sefusatti E., Crocce M., Scoccimarro R., Couchman H. M. P., 2016, MNRAS, 460, 3624
  * Shaver et al. (1999) Shaver P. A., Windhorst R. A., Madau P., de Bruyn A. G., 1999, A&A, 345, 380
  * Sinigaglia et al. (2021) Sinigaglia F., Kitaura F.-S., Balaguera-Antolínez A., Nagamine K., Ata M., Shimizu I., Sánchez-Benavente M., 2021, ApJ, 921, 66
  * Smith et al. (2022) Smith A., Kannan R., Garaldi E., Vogelsberger M., Pakmor R., Springel V., Hernquist L., 2022, MNRAS, 512, 3243
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * The HERA Collaboration et al. (2022) The HERA Collaboration et al., 2022, arXiv e-prints, p. arXiv:2210.04912
  * Tozzi et al. (2000) Tozzi P., Madau P., Meiksin A., Rees M. J., 2000, ApJ, 528, 597
  * Trott et al. (2020) Trott C. M., et al., 2020, MNRAS, 493, 4711
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro F., et al., 2018, ApJ, 867, 137
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Reviews Physics, 2, 42
  * Watkinson et al. (2017) Watkinson C. A., Majumdar S., Pritchard J. R., Mondal R., 2017, MNRAS, 472, 2436
  * Weinberger et al. (2020) Weinberger R., Springel V., Pakmor R., 2020, ApJS, 248, 32
  * Zaroubi (2013) Zaroubi S., 2013, The Epoch of Reionization. Springer Berlin Heidelberg, Berlin, Heidelberg, pp 45–101, doi:10.1007/978-3-642-32362-1_2, https://doi.org/10.1007/978-3-642-32362-1_2
  * Zel’dovich (1970) Zel’dovich Y. B., 1970, A&A, 5, 84


# End-to-end View Synthesis via NeRF Attention

Zelin Zhao
The Chinese University of Hong Kong
<EMAIL_ADDRESS>
&Jiaya Jia
The Chinese University of Hong Kong & SmartMore
<EMAIL_ADDRESS>

###### Abstract

In this paper, we present a simple seq2seq formulation for view synthesis
where we take a set of ray points as input and output colors corresponding to
the rays. Directly applying a standard transformer on this seq2seq formulation
has two limitations. First, the standard attention cannot successfully fit the
volumetric rendering procedure, and therefore high-frequency components are
missing in the synthesized views. Second, applying global attention to all
rays and pixels is extremely inefficient. Inspired by the neural radiance
field (NeRF), we propose the NeRF attention (NeRFA) to address the above
problems. On the one hand, NeRFA considers the volumetric rendering equation
as a soft feature modulation procedure. In this way, the feature modulation
enhances the transformers with the NeRF-like inductive bias. On the other
hand, NeRFA performs multi-stage attention to reduce the computational
overhead. Furthermore, the NeRFA model adopts the ray and pixel transformers
to learn the interactions between rays and pixels. NeRFA demonstrates superior
performance over NeRF and NerFormer on four datasets: DeepVoxels, Blender,
LLFF, and CO3D. Besides, NeRFA establishes a new state-of-the-art under two
settings: the single-scene view synthesis and the category-centric novel view
synthesis. The code will be made publicly available.

##  1 Introduction

The view synthesis task is to generate a novel image from an arbitrary
viewpoint [23, 61]. We consider a simple seq2seq formulation: the input
includes ray points (points sampled on rays), and the output is the colors
contributed by the corresponding rays in the image plane (as shown in Figure
1). This formation does not require rendering and does not require the 3D
model of objects. With this simple formulation, we are heading towards a
uniform end-to-end machine learning model [1].

One possible solution for this problem formulation is the transformer [47]
which has become the dominant neural model in many computer vision and natural
language processing tasks [47, 19, 43, 2]. We first try the standard
transformer under the seq2seq formulation. It is found that the vanilla
transformer yields poor performance (refer to Figure 4). We hypothesize the
following reasons for this fact. On the one hand, the standard global
attention lacks the inductive bias of the volumetric rendering process, and
the resulting images tend to be over-smoothed. On the other hand, global
attention is computationally heavy, so the training is extremely inefficient.

To address these issues, we propose the NeRF attention (NeRFA), inspired by
neural radiance fields (NeRF [23]). We turn the volumetric rendering process
into a feature modulation procedure to empower the transformer with the latent
rendering ability. Unlike NeRF, we do not explicitly predict the mid-way
colors and densities and instead conduct latent rendering via the attention
mechanism. Therefore, NeRFA preserves the attention mechanism’s representation
potential while exploiting the inductive bias inspired by the volumetric
rendering equation. Besides, the NeRFA computes ray attention and pixel
attention separately to save computational costs [19].

Importantly, NeRFA is the first to provide a full-attention solution to the
novel view synthesis task. The most related work, NerFormer [35] is published
accompanied by the CO3D dataset. As shown in Figure 1, NerFormer replaces the
MLP in NeRF with the transformers and consumes multi-view input to generate
color and density predictions. It does not eliminate the volumetric rendering
and is of different model architecture from NeRFA.

NeRFA generates photo-realistic images without explicit 3D object models. It
inherits the power of NeRF by learning latent volumetric representation in
continuous 3D space. In addition, NeRFA adopts transformers to learn
interaction between ray points and builds connections regarding colors of
surrounding pixels, which helps generate sharper regions in complex scenes.

We conduct extensive experiments to validate the effectiveness of the NeRFA
model. Visual comparison between NeRF and NeRFA shows that NeRFA produces
better fine-grained details of objects. The quantitative results also manifest
that the NeRFA model consistently outperforms NeRF [23] and NerFormer [35] on
Blender [23], LLFF [22], and the DeepVoxels [39] datasets for the single-scene
view synthesis. Furthermore, NeRFA also surpasses the NeRF [23] and NerFormer
[35] on the CO3D [35] dataset in the category-centric novel view synthesis.

Figure 1: Essential difference between NeRF [23], NerFormer [35] and NeRFA.
NeRF first predicts the color and densities of the ray points in 3D space,
then renders them into a 2D image. NerFormer replaces the MLP with 3D
transformers and builds color and densities from multiple source views. We
instead address the view synthesis problem in a seq2seq manner where the NeRFA
model transfers ray points to pixel colors in the target image.

##  2 Related work

#### View synthesis

Early graphics work, such as the Lumigraph [11], reconstructs novel views via
reconstructed light field. Later work operates on mesh-based representation to
assign color [52] or add texture [48, 58, 57] to 3D objects. In [61], a CNN
predicts the appearance flow from the input to the output view. [51] presents
an end-to-end point cloud-based framework Synsin for single-image view
synthesis. This setting is different from our multi-frame view synthesis task.

Another line of work leverages the volumetric representation to synthesize
high-quality views from observed images. DeepView [7] and the follow-up work
[41, 22] synthesize photo-realistic images via multi-plane images (MPIs). [13,
29, 46] use learnable volumetric rendering to render views. They can hardly
generate photo-realistic natural images.

#### Neural radiance field

The neural radiance field (NeRF [23]) provides a simple yet effective solution
to the view synthesis problem. The NerFormer model is presented along with the
CO3D dataset in [35]. It replaces the MLP module with a transformer encoder
and does not change the ray march module. We compare the NerFormer thoroughly
in the experiments. In [36], a transformer architecture for the single-view
synthesis task was proposed. Their setting is different from ours. The VQGAN
[6] structure is adopted to address the ambiguity issue in single views.
Neural sparse voxel fields (NSVF [18]) use voxel octrees to boost inference in
learning scene representations. Other recent work speeds up training [3] and
inference [17, 10, 34] of NeRF. The work of [27, 30, 8, 32] reconstructs a
non-rigid deformable scene from captured images. Further, GRAF [37] and
Giraffe [26] propose generative networks based on NeRF to synthesize new
images. IBRNet [50] uses a ray transformer to summarize densities from
neighborhood views. Another important application of NeRF is scene
reconstruction [49, 54].

#### Neural rendering

Researchers made extensive efforts to pursue learnable photo-realistic
rendering models. RenderNet [25] is a convolutional neural network learned to
perform projection and shading. The visual object networks [62] learn to
generate realistic images from learned disentangled 3D representations.
Textured neural avatars [38] generate full-body rendering of a person via
image-to-image deep networks. In [45], neural textures improve the rendering
quality. These methods rely on the 3D models of objects and therefore cannot
be applied to general natural images directly. The concurrent work ViewFormer
[16] proposes a NeRF-free neural rendering, which addresses the view synthesis
problem in 2D space, so their synthesis quality is limited.

#### Attention mechanisms and transformers

The attention mechanism was proposed in machine translation [47]. Transformers
then demonstrate power in many different domains, such as cross-modal
reasoning [43, 59], large-scale unsupervised pretraining [4, 42, 53], and
multi-task representation learning [31, 24, 60]. Recently, transformers become
competitive models in fundamental vision tasks, such as image classification
[5, 55], object detection [2], and segmentation [19, 9]. Transformers
demonstrate high generalization [28] ability and serve as foundation models
[1].

##  3 Problem formation and approach

In Sec. 3.1, we present a simple seq2seq formulation for the view synthesis
problem. We then introduce the vanilla attention model in Sec. 3.2.
Subsequently, we introduce the NeRF attention in Sec. 3.3, which addresses key
limitations of the vanilla attention model.

###  3.1 End-to-end view synthesis

As shown in Figure 1, we uniformly sample sequences of points on rays in the
3D space, and we call them ray points (they can also be called query points
[23]). The input to the model is a set of ray points in the 3D space, and the
output is a set of pixel colors corresponding to those ray points. Formally,
the input ray points are represented by a matrix $\mathbf{P}$ of shape
$N_{p}\times N_{r}\times 6$ where $N_{r}$ is the number of ray points
corresponding to one ray, and $N_{p}$ is the number of rays sampled in one
image. Since each ray corresponds to a projected pixel in the image plane,
$N_{p}$ is also the number of sampled pixels. Here one ray point is
represented by a 6D tuple $(x,y,z,d_{x},d_{y},d_{z})$ where $(x,y,z)$ is the
3D location and $(d_{x},d_{y},d_{z})$ is the unit vector representing view
direction. The output is the RGB color $\mathbf{C}\in\mathbb{R}^{N_{p}\times
3}$. Given the ground-true RGB color
$\mathbf{C}^{gt}\in\mathbb{R}^{N_{p}\times 3}$, we can train the model via a
simple L2 loss of
$\mathcal{L}=\left\|\mathbf{C}^{gt}-\mathbf{C}\right\|_{2}^{2}$. This simple
end-to-end framework eliminates the explicit rendering procedure and the
hierarchical volume sampling proposed in NeRF [23].

Figure 2: Illustration of the vanilla attention model (Sec. 3.2) and our
proposed NeRF attention model (Sec. 3.3). Left: The vanilla attention model
simply adopting self-attention, and the pooling and fully connected (FC)
layers to get the predicted color results. Right: NeRFA proposes a feature
modulation procedure to conduct latent rendering. It uses a ray transformer
and a pixel transformer to learn interactions of rays and pixels. A FC layer
then summarizes the results and predicts the final color.

###  3.2 Vanilla attention

Based on the seq2seq framework, we first introduce the vanilla attention
(VaniA) solution. It is illustrated on the left side of Figure 2. The first
step is to transfer the input ray points $\mathbf{P}$ to positional embedding,
essential for view synthesis. The embedder of NeRF [23] is expressed as

$\mathbf{F_{emb}}=\mbox{MLP}(\mathbf{P})\circ\mbox{CosEmbed}(\mathbf{P}),$ (1)

where $\mathbf{F_{emb}}$ denotes the positional embedding matrix with shape
$N_{p}\times N_{r}\times d$ where $d$ is the hidden dimension size.

We then send the embedding matrix $\mathbf{F_{emb}}$ to the transformer [47]
block to produce pixel-ray features $\mathbf{F_{pv}}\in\mathbb{R}^{N_{p}\times
N_{r}\times d}$. The transformer block is denoted as

$\mathbf{F_{pv}}=\mbox{GlobalTransformer}(\mathbf{F_{emb}})=\mbox{SelfAtten}(\mbox{LN}(\mathbf{F_{emb}}))+\mathbf{F_{emb}},$
(2)

where LN is layer normalization and SelfAtten represents the standard self-
attention [47] of

$\mbox{SelfAtten}(\mathbf{X})=\mbox{softmax}\left(\frac{\mathbf{X}\mathbf{X}^{T}}{\sqrt{d}}\right)\mathbf{X}.$
(3)

The transformer contains $L$-layer transformer blocks and each block contains
a $H$-head attention where $L$ and $H$ are two hyper-parameters. We omit the
multi-layer and multi-head implementation of transformers for simplicity in
all equations. We call this transformer in Eq. (2) “global transformer"
because the sequence length is $N_{p}N_{r}$ (that is, the embedding matrix is
reshaped to a two-dimensional tensor of shape $N_{p}N_{r}\times d$). We
accumulate the latent features along with rays via average pooling over the
ray dimension, and use fully connected layers (FC) to obtain the final color
prediction $\mathbf{C}$ as

$\mathbf{C}=\mathbf{FC}\left(\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\mathbf{F_{pv}}[:,i,:]\right).$
(4)

Our experiments (Sec. 4.5) show that vanilla attention does not satisfyingly
model the rendering process because the synthesized view misses a lot of high-
frequency information. This observation aligns with the recent finding of [23,
44]. Besides, the vanilla attention is computational heavy because of the
global transformer. The following section proposes the NeRF attention to
address these key issues.

###  3.3 NeRF attention

#### Positional embedding and ray transformer.

Refer to the right of Fig. 2, our NeRF attention first produces positional
embedding via Eq. (1). After that, the NeRF computes the ray attention to get
the pixel-ray features. Different from Eq. (3), the ray attention is only
computed for each ray independently. We write the pixel-ray feature
$\mathbf{F^{\prime}_{pv}}$ as

$\mathbf{F^{\prime}_{pv}}[p]=\mbox{RayTransformer}(\mathbf{F_{emb}}[p])=\mbox{SelfAtten}(\mbox{LN}(\mathbf{F_{emb}}[p]))+\mathbf{F_{emb}}[p],\forall
p\in[1,N_{p}].$ (5)

The global transformer in vanilla attention (Eq. (1)) has complexity
$O(N_{p}^{2}N_{r}^{2}d)$ since $N_{p}$ is much larger than $d$ [63]. In
contrast, the ray transformer in the NeRF attention (Eq. (3)) reduces the
complexity to $O(N_{p}N_{r}^{2}d)$, making training more efficient.

#### Feature modulation.

Inspired by NeRF [23], we modulate the pixel-ray feature to learn volumetric
rendering. In NeRF [23], volumetric rendering is expressed as

$C(\mathbf{r})=\sum_{i=1}^{N_{r}}\left(e^{-\sum_{j=1}^{i-1}\sigma_{j}\delta_{j}}(1-e^{-\sigma_{i}\delta_{i}})\mathbf{c_{i}}\right),$
(6)

where $C(\mathbf{r})$ is the color corresponding to a traced ray $\mathbf{r}$,
$\sigma_{j}$ denotes the predicted density, $\mathbf{c_{i}}$ is the predicted
color of a ray point, and $\delta_{i}$ is the distance between adjacent
samples. We construct the modulated feature matrix
$\mathbf{F_{mod}}\in\mathbb{R}^{N_{p}\times d}$ from the pixel-ray feature as

$\mathbf{F_{mod}}[p]=\sum_{i=1}^{N_{r}}\left(\mbox{exp}(-\sum_{j=1}^{i-1}\delta_{j}\mathbf{F^{\prime}_{pv}}[p][j])\circ\left(\mathbf{1}-\mbox{exp}(-\delta_{i}\mathbf{F^{\prime}_{pv}}[p][i])\right)\circ\mathbf{F^{\prime}_{pv}}[p][i]\right),$
(7)

where the $\mbox{exp}(\mathbf{v})$ is an element-wise exponential function for
a vector $\mathbf{v}$ and $\circ$ represents the element-wise product
operator. The modulation function in Eq. (7) differs from the NeRF equation
(Eq. (6)) because it does not rely on the explicit prediction of volume
density and color and it conducts rendering in the latent space. Our
experiments show that the modulated features help model the high-frequency
components of the scene.

#### Pixel transformer and output.

We adopt a pixel transformer to enable pixel-wise interactions. The pixel
transformer takes the modulated features $\mathbf{F_{mod}}$ as input, and
models interactions of pixel-level features. This process is formulated as

$\mathbf{F_{color}}=\mbox{PixelTransformer}(\mathbf{F_{mod}})=\mbox{SelfAtten}(\mbox{LN}(\mathbf{F_{mod}}))+\mathbf{F_{mod}},$
(8)

where $\mathbf{F_{color}}\in\mathbb{R}^{N_{p}\times d}$ is the color
embeddings. This operation in Eq. (8) yields complexity $O(N_{p}^{2}d)$.
Finally, we use a fully connected layer to get color
$\mathbf{C}=\mbox{FC}(\mathbf{F_{color}})$. We empirically find that the pixel
transformer contributes to the fine-grained details in the synthesized image
in experiments.

#### Overall complexity analysis.

The overall time complexity of the NeRFA model is
$O(N_{p}N_{r}^{2}d+N_{p}^{2}d)$, which is much lower than the vanilla
attention model $O(N_{p}^{2}N_{r}^{2}d)$. The NeRFA model can be trained
within three days on a single NVIDIA V100 GPU for 300k iterations, comparable
to NeRF. Note that NeRF training needs one or two days according to statistics
in [23]. The speed of NeRFA can be further optimized via recent techniques,
such as the fast NeRF [10] or the efficient transformer [15].

Figure 3: Single-scene view synthesis results on the test split of the Blender
[23] and LLFF [22]. Our approach is compared to SRN [40], NeRF [23], and
NerFormer [35]. We observe that NerFormer would deteriorate the NeRF but our
NeRFA model improves the NeRF and produces sharper details.

##  4 Experiments

In this section, we compare NeRFA to NeRF [23], NerFormer [35] and other
baselines on four datasets. We also conduct ablation studies to verify the
effectiveness of NeRFA components. A demo video of our approach is provided in
the supplementary material.

###  4.1 Datasets

#### DeepVoxels.

The DeepVoxels dataset [39] includes four objects: vase, pedestal, chair, and
cube. The objects are relatively simple in geometry. The training sets are
formed by rendering objects via 479 poses uniformly sampled in the northern
hemisphere. We reserve thirty images for validation. Another 1000 testing
views are rendered on an Archimedean spiral in the northern hemisphere. All
images are of $512\times 512$ resolution.

#### Blender.

The Blender dataset is presented in the original NeRF [23] paper. This dataset
contains eight objects (materials, drums, etc.) with complicated geometry and
non-Lambertian effect. Each object is rendered for 100 training views and 200
testing views. We reserve ten training views for a validation split. The
resolution is $800\times 800$ pixels.

#### LLFF.

The LLFF [22] dataset includes eight rendered scenes in total. All are
captured by a cellphone with $20-62$ images where $\frac{1}{8}$ of the images
are reserved for the test. We reserve another $\frac{1}{8}$ of all images to
form the validation split to tune model hyper-parameters. The resolution of
the rendered frames is $1008\times 756$.

Table 1: Experimental results on the test split of three datasets (DeepVoxels
[39], Blender [23] and LLFF [22]) under the single-scene view synthesis
setting. We highlight the best results in bold.

| DeepVoxels [39] | Blender [23] | LLFF [22]
---|---|---|---
| PSNR$\uparrow$ | SSIM$\uparrow$ | LPIPS$\downarrow$ | PSNR$\uparrow$ | SSIM$\uparrow$ | LPIPS$\downarrow$ | PSNR$\uparrow$ | SSIM$\uparrow$ | LPIPS$\downarrow$
SRN [40] | 33.15 | 0.947 | 0.069 | 22.15 | 0.843 | 0.164 | 22.78 | 0.654 | 0.369
NV [20] | 29.57 | 0.930 | 0.094 | 26.31 | 0.879 | 0.158 | - | - | -
LLFF [22] | 34.21 | 0.986 | 0.045 | 24.95 | 0.933 | 0.113 | 24.14 | 0.801 | 0.204
NeRF [23] | 40.14 | 0.990 | 0.021 | 32.44 | 0.956 | 0.076 | 26.49 | 0.812 | 0.235
NerFormer [35] | 35.51 | 0.975 | 0.030 | 28.53 | 0.941 | 0.095 | 25.61 | 0.821 | 0.218
NeRFA (Ours) | 40.31 | 0.992 | 0.020 | 33.11 | 0.959 | 0.072 | 28.14 | 0.832 | 0.198

Table 2: Experimental results on the CO3D [35] under the category-centric
novel view synthesis setting. Baseline numbers are taken from [35]. The best
results are highlighted in bold.

| Train-unseen | Test-unseen
---|---|---
|  | Average
---
statistics
| PSNR @
---
# source views
| PSNR @
---
target view difficulty
| Average
---
statistics
| PSNR @
---
# source views
| PSNR @
---
target view difficulty
| PSNR$\uparrow$ | LPIPS$\downarrow$ | 9 | 5 | 1 | easy | med. | hard | PSNR$\uparrow$ | LPIPS$\downarrow$ | 9 | 5 | 1 | easy | med. | hard
NV-WCE [20] | 12.3 | 0.34 | 12.5 | 12.3 | 12.0 | 12.4 | 12.0 | 13.6 | 11.6 | 0.35 | 11.7 | 11.6 | 11.6 | 11.7 | 11.2 | 12.0
P3DMesh [33] | 17.2 | 0.23 | 17.6 | 17.4 | 16.2 | 17.5 | 16.8 | 16.0 | 12.4 | 0.26 | 12.6 | 12.5 | 12.1 | 12.6 | 11.8 | 13.2
IPC-WCE [51] | 14.1 | 0.36 | 14.4 | 14.2 | 13.4 | 14.4 | 13.7 | 13.4 | 13.8 | 0.27 | 13.8 | 13.7 | 12.6 | 13.8 | 12.8 | 12.2
SRN-WCE [40] | 16.6 | 0.26 | 17.0 | 16.7 | 15.8 | 16.9 | 15.7 | 14.5 | 14.6 | 0.27 | 14.9 | 14.8 | 13.9 | 14.9 | 13.7 | 14.8
NeRF-WCE [12] | 14.3 | 0.27 | 14.3 | 14.9 | 14.2 | 12.1 | 13.6 | 13.6 | 13.8 | 0.27 | 14.3 | 14.9 | 14.2 | 14.4 | 13.0 | 13.0
NerFormer [35] | 17.9 | 0.26 | 19.3 | 18.3 | 15.6 | 18.9 | 15.5 | 14.6 | 17.6 | 0.27 | 18.9 | 18.1 | 15.1 | 18.6 | 14.9 | 14.7
NeRFA (Ours) | 18.6 | 0.25 | 20.2 | 19.1 | 16.3 | 19.7 | 16.1 | 14.9 | 18.7 | 0.26 | 19.7 | 18.6 | 15.3 | 19.4 | 15.2 | 14.9

#### CO3D.

The Common Objects in 3D (CO3D [35]) is a recent dataset containing 1.5
million multi-view images of around 19k objects. This dataset measure the
ability of novel view synthesis in the wild. The $18,619$ videos are divided
following official release [35] into four splits: train-known, train-unseen,
test-known, and test-unseen. The difference between train-unseen and test-
unseen is that the train-unseen contains unseen frames in seen videos while
test-unseen contains only unseen videos. We focus on the setting of category-
centric view synthesis. Under this setting, a single model is trained on the
train-known split and is evaluated on 1000 testing samples of train-unseen and
test-unseen. Evaluation uses $10$ object categories. We do not use the auto-
decoders [35] in all baselines for fairness.

###  4.2 Experimental details

We generate $N_{p}=128$ rays in a batch for each training iteration, and we
sample $N_{v}=64$ ray points along each ray. The optimizer is Adam [14] with
an initial learning rate of $5\times 10^{-4}$ and an exponential decay of
$5\times 10^{-5}$. The hidden dimension of the attention functions is $d=64$,
and the head number is $H=8$. On the CO3D [35] dataset, we use $L=3$ layers of
transformers, and on the other datasets, we only use $L=1$. The
hyperparameters are the same for all compared approaches. We use three metrics
to evaluate NeRFA and the baselines. They are peak signal-to-noise ratio
(PSNR), structural similarity index measure (SSIM), and LPIPS [56].

###  4.3 Baselines

#### NV [20], LLFF [22] and SRN [40].

These baselines are proposed earlier than NeRF [23]. The neural volumes (NV)
generate novel views by predicting a discretized voxel grid representing
RGB$\alpha$ values with $128^{3}$ samples. LLFF [22] uses a 3D neural network
to predict a multi-plane image (MPI) for each training view. LLFF generates
novel views via blending and alpha compositing. SRN [40] is a follow-up work
of DeepVoxels [39], representing a scene via an MLP mapping from each
coordinate to a feature vector. A recurrent neural network is built upon the
feature vector to predict the step size of the ray. The color vector is
decoded from the final feature vector. On the CO3D dataset, we use the WCE
[12] technique to augment these methods following [35].

#### NeRF [23].

NeRF uses an eight-layer MLP with ReLU activations (the hidden dimension is
$256$) to learn the volume density. Another FC layer ($128$ channels with a
ReLU activation) is added to generate the view-dependent color vector. NeRF
adopts the standard volumetric rendering [21] to render views based on the MLP
predictions. Besides, NeRF introduces a hierarchical volume sampling procedure
to avoid sampling occluded regions or free space.

#### NerFormer [35].

The NerFormer replaces the MLP in NeRF with 3D transformer modules and pooling
heads. NerFormer consumes multi-source inputs, uses a series of 3D
transformers to learn the ray interaction and uses pooling layers to aggregate
per-pixel information. It keeps other components of NeRF, such as the
volumetric rendering. In contrast, NeRFA has different input and is of a
different model architecture.

#### NeRF-WCE [12]

NeRF-WCE augments the NeRF model via warp-conditioned embedding (WCE), which
is formed by sampling a tensor of the source image feature. WCE helps capture
invariance of rigid scene misalignment. When multiple source views are given,
the NeRF-WCE adopts the concatenation of the mean and standard deviation of
the view-dependent WCEs as the aggregated WCE.

Besides these baselines, we also compare with other reconstruction-based
approaches, P3DMesh [33] and IPC-WCE [51], on the CO3D dataset.

Table 3: Ablation study on the validation split of three datasets (DeepVoxels
[39], Blender [23] and LLFF [22]). VaniA, FM, RT, and PT represent the vanilla
attention, feature modulation, ray transformer and pixel transformer
correspondingly.

| DeepVoxels [39] | Blender [23] | LLFF [22]
---|---|---|---
| PSNR$\uparrow$ | SSIM$\uparrow$ | LPIPS$\downarrow$ | PSNR$\uparrow$ | SSIM$\uparrow$ | LPIPS$\downarrow$ | PSNR$\uparrow$ | SSIM$\uparrow$ | LPIPS$\downarrow$
NeRFA | 40.30 | 0.991 | 0.024 | 33.17 | 0.941 | 0.074 | 28.15 | 0.833 | 0.192
VaniA | 23.42 | 0.783 | 0.213 | 22.35 | 0.761 | 0.224 | 21.49 | 0.731 | 0.235
w/o FM | 24.51 | 0.791 | 0.194 | 25.14 | 0.793 | 0.195 | 25.42 | 0.801 | 0.208
w/o RT | 31.28 | 0.941 | 0.075 | 28.04 | 0.811 | 0.193 | 26.38 | 0.803 | 0.191
w/o PT | 37.14 | 0.973 | 0.041 | 32.42 | 0.938 | 0.076 | 28.11 | 0.818 | 0.193

Figure 4: Visualization of the ablated models. We denote the vanilla
attention, feature modulation, ray transformer and pixel transformer as VaniA,
FM, RT, and PT correspondingly.

###  4.4 Comparison

#### Single-scene view synthesis.

We first present the results under the single-scene view synthesis setting.
This setting is the same as that of NeRF [23], where one model is trained for
each scene. Results are presented in Table 1 and Figure 3. Under this setting,
SRN [40] does not learn enough details and generates blurry results. Although
NerFormer [35] replaces the MLP in NeRF [23] with a group of transformers, the
performance is not boosted. In contrast, our NeRFA model successfully improves
NeRF by generating sharper details (e.g., the edges of the fern and text of
the horns).

#### Category-centric novel-view synthesis.

We report results under the category-centric novel-view synthesis setting. We
also test the models with varying numbers of views and difficulty of target
views. The numbers are in Table 2, and some visualizations are provided in the
supplementary. The NeRFA model outperforms NerFormer [35], NeRF-WCE [12] and
other baselines on the test-unseen split. NeRFA is also robust under different
numbers of source views and different target-view difficulties. The mesh-based
method P3D-mesh [33] yields competitive results on challenging scenarios of
the train-unseen split. However, it works poorly on unseen videos in the test-
unseen split.

###  4.5 Ablation study

In this subsection, we investigate the effects of each NeRFA component.
Besides the vanilla attention (VaniA), which is presented in Sec. 3.2, we
compare the following ablated variants of NeRFA. The unmentioned model
components are the same as the NeRFA model.

#### NeRFA w/o feature modulation (w/o FM).

We remove the feature modulation step from the NeRFA model. An average pooling
layer replaces the feature modulation layer. The pooled features
$\mathbf{F_{pool}}\in\mathbb{R}^{N_{p}\times d}$ is then sent to the pixel
transformer to get the color features $\mathbf{F_{color}}$.

#### NeRFA w/o ray transformer (w/o RT).

We remove the ray transformer, and the positional embedding is directly sent
to the feature modulation layer. In other words, we set the positional
embeddings as the pixel-ray features
$\mathbf{F^{\prime}_{pv}}=\mathbf{F_{emb}}$.

#### NeRFA w/o pixel transformer (w/o PT).

In this variant, the pixel transformer is removed from the NeRFA model, and
the modulated features are directly sent to the final FC layer. In other
words, the color features are set to the modulated features
$\mathbf{F_{color}}=\mathbf{F_{mod}}$.

All the ablated models share the same hyper-parameters and training
procedures. We present the results on the validation split of three datasets
(DeepVoxels [39], Blender [23] and LLFF [22]) in Table 3. The visualizations
are presented in Figure 4. We observe that the feature modulation effectively
helps learn the high-frequency components of the scene. The ray transformer
reduces the blurry and distorted renderings. We also observe that the pixel
transformer enhances sharp details.

##  5 Conclusions and future work

We have presented a simple end-to-end framework for view synthesis. Because
vanilla attention could not model the rendering procedure well, we proposed
the NeRF attention (NeRFA) model, which considered the NeRF rendering process
a soft feature generation procedure. The NeRFA model used ray and pixel
transformers to learn connections between rays and pixels. We empirically
verified that NeRFA significantly improved vanilla attention and helped learn
high-frequency components. NeRFA also demonstrated superior performance over
strong baselines such as NerFormer [35] and NeRF-WCE [12]. Our work suggests
multiple future work directions, such as exploring better transformer models
[19] under the seq2seq formulation and learning representations across
categories [36]. Learning to reconstruct surfaces [49] will also be our
target.

## References

  * [1] Bommasani, R., Hudson, D.A., Adeli, E., Altman, R., Arora, S., von Arx, S., Bernstein, M.S., Bohg, J., Bosselut, A., Brunskill, E., et al.: On the opportunities and risks of foundation models. arXiv preprint arXiv:2108.07258 (2021)
  * [2] Carion, N., Massa, F., Synnaeve, G., Usunier, N., Kirillov, A., Zagoruyko, S.: End-to-end object detection with transformers. In: European Conference on Computer Vision. pp. 213–229. Springer (2020)
  * [3] Cheng, S., Min, S., Hwann-Tzong, C.: Direct voxel grid optimization: Super-fast convergence for radiance fields reconstruction. arXiv preprint arXiv:2111.11215 (2021)
  * [4] Devlin, J., Chang, M.W., Lee, K., Toutanova, K.: Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:1810.04805 (2018)
  * [5] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.: An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint arXiv:2010.11929 (2020)
  * [6] Esser, P., Rombach, R., Ommer, B.: Taming transformers for high-resolution image synthesis. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 12873–12883 (2021)
  * [7] Flynn, J., Broxton, M., Debevec, P., DuVall, M., Fyffe, G., Overbeck, R., Snavely, N., Tucker, R.: Deepview: View synthesis with learned gradient descent. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 2367–2376 (2019)
  * [8] Gafni, G., Thies, J., Zollhofer, M., Nießner, M.: Dynamic neural radiance fields for monocular 4d facial avatar reconstruction. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 8649–8658 (2021)
  * [9] Gao, P., Lu, J., Li, H., Mottaghi, R., Kembhavi, A.: Container: Context aggregation network. CoRR abs/2106.01401 (2021), https://arxiv.org/abs/2106.01401
  * [10] Garbin, S.J., Kowalski, M., Johnson, M., Shotton, J., Valentin, J.: Fastnerf: High-fidelity neural rendering at 200fps. https://arxiv.org/abs/2103.10380 (2021)
  * [11] Gortler, S.J., Grzeszczuk, R., Szeliski, R., Cohen, M.F.: The lumigraph. In: Proceedings of the 23rd annual conference on Computer graphics and interactive techniques. pp. 43–54 (1996)
  * [12] Henzler, P., Reizenstein, J., Labatut, P., Shapovalov, R., Ritschel, T., Vedaldi, A., Novotny, D.: Unsupervised learning of 3d object categories from videos in the wild. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). pp. 4700–4709 (June 2021)
  * [13] Kar, A., Häne, C., Malik, J.: Learning a multi-view stereo machine. In: Guyon, I., Luxburg, U.V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., Garnett, R. (eds.) Advances in Neural Information Processing Systems. vol. 30. Curran Associates, Inc. (2017), https://proceedings.neurips.cc/paper/2017/file/9c838d2e45b2ad1094d42f4ef36764f6-Paper.pdf
  * [14] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [15] Kitaev, N., Kaiser, Ł., Levskaya, A.: Reformer: The efficient transformer. arXiv preprint arXiv:2001.04451 (2020)
  * [16] Kulhánek, J., Derner, E., Sattler, T., Babuška, R.: Viewformer: Nerf-free neural rendering from few images using transformers. arXiv preprint arXiv:2203.10157 (2022)
  * [17] Lindell, D., Martel, J., Wetzstein, G.: AutoInt: Automatic integration for fast neural volume rendering. https://arxiv.org/abs/2012.01714 (2020)
  * [18] Liu, L., Gu, J., Zaw Lin, K., Chua, T.S., Theobalt, C.: Neural sparse voxel fields. Advances in Neural Information Processing Systems 33, 15651–15663 (2020)
  * [19] Liu, Z., Lin, Y., Cao, Y., Hu, H., Wei, Y., Zhang, Z., Lin, S., Guo, B.: Swin transformer: Hierarchical vision transformer using shifted windows. arXiv preprint arXiv:2103.14030 (2021)
  * [20] Lombardi, S., Simon, T., Saragih, J., Schwartz, G., Lehrmann, A., Sheikh, Y.: Neural volumes: Learning dynamic renderable volumes from images. arXiv preprint arXiv:1906.07751 (2019)
  * [21] Max, N.: Optical models for direct volume rendering. IEEE Transactions on Visualization and Computer Graphics 1(2), 99–108 (1995)
  * [22] Mildenhall, B., Srinivasan, P.P., Ortiz-Cayon, R., Kalantari, N.K., Ramamoorthi, R., Ng, R., Kar, A.: Local light field fusion: Practical view synthesis with prescriptive sampling guidelines. ACM Transactions on Graphics (TOG) 38(4), 1–14 (2019)
  * [23] Mildenhall, B., Srinivasan, P.P., Tancik, M., Barron, J.T., Ramamoorthi, R., Ng, R.: Nerf: Representing scenes as neural radiance fields for view synthesis. In: European conference on computer vision. pp. 405–421. Springer (2020)
  * [24] Mohammad, S.M.: Sentiment analysis: Detecting valence, emotions, and other affectual states from text. In: Emotion measurement, pp. 201–237. Elsevier (2016)
  * [25] Nguyen-Phuoc, T.H., Li, C., Balaban, S., Yang, Y.: Rendernet: A deep convolutional network for differentiable rendering from 3d shapes. Advances in Neural Information Processing Systems 31 (2018)
  * [26] Niemeyer, M., Geiger, A.: Giraffe: Representing scenes as compositional generative neural feature fields. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 11453–11464 (2021)
  * [27] Park, K., Sinha, U., Barron, J.T., Bouaziz, S., Goldman, D.B., Seitz, S.M., Martin-Brualla, R.: Nerfies: Deformable neural radiance fields. ICCV (2021)
  * [28] Park, N., Kim, S.: How do vision transformers work? arXiv preprint arXiv:2202.06709 (2022)
  * [29] Penner, E., Zhang, L.: Soft 3d reconstruction for view synthesis. ACM Transactions on Graphics (TOG) 36(6), 1–11 (2017)
  * [30] Pumarola, A., Corona, E., Pons-Moll, G., Moreno-Noguer, F.: D-nerf: Neural radiance fields for dynamic scenes. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 10318–10327 (2021)
  * [31] Radford, A., Wu, J., Child, R., Luan, D., Amodei, D., Sutskever, I.: Language models are unsupervised multitask learners. OpenAI blog 1(8), 9 (2019)
  * [32] Raj, A., Zollhoefer, M., Simon, T., Saragih, J., Saito, S., Hays, J., Lombardi, S.: Pva: Pixel-aligned volumetric avatars. arXiv preprint arXiv:2101.02697 (2021)
  * [33] Ravi, N., Reizenstein, J., Novotny, D., Gordon, T., Lo, W.Y., Johnson, J., Gkioxari, G.: Accelerating 3d deep learning with pytorch3d. arXiv preprint arXiv:2007.08501 (2020)
  * [34] Reiser, C., Peng, S., Liao, Y., Geiger, A.: Kilonerf: Speeding up neural radiance fields with thousands of tiny mlps (2021)
  * [35] Reizenstein, J., Shapovalov, R., Henzler, P., Sbordone, L., Labatut, P., Novotny, D.: Common objects in 3d: Large-scale learning and evaluation of real-life 3d category reconstruction. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 10901–10911 (2021)
  * [36] Rombach, R., Esser, P., Ommer, B.: Geometry-free view synthesis: Transformers and no 3d priors. In: Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV). pp. 14356–14366 (October 2021)
  * [37] Schwarz, K., Liao, Y., Niemeyer, M., Geiger, A.: Graf: Generative radiance fields for 3d-aware image synthesis. Advances in Neural Information Processing Systems 33, 20154–20166 (2020)
  * [38] Shysheya, A., Zakharov, E., Aliev, K.A., Bashirov, R., Burkov, E., Iskakov, K., Ivakhnenko, A., Malkov, Y., Pasechnik, I., Ulyanov, D., et al.: Textured neural avatars. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 2387–2397 (2019)
  * [39] Sitzmann, V., Thies, J., Heide, F., Niessner, M., Wetzstein, G., Zollhofer, M.: Deepvoxels: Learning persistent 3d feature embeddings. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) (June 2019)
  * [40] Sitzmann, V., Zollhöfer, M., Wetzstein, G.: Scene representation networks: Continuous 3d-structure-aware neural scene representations. Advances in Neural Information Processing Systems 32 (2019)
  * [41] Srinivasan, P.P., Tucker, R., Barron, J.T., Ramamoorthi, R., Ng, R., Snavely, N.: Pushing the boundaries of view extrapolation with multiplane images. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) (June 2019)
  * [42] Sun, C., Myers, A., Vondrick, C., Murphy, K., Schmid, C.: Videobert: A joint model for video and language representation learning. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 7464–7473 (2019)
  * [43] Tan, H., Bansal, M.: Lxmert: Learning cross-modality encoder representations from transformers. In: Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP). pp. 5103–5114 (2019)
  * [44] Tancik, M., Srinivasan, P., Mildenhall, B., Fridovich-Keil, S., Raghavan, N., Singhal, U., Ramamoorthi, R., Barron, J., Ng, R.: Fourier features let networks learn high frequency functions in low dimensional domains. Advances in Neural Information Processing Systems 33, 7537–7547 (2020)
  * [45] Thies, J., Zollhöfer, M., Nießner, M.: Deferred neural rendering: Image synthesis using neural textures. ACM Transactions on Graphics (TOG) 38(4), 1–12 (2019)
  * [46] Tulsiani, S., Zhou, T., Efros, A.A., Malik, J.: Multi-view supervision for single-view reconstruction via differentiable ray consistency. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 2626–2634 (2017)
  * [47] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, L., Polosukhin, I.: Attention is all you need. arXiv preprint arXiv:1706.03762 (2017)
  * [48] Waechter, M., Moehrle, N., Goesele, M.: Let there be color! large-scale texturing of 3d reconstructions. In: European conference on computer vision. pp. 836–850. Springer (2014)
  * [49] Wang, P., Liu, L., Liu, Y., Theobalt, C., Komura, T., Wang, W.: Neus: Learning neural implicit surfaces by volume rendering for multi-view reconstruction. NeurIPS (2021)
  * [50] Wang, Q., Wang, Z., Genova, K., Srinivasan, P.P., Zhou, H., Barron, J.T., Martin-Brualla, R., Snavely, N., Funkhouser, T.: Ibrnet: Learning multi-view image-based rendering. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 4690–4699 (2021)
  * [51] Wiles, O., Gkioxari, G., Szeliski, R., Johnson, J.: Synsin: End-to-end view synthesis from a single image. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 7467–7477 (2020)
  * [52] Wood, D.N., Azuma, D.I., Aldinger, K., Curless, B., Duchamp, T., Salesin, D.H., Stuetzle, W.: Surface light fields for 3d photography. In: Proceedings of the 27th annual conference on Computer graphics and interactive techniques. pp. 287–296 (2000)
  * [53] Yang, Z., Dai, Z., Yang, Y., Carbonell, J., Salakhutdinov, R., Le, Q.V.: Xlnet: Generalized autoregressive pretraining for language understanding. arXiv preprint arXiv:1906.08237 (2019)
  * [54] Yariv, L., Gu, J., Kasten, Y., Lipman, Y.: Volume rendering of neural implicit surfaces. In: NeurIPS (2021)
  * [55] Zhang, Q., Yang, Y.: Rest: An efficient transformer for visual recognition. CoRR abs/2105.13677 (2021), https://arxiv.org/abs/2105.13677
  * [56] Zhang, R., Isola, P., Efros, A.A., Shechtman, E., Wang, O.: The unreasonable effectiveness of deep features as a perceptual metric. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 586–595 (2018)
  * [57] Zhao, Z., Gan, C., Wu, J., Guo, X., Tenenbaum, J.: Augmenting policy learning with routines discovered from a demonstration. arXiv preprint arXiv:2012.12469 (2020)
  * [58] Zhao, Z., Peng, G., Wang, H., Fang, H.S., Li, C., Lu, C.: Estimating 6d pose from localizing designated surface keypoints. arXiv preprint arXiv:1812.01387 (2018)
  * [59] Zhao, Z., Samel, K., Chen, B., song, l.: Proto: Program-guided transformer for program-guided tasks. In: Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., Vaughan, J.W. (eds.) Advances in Neural Information Processing Systems. vol. 34, pp. 17021–17036. Curran Associates, Inc. (2021), https://proceedings.neurips.cc/paper/2021/file/8d34201a5b85900908db6cae92723617-Paper.pdf
  * [60] Zhao, Z., Wu, Z., Zhuang, Y., Li, B., Jia, J.: Tracking objects as pixel-wise distributions (2022). https://doi.org/10.48550/ARXIV.2207.05518, https://arxiv.org/abs/2207.05518
  * [61] Zhou, T., Tulsiani, S., Sun, W., Malik, J., Efros, A.A.: View synthesis by appearance flow. In: European conference on computer vision. pp. 286–301. Springer (2016)
  * [62] Zhu, J.Y., Zhang, Z., Zhang, C., Wu, J., Torralba, A., Tenenbaum, J., Freeman, B.: Visual object networks: Image generation with disentangled 3d representations. Advances in neural information processing systems 31 (2018)
  * [63] Zhu, X., Su, W., Lu, L., Li, B., Wang, X., Dai, J.: Deformable detr: Deformable transformers for end-to-end object detection. arXiv preprint arXiv:2010.04159 (2020)