# NerfAcc: Efficient Sampling Accelerates NeRFs

Ruilong Li
UC Berkeley
<EMAIL_ADDRESS>Hang Gao
UC Berkeley
<EMAIL_ADDRESS>Matthew Tancik
UC Berkeley
<EMAIL_ADDRESS>Angjoo Kanazawa
UC Berkeley
<EMAIL_ADDRESS>
###### Abstract

Optimizing and rendering Neural Radiance Fields is computationally expensive
due to the vast number of samples required by volume rendering. Recent works
have included alternative sampling approaches to help accelerate their
methods, however, they are often not the focus of the work. In this paper, we
investigate and compare multiple sampling approaches and demonstrate that
improved sampling is generally applicable across NeRF variants under an
unified concept of transmittance estimator. To facilitate future experiments,
we develop NerfAcc, a Python toolbox that provides flexible APIs for
incorporating advanced sampling methods into NeRF related methods. We
demonstrate its flexibility by showing that it can reduce the training time of
several recent NeRF methods by 1.5× to 20× with minimal modifications to the
existing codebase. Additionally, highly customized NeRFs, such as Instant-NGP,
can be implemented in native PyTorch using NerfAcc. Our code are open-sourced
at https://www.nerfacc.com.

Figure 1: NerfAcc Toolbox. Our proposed toolbox, _NerfAcc_ , integrates
advanced efficient sampling techniques that lead to significant speedups in
training various recent NeRF papers with minimal modifications to existing
codebases.

##  1 Introduction

Neural volume rendering has revolutionized the inverse rendering problem, with
the Neural Radiance Field (NeRF) [16] being a key innovation. The continuous
radiance field representation allows for rendering novel views of a scene from
any camera position. However, the optimization of a NeRF can be
computationally expensive due to the neural representation of radiance field
and the large number of samples required by volume rendering. These challenges
have limited practical applications of NeRF-based optimization and rendering.

Several recent works have successfully reduced the computational cost of
neural volume rendering by proposing more efficient radiance field
representations [17, 35, 36, 4, 24, 7]. While there are differences in the
specific radiance field representations and their applications, most of these
methods share a similar volume rendering pipeline which involves creating
samples along the ray and accumulating them through alpha-composition.

However, compared to the considerable efforts focused on developing efficient
radiance field representations, there has been limited attention devoted to
reducing the computational cost of neural volume rendering through efficient
sampling. While a few recent works have included alternative sampling
approaches to accelerate their methods [17, 35, 24, 2], these methods are
often not the main focus of the paper. Moreover, implementing advanced
sampling approaches typically requires non-trivial efforts. For example,
Instant-NGP [17] and Plenoxels [35] both employ highly customized CUDA
implementations to achieve spatial skipping during ray marching, which are
tightly coupled with their respective radiance field implementations.
Consequently, it can be challenging for researchers to benefit from these
advanced sampling approaches in their own research.

In this paper, we investigate and compare several advanced sampling approaches
from the literature and provide mathematical proofs demonstrating that they
can all be viewed in a unified way of creating an _estimation of
transmittance_ for importance sampling. Our analysis shows that by
understanding the spectrum of sampling through transmittance estimator, novel
sampling strategies can be explored. Based on this, we decouple the sampling
procedure from the neural volumetric rendering pipeline and demonstrate that
improved sampling is generally applicable across different variants of NeRF.
Furthermore, to facilitate future experiments, we propose NerfAcc, a plug-and-
play toolbox that provides a flexible Python API for integrating advanced
sampling approaches into NeRF-related methods, ready for researchers to
incorporate into their own codebases. We demonstrate that with less than $100$
lines of code change using NerfAcc, various NeRF methods [16, 4, 5, 7, 14, 21,
17] can enjoy 1.5$\times$ to 20$\times$ training speedup with better
performance. Notably, using the NerfAcc library, one can train an Instant-NGP
[17] model with pure Python code in the same amount of time, and achieve
slightly better performance ($+0.2$dB) than the official pure CUDA
implementation.

Our paper presents a unique contribution to the community. Unlike other papers
that introduce novel algorithms, our work sheds light on the intricacies of
various sampling approaches, which are often overlooked despite their
significance. As far as we know, this work is the first paper that dives deep
into this crucial aspect in the context of neural radiance field, and offer a
novel, unified concept that allows researchers to view existing sampling
algorithms in a fresh perspective. In addition, we provide a toolbox that
facilitates research and development in this area. We believe translation of
mathematical ideas into efficient, easy-to-modify implementation is
fundamental to the research development. Overall, we hope that our proposed
concept, along with the toolbox, can inspire new research ideas and accelerate
progress in this field.

##  2 Related Works

#### NeRF Codebases.

The recent explosion of NeRF-related research has led to numerous papers, many
of which have released their own codebases [1, 2, 4, 7, 13, 16, 19, 37, 9].
These codebases address various tasks related to NeRF, including surface
reconstruction [18, 32, 34], radiance field representation [1, 2, 17, 24, 35],
dynamic modeling [5, 7, 13, 21], and camera optimization [14, 33]. However,
each codebase is tailored to a specific task and supports only a single
approach. While most of these methods share the same volume rendering pipeline
of accumulating samples along the ray, transferring the implementation from
one codebase to another requires non-trivial efforts. In this work, we address
this common sampling problem by introducing NerfAcc, a plug-and-play toolbox
that supports all the aforementioned tasks and methods, making it easy to
integrate into any existing codebase.

#### NeRF Frameworks.

Recently, several projects have been introduced to integrate different NeRF
variants into a single framework, such as NeRF-Factory [11], Nerfstudio [29],
and Kaolin-Wisp [28]. These frameworks have made significant progress in
facilitating future NeRF-related research. NeRF-Factory offers a collection of
NeRF variants [16, 39, 24, 35, 1, 2, 30] with original implementations and
focuses on comprehensive benchmarking. Nerfstudio consolidates critical
techniques introduced in existing literature [16, 1, 2, 17, 15] and provides
modular components for the community to easily build on. Kaolin-Wisp builds
upon the Kaolin [8] framework and implements a set of voxel-based NeRF papers
[27, 17, 26]. However, these frameworks are designed to encourage researchers
to develop within the framework, and do not benefit users working on their own
codebases. Moreover, due to the high activity in the NeRF-related research, it
is almost impossible to keep track of the latest developments and integrate
advanced techniques into a single framework. Therefore, NerfAcc is designed as
a standalone library that can be plugged into any codebase. It supports a wide
range of NeRF related methods and can be easily maintained as new methods
emerge.

##  3 Importance Sampling via Transmittance

Figure 2: Illustration of Sampling via Transmittance Estimator. Although
spatial skipping approaches (e.g., Occupancy Grid in Instant-NGP [17]) and PDF
approaches (e.g., Proposal Network in Mip-NeRF 360 [2]) appear distinct from
each other, they can both be viewed as constructing a transmittance estimator
from which samples can be created via importance sampling.

Several advanced sampling approaches exist in the literature. For instance,
Plenoxels [35] uses a sparse grid, Instant-NGP [17] uses an occupancy grid,
NeRF [16] employs a coarse-to-fine strategy, Mip-NeRF 360 [2] proposes
proposal networks. However, creating a toolbox that can support all of these
approaches is not an easy task since they operate in completely different
ways. In this section, we demonstrate that these approaches can all be viewed
in a unified way of importance sampling. We also provide a mathematical proof
that _transmittance is all you need for importance sampling_. Thus, each
method essentially has its own way of creating an estimation of the
transmittance along the ray, which we refer to as the _transmittance
estimator_. This observation enables us to unify different types of sampling
approaches under the same framework and leads to our NerfAcc toolbox.

###  3.1 Formulation

Efficient sampling is a well-explored problem in Graphics [6], wherein the
emphasis is on identifying regions that make the most significant contribution
to the final rendering. This objective is generally accomplished through
importance sampling, which aims to distribute samples based on the probability
density function (PDF), denoted as $p(t)$, between the range of
$[t_{n},t_{f}]$. By computing the cumulative distribution function (CDF)
through integration, i.e., $F(t)=\int_{t_{n}}^{t}p(v)\,dv$, samples are
generated using the inverse transform sampling method:

$t=F^{-1}(u)\quad\text{where}\quad u\sim\mathcal{U}[0,1].$ (1)

In volumetric rendering, the contribution of each sample to the final
rendering is expressed by the accumulation weights $T(t)\sigma(t)$:

$\begin{gathered}C(\mathbf{r})=\int_{t_{n}}^{t_{f}}T(t)\,\sigma(t)\,c(t)\,dt\\\
\quad\text{where}\quad
T(t)=\exp\left(-\int_{t_{n}}^{t}\,\sigma(s)\,ds\right).\end{gathered}$ (2)

Hence, the PDF for volumetric rendering is $p(t)=T(t)\sigma(t)$ and the CDF
$F(t)=\int_{t_{n}}^{t}p(v)\,dv$ can be derived as a function of transmittance
$T(t)$:

$\displaystyle F(t)$ $\displaystyle=\int_{t_{n}}^{t}T(v)\sigma(v)\,dv$ (3)
$\displaystyle=\int_{t_{n}}^{t}\frac{d}{dv}\left[-\exp\left(-\int_{t_{n}}^{v}\sigma(s)ds\right)\right]dv$
$\displaystyle=1-\exp\left(-\int_{t_{n}}^{t}\sigma(s)ds\right)$
$\displaystyle=1-T(t).$

Sampling Methods |  PDF Estimator $p(t_{i})$ |  Transmittance Estimator $T(t_{i})$ |  Update $\mathcal{F}$ |  Instantiation
---|---|---|---|---
Uniform |  Constant |  $1-\frac{t-t_{n}}{t_{f}-t_{n}}$ |  - |  Coarse MLP NeRF [16]
Spatial Skipping |  $\frac{\mathbb{1}\big{[}\sigma(t_{i})>\tau\big{]}}{\sum_{j=1}^{n}\mathbb{1}\big{[}\sigma(t_{i})>\tau\big{]}}$ |  $1-\frac{\sum_{j=1}^{i-1}\mathbb{1}\big{[}\sigma(t_{i})>\tau\big{]}}{\sum_{j=1}^{n}\mathbb{1}\big{[}\sigma(t_{i})>\tau\big{]}}$ |  EMA |  Instant-NGP [17]
|  |  |  SGD |  Plenoxels [35]
PDF Approaches |  $\sigma(t_{i})\exp(-\sigma(t_{i})\,dt)$ |  $\exp(-\sum_{j=1}^{i-1}\sigma(t_{i})\,dt)$ |  SGD |  Fine MLP NeRF [16]
|  |  |  SGD |  Mip-NeRF 360 [2]
Table 1:  Mathematical Formulations of Different Sampling Approaches. We
outlines the significant mathematical distinctions among each sampling
approach, under the perspective of transmittance estimator. See Section 3.1
for notations.

Therefore, inverse sampling the CDF $F(t)$ is equivalent to inverse sampling
the transmittance $T(t)$. Thus a transmittance estimator is sufficient to
determine the optimal samples. Intuitively, this suggests to put more samples
around regions where the transmittance changes rapidly – and that is exactly
what happens when a ray hit a surface. Implementation wise, this observation
enables us to compute the CDF directly using $1-T(t)$, instead of the
computationally expensive integral $\int_{t_{n}}^{t}T(v)\sigma(v)\,dv$, which
is the standard implementation adopted by many popular codebases [16, 1, 2,
29, 7].

While advanced transmittance estimator techniques, such as delta tracking [31]
with Monte Carlo sampling, are utilized in production-level volumetric
rendering in Graphics [6], NeRFs operate in a distinct setting in which the
scene geometry is not predefined but optimized on the fly. During NeRF
optimization, the radiance field changes between iterations, necessitating the
dynamic update of the transmittance estimator at each step $k$:

$\mathcal{F}:T(t)^{k-1}\mapsto T(t)^{k}.$ (4)

This introduces additional challenges to efficient sampling because it becomes
more difficult to accurately estimate the transmittance from a radiance field
that is constantly changing. Current approaches employ either exponential
moving average (EMA) or stochastic gradient descent (SGD) as the update
function $\mathcal{F}$. However, we note that there may be other update
functions that could be explored.

With these concepts in mind, let us now examine some of the existing
approaches towards efficient sampling.

#### Uniform.

If a transmittance estimator is not available, the only assumption we can make
is that every point along the ray contributes equally to the final rendering.
Mathematically, this assumption translates to a constant PDF and a linearly
decaying transmittance $T(t)=1-(t-t_{n})/(t_{f}-t_{n})$. In this case, the
sampling process is equivalent to uniformly sampling along the ray, i.e.,
$t_{i}=t_{n}+(t_{f}-t_{n})\cdot u_{i}$. It is worth noting that every NeRF
model that uses uniform sampling inherently assumes this linear transmittance
decay, such as the coarse level in vanilla NeRF [16]. See Fig. 2(a) for the
illustration.

#### Spatial Skipping.

A more sophisticated approach to improve uniform sampling is to identify empty
regions and skip them during sampling, as proposed in Instant-NGP’s Occupancy
Grid [17] and Plenoctrees’ Sparse Grid [36]. This technique binarizes the
density along the ray with a conservative threshold $\tau$:
$\hat{\sigma}(t_{i})=\mathbb{1}\big{[}\sigma(t_{i})>\tau\big{]}$.
Consequently, the piece-wise constant PDF can be expressed as
$p(t_{i})=\hat{\sigma}(t_{i})/\sum_{j=1}^{n}\hat{\sigma}(t_{j})$, and the
piece-wise linear transmittance estimator is
$T(t_{i})=1-\sum_{j=1}^{i-1}\hat{\sigma}(t_{j})/\sum_{j=1}^{n}\hat{\sigma}(t_{j})$.
To update this estimator during optimization, Instant-NGP [17] directly
updates the cached density with exponential moving average (EMA) over
iteration $k$:
$\sigma(t_{i})^{k}=\gamma\cdot\sigma(t_{i})^{k-1}+(1-\gamma)\cdot\sigma(t_{i})^{k}$.
Meanwhile, Plenoxels [35] updates the density via the gradient descent through
the rendering loss. See Fig. 2(b) for an illustration.

#### PDF Approaches.

Another type of approach is to directly estimate the PDF along the ray with
discrete samples. In vanilla NeRF [16], the coarse MLP is trained using
volumetric rendering loss to output a set of densities ${\sigma(t_{i})}$. This
allows for the creation of a piece-wise constant PDF:
$p(t_{i})=\sigma(t_{i})\exp(-\sigma(t_{i})\,dt)$, and a piece-wise linear
transmittance estimator $T(t_{i})=\exp(-\sum_{j=1}^{i-1}\sigma(t_{i})\,dt)$.
This approach was further improved in Mip-NeRF 360 [16] with a PDF matching
loss, which allows for the use of a much smaller MLP in the coarse level,
namely Proposal Network, to speedup the PDF construction. In both cases, the
transmittance estimator is updated through gradient descent. See Fig. 2(c) for
an illustration.

Table 1 provides a mathematical summary and comparison of these approaches.
Additionally, we present an illustration in Figure 2 to provide an intuitive
comparison of these approaches from the perspective of PDF (second row) and
transmittance (third row), as well as how the samples can be created from the
transmittance estimator via importance sampling (last row). This visualization
also reveals some pros and cons for each approach, which we will discuss in in
Section. 3.3.

###  3.2 Design Spaces

#### Choice of Representations.

The transmittance estimator can use either an explicit voxel [17, 35, 4], an
MLP [1, 2, 16], or a hybrid representation [29]. Depending on whether the
estimator is explicit or not, it can be updated with either rule-based EMA
[17, 4] or gradient descent with some supervision [35, 2, 16, 29]. Generally,
voxel-based estimators are faster than implicit (e.g., MLP-based) estimators
but suffer more from aliasing issues. It is worth noting that the
transmittance estimator’s representation can significantly benefit from
advances in radiance field representation. For example, the Nerfacto model
[29] uses the most recent hybird-representation HashEncoding [17] for both the
radiance field and the sampling module, achieving the best quality-speed
tradeoff in in-the-wild settings.

#### Handling Unbounded Scenes.

So far, we have only discussed sampling within a bounded area $[t_{n},t_{f}]$.
For unbounded scenes, it is impossible to densely sample along the ray.
Similar to the mipmaps used in graphics rendering, a general solution is to
sample more coarsely as the ray goes further, as objects farther away appear
in fewer pixels on the image plane. This can be achieved by creating a
bijective mapping function $\Phi:s\in[s_{n},s_{f}]\mapsto
t\in[t_{n},+\infty]$, and performing sampling in the $s$-space instead of the
$t$-space. Several papers [2, 39, 23] that work on unbounded scenes have
introduced different mapping functions $\Phi$, to which we refer our readers
for details.

###  3.3 Discussions

#### Pros and Cons.

Sampling with uniform assumption is the easiest one to implement but with
lowest efficiency in most cases. Spatial skipping is a more efficient
technique since most of the 3D space is empty, but it still samples uniformly
within occupied but occluded areas that contribute little to the final
rendering (e.g., the last sample in Figure 2(b)). PDF-based estimators
generally provide more accurate transmittance estimation, enabling samples to
concentrate more on high-contribution areas (e.g., surfaces) and to be more
spread out in both empty and occluded regions. However, this also means that
samples are always spread out throughout the entire space without any
skipping, as shown in Figure 2(c). Moreover, the current approaches all
introduce aliasing effects to the volumetric rendering due to either (1) the
piece-wise linearity assumption for estimating transmittance along the ray [1,
2, 16, 35, 17] as illustrated in Figure 2, or (2) the underlying voxel
representation for the transmittance estimator [35, 17] as discussed in
Section 3.2. A recent work, Zip-NeRF [3], addresses the aliasing problem tied
to these two exact issues (called “z-aliasing” and “xy-aliasing” in their
work), which are naturally revealed under our unified framework.

#### Implementation Difficulties.

The current implementations for efficient sampling are all highly customized
and tightly integrated with the specific radiance field proposed in each
paper. For instance, spatial skipping is implemented with customized CUDA
kernels in Instant-NGP [17] and Plenoxels [35]. Mip-NeRF 360 [2], K-planes
[7], and Nerfacto [29] have a proposal network implemented but it is closely
integrated with their repositories and can only support limited types of
radiance fields that come with the repository. However, as shown before, the
sampling process is independent of the radiance field representation, thus it
should be easily transferable across different NeRF variants. Due to the
various implementation details, it typically requires significant effort to
correctly implement an efficient sampling approach from scratch. Therefore,
having an implementation that is easily transferable from repository to
repository would be valuable in supporting future research on NeRF.

#### Insights from Unified Formulation.

Comprehending the sampling spectrum through the lens of Transmittance
Estimator paves the way for investigating novel sampling strategies. For
example, our framework reveals that the Occupancy Grid from Instant-NGP [17]
and the Proposal Network from Mip-NeRF 360 [2] are not mutually exclusive but
complementary, as both aim to estimate the transmittance along the ray.
Therefore, combining them becomes straight-forward: one can first compute the
transmittance with the occupancy grid, and then refine the estimated
transmittance with a proposal network. This enables both skipping on the empty
space and concentrating the samples onto the surface. We explore this approach
in Section 4.4 and demonstrate that it overcomes the limitation of the
proposal network approach, which always samples the entire space. Furthermore,
this formulation could potentially shed light on questions such as how to
enhance the sampling procedure with depth information or other priors, which
we encourage readers to investigate further.

##  4 NerfAcc Toolbox

Algorithm 1 : NerfAcc Rendering Pipeline. In NerfAcc, the samples are created
from a transmittance estimator, which can be updated during the NeRF training.
See Section 4.2 for details.

⬇

# nerf: a radiance field model.

# r_o: ray origins. (n_rays, 3)

# r_d: ray normalized directions. (n_rays, 3)

# e.g., Prop Net, Occ Grid.

estimator = nerfacc.TransmittanceEstimator()

def density_fn(t0, t1, r_id):

"""Query␣density␣from␣nerf."""

return nerf.density(r_o[r_id], r_d[r_id], t0, t1)

def rgb_density_fn(t0, t1, r_id):

"""Query␣rgb␣and␣density␣from␣nerf."""

return nerf.forward(r_o[r_id], r_d[r_id], t0, t1)

# Efficient sampling.

# (t0, t1, r_id): packed samples. (all_samples,)

t0, t1, r_id = nerfacc.sampling(

r_o, r_d, estimator, density_fn=density_fn

)

# Differentiable volumetric rendering.

color, opacity, depth, aux = nerfacc.rendering(

t0, t1, r_id, rgb_density_fn=rgb_density_fn

)

# Update the transmittance estimator.

estimator.update_every_n_steps(t0, t1, r_id, aux)

# nerf, r_o and r_d all receive gradients.

F.mse_loss(color, color_gt).backward()

In this paper, we present _NerfAcc_ toolbox, designed for Neural radiance
field Acceleration. It provides efficient sampling for volumetric rendering,
that is universally applicable and easily integrable for a diverse range of
radiance fields [16, 4, 17, 14, 32]. In this section, we first introduce the
design principles of this toolbox, along with critical implementation details.
To demonstrate its flexibility, we further show that it can significantly
speedup the training of various NeRF-related papers by $1.5\times$ to
$20\times$ with only minor modifications to existing codebases.

###  4.1 Design Principles

This library is designed with these goals in mind:

  * •

Plug-and-play. Our primary objective is to ease the challenges of implementing
an efficient volumetric sampling technique for the research community.
Therefore, NerfAcc is designed as a standalone library that can be easily
installed from PyPI on both Windows and Linux platforms, and seamlessly
integrated into any PyTorch codebase.

  * •

Efficiency & Flexibility. To maximize the speed of the code, we fuse the
operations into CUDA kernels as much as possible, while exposing flexible
Python APIs to the users.

  * •

Radiance Field Complexity. We target on supporting any radiance fields that
are designed for per-scene optimization, including both density-based and SDF-
based fields, for both static and dynamic scenes.

Original |  |  |  |  |
---|---|---|---|---|---
w/ NerfAcc |  |  |  |  |
|  TensoRF [4] |  Instant-NGP [17] |  T-NeRF [16] |  TiNeuVox [5] |  BARF [14]
|  |  |  |  |
|  Tanks&Temples Dataset |  360 Dataset |  D-NeRF Dataset |  HyperNeRF Dataset |  NeRF-Synthetic Dataset
Figure 3:  Qualitative Results. NerfAcc is able to significantly reduce the
training time of various NeRF-related methods across multiple datasets, while
also yielding slightly improved quality in the majority of cases. The training
time and test PSNR are shown on the left corner of each image.

Methods |  Dataset |  Speedup |  $\mathbf{T}$ (min) $\downarrow$ |  PSNR $\uparrow$ |  LPIPS $\downarrow$
---|---|---|---|---|---
TensoRF [4] |  T&T |  $\mathbf{1.5\times}$ |  18.3 |  28.13 |  0.143
\+ nerfacc (occ) |  |  |  12.6 |  28.10 |  0.150
TensoRF [4] |  Syn. |  $\mathbf{1.6\times}$ |  10.6 |  32.52 |  0.047
\+ nerfacc (occ) |  |  |  6.5 |  32.51 |  0.044
NeRF [16] |  Syn. |  $\mathbf{20\times}$ |  $>$1000 |  31.00 |  0.081
\+ nerfacc (occ)† |  |  |  52.0 |  31.55 |  0.072
Instant-NGP [17] |  Syn. |  $\mathbf{1.0\times}$ |  4.4 |  32.35 |  -
\+ nerfacc (occ)† |  |  |  4.4 |  32.55 |  0.056
\+ nerfacc (prop)† |  |  |  5.2 |  31.40 |  0.064
Instant-NGP [17] |  360 |  $\mathbf{1.1\times}$ |  5.3 |  25.93 |  -
\+ nerfacc (occ)† |  |  |  5.0 |  26.41 |  0.353
\+ nerfacc (prop)† |  |  |  4.9 |  27.58 |  0.292

(a) Static NeRFs.

Methods |  Dataset |  Speedup |  $\mathbf{T}$ (min) $\downarrow$ |  PSNR $\uparrow$ |  LPIPS $\downarrow$
---|---|---|---|---|---
TiNeuVox [20] |  Hyper. |  $\mathbf{1.7\times}$ |  56.3 |  24.19 |  0.425
\+ nerfacc (occ) |  |  |  33.0 |  24.19 |  0.434
\+ nerfacc (prop) |  |  |  34.3 |  24.26 |  0.398
TiNeuVox [5] |  D-NeRF |  $\mathbf{2.8\times}$ |  11.8 |  31.14 |  0.050
\+ nerfacc (occ) |  |  |  4.2 |  31.75 |  0.038
K-Planes [7] |  D-NeRF |  $\mathbf{1.6\times}$ |  63.9 |  30.28 |  0.043
\+ nerfacc (occ) |  |  |  38.8 |  30.35 |  0.042
T-NeRF [21] |  D-NeRF |  $\mathbf{20\times}$ |  $>$1000 |  28.78 |  0.069
\+ nerfacc (occ)† |  |  |  58.0 |  32.22 |  0.040

(b) Dynamic NeRFs.

Methods |  Dataset |  Speedup |  $\mathbf{T}$ (min) $\downarrow$ |  PSNR $\uparrow$ |  LPIPS $\downarrow$ |  $\mathbf{E_{R}}$ / $\mathbf{E_{T}}$ $\downarrow$
---|---|---|---|---|---|---
BARF [20] |  Syn. |  $\mathbf{4.5\times}$ |  586 |  28.83 |  0.054 |  0.19 / 0.74
\+ nerfacc (occ) |  |  |  130 |  30.11 |  0.044 |  0.07 / 0.35

(c) NeRFs for Camera Optimization. Camera rotational / transnational errors
are denoted as $E_{R}/E_{T}$. $E_{T}$ is scaled by 100.

Table 2: Improvements on Various NeRFs. Experiments are conducted by replacing
necessary code in the official repositories, except for those marked by † that
are based on our re-implementations. Figure 4: Plug-and-play Example in BARF
[14]’s Repository. With around $50$ lines of code change, both training speed
and performance can be greatly improved with efficient sampling equipped in
NerfAcc, as shown in Table 2(b).

###  4.2 Implementation Details

NerfAcc incorporates two advanced sampling methods that can be decoupled from
the radiance field representation, namely the Occupancy Grid from Instant-NGP
[17] and the Proposal Network from Mip-NeRF 360 [2]. Algorithm 1 presents the
pseudo code for volumetric rendering using the NerfAcc toolbox. In this
section, we do not dive into the details of each algorithm, as we basically
follow the original paper’s implementation. Instead, we discuss the
implementation designs that are crucial to maintain high efficiency and
flexibility of this toolbox.

#### Sample as Interval.

In the NerfAcc toolbox, instead of representing each sample with a coordinate
$x$, we represent the sample as an interval along the ray $(t_{0},t_{1},r)$,
where $t_{0}$ and $t_{1}$ are the start and end of the interval along the
$r$-th ray. This interval-based representation offers three key advantages.
Firstly, representing a sample as an interval instead of a single point allows
us to support research based on cone-based rays for anti-aliasing, such as
Mip-NeRF [1] and Mip-NeRF 360 [2]. Secondly, since in almost all cases $t_{i}$
does not require gradients, using $(t_{0},t_{1},r)$ instead of
$(\mathbf{x}_{0},\mathbf{x}_{1})$ to represent the interval allows for the
detachment of the sampling process from the differentiable computational
graph, thereby maximizing its speed. Lastly, the ray id $r$ attached to each
sample enables support for various numbers of samples across rays with a
packed tensor, which we will discuss in the next paragraph. A similar
representation has been adopted in Nerfstudio [29] to support various radiance
fields.

#### Packed Tensor.

To support sampling with spatial skipping, it is necessary to consider that
each ray may result in a different number of valid samples. Storing the data
as a tensor with shape (n_rays, n_samples, ...) and an extra mask with shape
(n_rays, n_samples) to indicate which samples are valid leads to significant
inefficient memory usage when large portions of space are empty. To address
this, in NerfAcc, we represent samples as “packed tensors” with shape
(all_samples, ...), in which only valid samples are stored (see Algo. 1). To
keep track of the associated rays for each sample, we also host an integer
tensor with shape (n_rays, 2), which stores the start index in the packed
tensors and the number of samples on this ray. This approach is similar to
that used in Instant-NGP [17] and PyTorch3D [22].

#### No Gradient Filtering.

After importance sampling, inaccurate transmittance estimations can result in
some samples lying in empty or occluded spaces, particularly in spatial
skipping methods like Occupancy Grid. These samples can be filtered before
being included in PyTorch’s differentiable computation graph by evaluating
their transmittance using the radiance field with gradients disabled. As
backward passes are not required during filtering, this is much faster
(${\sim}10\times$) than keeping all samples in the computation graph. In
practice, samples with transmittance below $10^{-4}$ are disregarded during
this process with almost no influence on the rendering quality. Note that this
strategy is inspired from Instant-NGP [17]’s implementation.

###  4.3 Case Studies

We showcase the flexibility of NerfAcc on three types of NeRFs across seven
papers: static NeRFs (NeRF [16], TensoRF [4], Instant-NGP [17]); dynamic NeRFs
(D-NeRF [21], K-Planes [7] TiNeuVox [5]); and a NeRF variation for camera
optimization (BARF [14]). Although many of these methods, e.g., Instant-NGP,
TensoRF, TiNeuVox and K-Planes, have already been highly optimized for
efficiency, we are still able to accelerate their training by a large margin
and achieve slightly better performance on nearly all cases. It is worth
mentioning that experiments with TensoRF, TiNeuVox, K-Planes and BARF are
conducted by integrating NerfAcc into the official codebase with around 100
lines of code change. The results of our experiments, including those of our
baselines, are presented in Table 2(a),2(b), and2(c), all of which were
conducted under the same physical environment using a single NVIDIA RTX A5000
GPU to facilitate comparison, as per [25]. Aside from the experiments reported
in this paper, NerfAcc has also been integrated into a few popular open-source
projects such as nerfstudio [29] for density-based NeRFs, as well as sdfstudio
[38] and instant-nsr-pl [10] for SDF-based NeRFs.

#### Static NeRFs.

In this task, we experiment with three NeRF vaiants, including the original
MLP-based NeRF [16], TensoRF [4] and Instant-NGP [17]. We show that NerfAcc
works with both MLP-based and Voxel-based radiance fields, on both bounded
(NeRF-Synthetic dataset [16], Tank&Template dataset [12]) and unbounded scenes
(360 Dataset [2]). It is worth to note that with NerfAcc, one can train an
Instant-NGP model with pure Python code and achieve slightly better
performance than the official pure CUDA implementation, as shown in Table
2(a).

#### Dynamic NeRFs.

In this task, we apply the NerfAcc toolbox to T-NeRF [21], K-Planes [7] and
TiNeuVox [5], covering both the synthetic (D-NeRF [21]) and “in-the-wild”
captures111These datasets teleport cameras and do not represent real captures
[9]. (that accompany HyperNeRF [20]). When applying the occupancy grid
approach to accelerate those dynamic methods, instead of representing a static
scene with it, we share the occupancy grid across all frames. In other words,
instead of using it to indicate the opacity of an area at a single timestamp,
We use it to indicate the _maximum opacity at this area over all the
timestamps_. This is not optimal but still makes the rendering very efficient,
due to the fact there are limited movements in these datasets.

#### NeRFs for Camera Optimization.

In this task, we employed the NerfAcc toolbox to BARF [14] on the NeRF-
Synthetic dataset with perturbed cameras. The goal is to jointly optimize the
radiance field and camera extrinsic for multi-view images. We observed that
the spatial skipping sampling provided by NerfAcc facilitated faster training
and significantly improved both image quality and camera pose reconstruction.
These improvements could be attributed to the sparsity enforced in our
sampling procedure. This finding may provide interesting avenues for future
research.

Figure 5: Comparison between Different Sampling Approaches in terms of
Training Time and Rendering FPS. We sweep the hyper-parameters space for each
sampling approach, and find out that occupancy grid approach gives the best
performance-speed trade-off on the NeRF-Synthetic dataset, while the proposal
network approach performs the best on the Mip-NeRF 360 dataset. Note T (min)
denotes for _training_ time and FPS is for _rendering_ frames per second. All
experiments use the HashEncoding from Instant-NGP [17] as the radiance field
representation. Please see the supplementary materials for the hyperparameter
space that we explored.

#### Analysis of Different Sampling Approaches.

Results in Table 2(a) show that the choice between occupancy grid and proposal
network sampling can noticeably affect run-time and performance on different
datasets. As each approach relies on a distinct set of hyperparameters, a
systematic comparison between the two methods is crucial by sweeping the
hyperparameter space. We varied the resolution and marching step size for
occupancy grid and the number of samples and size of the proposal network for
the proposal network approach. We plot pareto curves for each approach for
both the NeRF-Synthetic and Mip-NeRF 360 datasets in Fig.5. This analysis
indicates that occupancy grid sampling is suitable for the NeRF-Synthetic
dataset, whereas the proposal network approach performs better on the 360
dataset. This is likely because the NeRF-Synthetic dataset contains more empty
space that can be skipped effectively using the occupancy grid approach.
However, in the case of real, unbounded data, the use of the occupancy grid
approach is limited by the bounding box and the lack of empty space to skip,
making the proposal network approach more effective. These experiments used
the radiance field from Instant-NGP [17], with the same training recipes.

###  4.4 Combined Sampling

Figure 6: Results of Combined Sampling on NeRF-Synthetic Dataset. Benefit from
our unified concept of transmittance estimator, we find it’s straight-forward
to combine the two distinct sampling approaches.

A benefit from the unified concept of transmittance estimator introduced in
Section 3, is that it’s straight-forward to combine the two distinct sampling
approaches, as both of them essentially provide an estimation of the
transmittance that can be used for importance sampling. For example, we find
that simply stacking an occupancy grid on top of the proposal network, can
significantly reduce the number of rays and shrink the near-far plane for the
remaining rays on the NeRF-Synthetic dataset. This leads to improvements in
quality, from $31.40\text{dB}$ to $32.35\text{dB}$, and a reduction in
training time, from $5.2\text{min}$ to $4.3\text{min}$, compared to using only
the proposal network for importance sampling. Figure 6 shows an example with
the Ficus scene, where the floaters are cleaned out with the combined
sampling. This experiment is conducted using the HashEncoding from Instant-NGP
[17] as the radiance field representation.

##  5 Conclusions

In conclusion, this paper highlights the significant impact of advanced
sampling approaches on improving the efficiency of Neural Radiance Fields
(NeRF) optimization and rendering. We demonstrate that advanced sampling can
significantly speed up the training of various recent NeRF papers, while
maintaining high-quality results. The development of NerfAcc, a flexible
Python toolbox, enables researchers to incorporate advanced sampling methods
into NeRF-related methods easily. The exploration and comparison of advanced
sampling methods are important steps towards developing more efficient and
accessible NeRF-based methods. The presented results also demonstrate the
potential for further research to improve the performance of NeRF and other
related techniques through advanced sampling strategies.

## Acknowledgement

This project was supported in part by the Bakar Fellows Program and the
BAIR/BDD sponsors.

## References

  * [1] Jonathan T Barron, Ben Mildenhall, Matthew Tancik, Peter Hedman, Ricardo Martin-Brualla, and Pratul P Srinivasan.  Mip-nerf: A multiscale representation for anti-aliasing neural radiance fields.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5855–5864, 2021.
  * [2] Jonathan T Barron, Ben Mildenhall, Dor Verbin, Pratul P Srinivasan, and Peter Hedman.  Mip-nerf 360: Unbounded anti-aliased neural radiance fields.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5470–5479, 2022.
  * [3] Jonathan T. Barron, Ben Mildenhall, Dor Verbin, Pratul P. Srinivasan, and Peter Hedman.  Zip-nerf: Anti-aliased grid-based neural radiance fields.  arXiv, 2023.
  * [4] Anpei Chen, Zexiang Xu, Andreas Geiger, Jingyi Yu, and Hao Su.  Tensorf: Tensorial radiance fields.  arXiv preprint arXiv:2203.09517, 2022.
  * [5] Jiemin Fang, Taoran Yi, Xinggang Wang, Lingxi Xie, Xiaopeng Zhang, Wenyu Liu, Matthias Nießner, and Qi Tian.  Fast dynamic radiance fields with time-aware neural voxels.  In SIGGRAPH Asia 2022 Conference Papers, 2022.
  * [6] Julian Fong, Magnus Wrenninge, Christopher Kulla, and Ralf Habel.  Production volume rendering: Siggraph 2017 course.  In ACM SIGGRAPH 2017 Courses, pages 1–79. 2017.
  * [7] Sara Fridovich-Keil, Giacomo Meanti, Frederik Warburg, Benjamin Recht, and Angjoo Kanazawa.  K-planes: Explicit radiance fields in space, time, and appearance.  arXiv preprint arXiv:2301.10241, 2023.
  * [8] Clement Fuji Tsang, Maria Shugrina, Jean Francois Lafleche, Towaki Takikawa, Jiehan Wang, Charles Loop, Wenzheng Chen, Krishna Murthy Jatavallabhula, Edward Smith, Artem Rozantsev, Or Perel, Tianchang Shen, Jun Gao, Sanja Fidler, Gavriel State, Jason Gorski, Tommy Xiang, Jianing Li, Michael Li, and Rev Lebaredian.  Kaolin: A pytorch library for accelerating 3d deep learning research.  https://github.com/NVIDIAGameWorks/kaolin, 2022.
  * [9] Hang Gao, Ruilong Li, Shubham Tulsiani, Bryan Russell, and Angjoo Kanazawa.  Monocular dynamic view synthesis: A reality check.  In Advances in Neural Information Processing Systems, 2022.
  * [10] Yuan-Chen Guo.  Instant neural surface reconstruction, 2022.
  * [11] Yoonwoo Jeong, Seungjoo Shin, and Kibaek Park.  Nerf-factory: An awesome pytorch nerf collection, 2022.
  * [12] Arno Knapitsch, Jaesik Park, Qian-Yi Zhou, and Vladlen Koltun.  Tanks and temples: Benchmarking large-scale scene reconstruction.  ACM Transactions on Graphics, 36(4), 2017.
  * [13] Ruilong Li, Julian Tanke, Minh Vo, Michael Zollhofer, Jurgen Gall, Angjoo Kanazawa, and Christoph Lassner.  Tava: Template-free animatable volumetric actors.  arXiv preprint arXiv:2206.08929, 2022.
  * [14] Chen-Hsuan Lin, Wei-Chiu Ma, Antonio Torralba, and Simon Lucey.  Barf: Bundle-adjusting neural radiance fields.  In IEEE International Conference on Computer Vision (ICCV), 2021\.
  * [15] Ricardo Martin-Brualla, Noha Radwan, Mehdi SM Sajjadi, Jonathan T Barron, Alexey Dosovitskiy, and Daniel Duckworth.  Nerf in the wild: Neural radiance fields for unconstrained photo collections.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7210–7219, 2021.
  * [16] Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng.  Nerf: Representing scenes as neural radiance fields for view synthesis.  Communications of the ACM, 65(1):99–106, 2021.
  * [17] Thomas Müller, Alex Evans, Christoph Schied, and Alexander Keller.  Instant neural graphics primitives with a multiresolution hash encoding.  arXiv preprint arXiv:2201.05989, 2022.
  * [18] Michael Oechsle, Songyou Peng, and Andreas Geiger.  Unisurf: Unifying neural implicit surfaces and radiance fields for multi-view reconstruction.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5589–5599, 2021.
  * [19] Keunhong Park, Utkarsh Sinha, Jonathan T Barron, Sofien Bouaziz, Dan B Goldman, Steven M Seitz, and Ricardo Martin-Brualla.  Nerfies: Deformable neural radiance fields.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5865–5874, 2021.
  * [20] Keunhong Park, Utkarsh Sinha, Peter Hedman, Jonathan T Barron, Sofien Bouaziz, Dan B Goldman, Ricardo Martin-Brualla, and Steven M Seitz.  Hypernerf: A higher-dimensional representation for topologically varying neural radiance fields.  arXiv preprint arXiv:2106.13228, 2021.
  * [21] Albert Pumarola, Enric Corona, Gerard Pons-Moll, and Francesc Moreno-Noguer.  D-nerf: Neural radiance fields for dynamic scenes.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10318–10327, 2021.
  * [22] Nikhila Ravi, Jeremy Reizenstein, David Novotny, Taylor Gordon, Wan-Yen Lo, Justin Johnson, and Georgia Gkioxari.  Accelerating 3d deep learning with pytorch3d.  arXiv preprint arXiv:2007.08501, 2020.
  * [23] Christian Reiser, Richard Szeliski, Dor Verbin, Pratul P Srinivasan, Ben Mildenhall, Andreas Geiger, Jonathan T Barron, and Peter Hedman.  Merf: Memory-efficient radiance fields for real-time view synthesis in unbounded scenes.  arXiv preprint arXiv:2302.12249, 2023.
  * [24] Cheng Sun, Min Sun, and Hwann-Tzong Chen.  Direct voxel grid optimization: Super-fast convergence for radiance fields reconstruction.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5459–5469, 2022.
  * [25] Cheng Sun, Min Sun, and Hwann-Tzong Chen.  Improved direct voxel grid optimization for radiance fields reconstruction.  arXiv preprint arXiv:2206.05085, 2022.
  * [26] Towaki Takikawa, Alex Evans, Jonathan Tremblay, Thomas Müller, Morgan McGuire, Alec Jacobson, and Sanja Fidler.  Variable bitrate neural fields.  In ACM SIGGRAPH 2022 Conference Proceedings, pages 1–9, 2022.
  * [27] Towaki Takikawa, Joey Litalien, Kangxue Yin, Karsten Kreis, Charles Loop, Derek Nowrouzezahrai, Alec Jacobson, Morgan McGuire, and Sanja Fidler.  Neural geometric level of detail: Real-time rendering with implicit 3d shapes.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 11358–11367, 2021.
  * [28] Towaki Takikawa, Or Perel, Clement Fuji Tsang, Charles Loop, Joey Litalien, Jonathan Tremblay, Sanja Fidler, and Maria Shugrina.  Kaolin wisp: A pytorch library and engine for neural fields research.  https://github.com/NVIDIAGameWorks/kaolin-wisp, 2022.
  * [29] Matthew Tancik, Ethan Weber, Evonne Ng, Ruilong Li, Brent Yi, Justin Kerr, Terrance Wang, Alexander Kristoffersen, Jake Austin, Kamyar Salahi, et al.  Nerfstudio: A modular framework for neural radiance field development.  arXiv preprint arXiv:2302.04264, 2023.
  * [30] Dor Verbin, Peter Hedman, Ben Mildenhall, Todd Zickler, Jonathan T Barron, and Pratul P Srinivasan.  Ref-nerf: Structured view-dependent appearance for neural radiance fields.  In 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 5481–5490. IEEE, 2022.
  * [31] John Von Neumann.  13\. various techniques used in connection with random digits.  Appl. Math Ser, 12(36-38):3, 1951.
  * [32] Peng Wang, Lingjie Liu, Yuan Liu, Christian Theobalt, Taku Komura, and Wenping Wang.  Neus: Learning neural implicit surfaces by volume rendering for multi-view reconstruction.  arXiv preprint arXiv:2106.10689, 2021.
  * [33] Zirui Wang, Shangzhe Wu, Weidi Xie, Min Chen, and Victor Adrian Prisacariu.  Nerf–: Neural radiance fields without known camera parameters.  arXiv preprint arXiv:2102.07064, 2021.
  * [34] Lior Yariv, Jiatao Gu, Yoni Kasten, and Yaron Lipman.  Volume rendering of neural implicit surfaces.  Advances in Neural Information Processing Systems, 34:4805–4815, 2021.
  * [35] Alex Yu, Sara Fridovich-Keil, Matthew Tancik, Qinhong Chen, Benjamin Recht, and Angjoo Kanazawa.  Plenoxels: Radiance fields without neural networks.  arXiv preprint arXiv:2112.05131, 2021.
  * [36] Alex Yu, Ruilong Li, Matthew Tancik, Hao Li, Ren Ng, and Angjoo Kanazawa.  Plenoctrees for real-time rendering of neural radiance fields.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5752–5761, 2021.
  * [37] Alex Yu, Vickie Ye, Matthew Tancik, and Angjoo Kanazawa.  pixelnerf: Neural radiance fields from one or few images.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4578–4587, 2021.
  * [38] Zehao Yu, Anpei Chen, Bozidar Antic, Songyou Peng Peng, Apratim Bhattacharyya, Michael Niemeyer, Siyu Tang, Torsten Sattler, and Andreas Geiger.  Sdfstudio: A unified framework for surface reconstruction, 2022.
  * [39] Kai Zhang, Gernot Riegler, Noah Snavely, and Vladlen Koltun.  Nerf++: Analyzing and improving neural radiance fields.  arXiv preprint arXiv:2010.07492, 2020.

##  Appendix A Implementation Details for Case Studies

This section describes the usage of NerfAcc to enhance multiple NeRFs, as
discussed in Section 4.3 of the main paper. We provide detailed information on
how to adapt NerfAcc to various NeRF models.

###  A.1 Static NeRFs.

#### Vanilla NeRF [16].

We trained an 8-layer MLP with the same structure as the vanilla NeRF paper.
The paper employs the PDF-based efficient sampling approach with two MLPs for
course-to-fine sampling (64 + 128). Although the PDF-based approach
effectively concentrates samples around the surface, the MLP itself is
computationally slow. To accelerate the process, we used the spatial skipping
approach with Occupancy Grid provided by NerfAcc. Due to the considerable
amount of empty space in NeRF-Synthetic dataset, we are able to increase the
number of samples per ray to 1024 without any memory issue. As a result, we
achieved a $20\times$ speedup and improved image quality by $+0.5\text{dB}$.

#### TensoRF [4].

Interestingly, incorporating NerfAcc into its own repository results in a
$1.5\times$ speedup for TensoRF [4], a more recent voxel-based method, on both
the NeRF-Synthetic and Tank-and-Temple datasets. We conjecture that this
improvement is due to our update to their transmittance estimator. Although
not discussed in its paper, TensoRF initially employs a skip-based
transmittance estimator when sampling, similar to our occupancy grid. However,
unlike our implementation, TensoRF’s estimator does not maintain an accurate
estimation of the main density field and only updates twice throughout the
entire training process. Our approach, on the other hand, distills a running
average of the density field every few training steps.

#### Instant-NGP [17].

Furthermore, we were able to reproduce Instant-NGP [17] on the NeRF-Synthetic
dataset with the same training speed and slightly better performance
($+0.2\text{dB}$) using the Occupancy Grid. We also achieved significantly
better results on the Mip-NeRF 360 dataset using the proposal network
approach. It is worth noting that the original Instant-NGP implementation is
in pure CUDA with all operations fused into the CUDA kernels. With NerfAcc
handling the underlying sampling logic, the entire training pipeline can be
implemented in Python on top of PyTorch. We believe that our toolbox strikes a
balance between high-performance computing and simplicity, benefiting rapid
research development.

###  A.2 Dynamic NeRFs.

#### T-NeRF [21]

We trained a T-NeRF model described in the D-NeRF paper [21] on the D-NeRF
dataset. By utilizing the occupancy grid in NerfAcc, we achieved a remarkable
speedup of $20\times$ with a corresponding $+3.5\text{dB}$ improvement in
performance. It is worth noting that we performed this experiment using our
re-implementation of T-NeRF model and training recipe.

#### K-Planes [7].

For K-Planes, we achieve around 1.6$\times$ speedup on D-NeRF dataset by
replacing its original proposal-based transmittance estimator with the
occupancy grid. We believe that the sparsity in the D-NeRF dataset makes it
more suitable for spatial skipping approaches.

#### TiNeuVox [5].

Our experiments with TiNeuVox on both datasets suggest a universal speedup by
using our toolbox. Specifically, NerfAcc speeds up training by $2.8\times$ on
D-NeRF and $1.7\times$ on HyperNeRF datasets, respectively. The original
TiNeuVox implementation extends upon DVGO [24] and directly marches rays
within the main density grid, which can be wasteful since space skipping does
not require precise geometry. Our occupancy grid overcomes this issue and has
been shown to be more effective. Additionally, we found that a time-
conditioned proposal network can achieve similar performance on the real
bounded scenes in the HyperNeRF dataset, which aligns with what we observed in
the Instant-NGP experiments with the proposal network.

###  A.3 Camera Optimization.

#### BARF [14].

In the original implementation, BARF uniformly samples within the near-far
range for each ray, analogous to the course level in the vanilla NeRF [16]. We
integrated the occupancy grid from our toolbox for spatial skipping, resulting
in a 4$\times$ speedup during training. Additionally, this led to improved
image quality ($+1.3\text{dB}$) and roughly $2\times$ lower camera
registration error.

##  Appendix B Hyper-parameter Space

Figure 5 in the main paper demonstrates the significant impact of hyper-
parameters on the performance of each sampling approach. The occupancy grid
approach relies on several primary hyper-parameters, such as the binary
threshold $\tau$, the grid resolution $L^{3}$, and the marching step size
$\Delta t$. Meanwhile, the proposal network approach involves hyper-parameters
related to the setup of the proposal network(s) $F_{\Theta}$, as well as the
number of samples $N$ to be taken along each ray. To ensure the robustness of
our experiments, each data point in Figure 5 is trained with a random
combination of hyper-parameters drawn from a reasonable range for each
dataset.


Multi-multifractality and dynamic scaling in stochastic porous lattice

Tushar Mitra and Md. Kamrul Hassan

Theoretical Physics Group, Department of Physics, University of Dhaka, Dhaka 1000, Bangladesh

In this article, we extend the idea of stochastic dyadic Cantor set to weighted planar stochastic
lattice that leads to a stochastic porous lattice. The process starts with an initiator which
we choose to be a square of unit area for convenience. We then define a generator that
divides the initiator or one of the blocks, picked preferentially with respect to their areas,
to divide it either horizontally or vertically into two rectangles of which
one of them is removed with probability $q=1-p$. We find that the remaining number of blocks and
their mass varies with time as $t^{p}$ and $t^{-q}$ respectively. Analytical solution shows that the
dynamics of this process is governed by infinitely many hidden conserved quantities each
of which is a multifractal measure with porous structure as it contains missing blocks of
various different sizes. The support where these measures are distributed is fractal with fractal
dimension $2p$ provided $0<p<1$. We find that if the remaining blocks are characterized
by their respective area then the corresponding block size distribution function obeys dynamic scaling.

61.43.Hv, 64.60.Ht, 68.03.Fg, 82.70Dd

§ INTRODUCTION

We live in an Euclidean space of dimension three and the geometry of objects that we see or interact with are also considered to be Euclidean. Geometry has always been an exceptional tool to describe natural and man-made objects using lines, squares, circles etc. However, not all objects are Euclidean. Rather
there are objects in nature that are accompanied by such an extent of complexity that they
cannot be described by Euclidean geometry and till 1975 they were called geometric
monsters. Benoit B. Mandelbrot conceived, introduced and disseminated the idea of fractal to describe these geometrically complex objects [1, 2]. It has revolutionized
the idea of geometry and has been found useful in almost every branch of academia. There is no
well defined and well accepted definition of fractal. Typically, an object that have
non-integer dimension and have self-similar property is said to be a fractal. However,
recently we have defined it as an object whose dimension $d_f$ is less than that of its embedding
space $d$. In the case, when $d=d_f$ the constituents of the object is uniformly distributed
and as $d_f<d$ it guarantees that the density of constituents is not constant rather decays
following a power-law with the linear system size $L$. It also implies that the smaller the
value of $d_f$ in comparison to $d$, the higher the inhomogeneity of the constituents.
Natural objects like broccoli, snowflake, coastline of an island, thunder etc. can be best described as fractals albeit they are not strictly self-similar rather they are
statistically self-similar. Fractals have been relevant in various subjects like physics, chemistry, biology, earth science, economics, arts etc.  [3, 4, 5]. Recently a group of researchers have created electronic quantum fractal where electron wave functions are delocalized over the Sierpiński structure [6].

One of the simplest text-book example of fractal is definitely the triadic Cantor set.
The process starts with an initiator of closed interval [0,1]. A generator divides it
into three equal parts and deletes the open middle third. The generator is then applied
on the remaining intervals ad infinitum. The resulting set has a non-integer
Hausdorff-Bessicovitch dimension $d_f=\log 2/ \log 3$ which is less than that of the
dimension of the space $d=1$ where the set is embedded. The Cantor set problem has been
relevant in many pedagogical and practical purposes [8, 9].
However, the simple construction of Cantor set lacks in two ways. Firstly, it does not
appear through evolution in time. Secondly, it is not governed by any sort of randomness
throughout its construction process. This two properties are the most common features of
natural fractals. To overcome these Hassan et al. proposed different variants of the
Cantor set like dyadic Cantor set, kinetic and stochastic dyadic Cantor set including
randomness and evolution with time, which has been useful in many disparate systems
[10]. Recently one of us have shown that the conservation law of stochastic
fractals and symmetry is in accordance with Noether's theorem [11, 12].

Constituents that make objects are not always distributed uniformly like we see in
Euclidean and fractal objects. In reality, many quantities fluctuates wildly. For instance, gold, diamond
and many other precious materials are not found at equal concentration everywhere on the Earth.
Rather it is found in high concentrations at only a few places, in low concentrations at many places,
and in very low concentration at almost everywhere on Earth.
To measure the fractal dimension of an object that constitute a set
$S$ we usually partition the embedding space into hyper-cubic cells with
size of their side $\delta$. We then count the number of hyper-cubes $N(\delta)$ that contain at
least one point of the set $S$, varying the size of $\delta$ and measuring the corresponding $N(\delta)$. Slope of the straight line obtained by plotting these data in the $\log-\log$ scale gives the fractal dimension. The important thing that is missing no matter how many points of the set each cell
contain as long as they contain at least one point all are counted with equal weight.
That is, this method does not give any information about detailed nature of the distribution
and hence it is the crudest form of measure. The question is: Is there a way to give a higher weight
to hyper-cubic boxes with higher population of points and a lower weight to those which contain
lower population of points? Multifractal formalism is an answer to this question. In 2010
Hassan et al. proposed a weighted planar stochastic lattice(WPSL) which has been
shown to have many non-trivial properties [13, 14]. One of the properties
is that it is not only a multifractal but infinitely many multifractal spectrum. A multifractal
system is a generalization of a fractal system where a single exponent (the fractal dimension)
is not enough to describe the system, rather a continuous spectrum of exponents are needed [15]. Multifractality has been observed in events like Anderson localization, quantum hall transition and has been useful in many different fields that includes turbulence in fluid dynamics, voltage and current distribution in random resistor networks, heartbeat dynamics, growth by diffusion-limited aggregation (DLA), collision cascade to name a few  [16, 17, 18]. Recently, we have proposed another close variant of WPSL that we named as WPSL1. The lattice provides different properties than that of previously proposed WPSL and is a small-world [20].

In this article, we extend the idea of stochastic dyadic Cantor set to two dimensional case and see how it
differs from its one dimensional counterpart. We start the process with a square of unit area.
To apply on this initiator over and over again we define a generator that
divides a block either horizontally or vertically at random into two smaller block and remove one of them
with probability $1-p$. It is then applied to only block at each time step after picking it
preferentially according its area. Due to removal of blocks we actually create structure that
can be treated as porous media since it contains pores or void of different sizes due to removal of mass.
We can use such porous structure as a backbone while solving theoretical models to gain deeper
insight into various physical phenomena such as percolation. It can also be considered
as fragmentation of two dimensional particles with mass loss. It has many interesting and non-trivial properties.
We show that it is governed by infinitely many conservation laws even though the conservation
of mass law is violated. It is shown that each of these conserved quantity can be regarded
as multifractal measure and hence the structure is not only a multifractal rather infinitely many multifractal. It has been shown that if the blocks of the structure is characterized by their respective
areas then the size distribution function obeys dynamic scaling and this is verified
by using the idea of data-collapse. It suggests that the structure is self-similar in time.

The organization of the rest of this paper is as follows. In section II, the connection between stochastic dyadic Cantor set and weighted planar stochastic lattice is developed. In section III, we propose a porous stochastic lattice by incorporating maass loss in WPSL1 and describe the algorithm. In section IV, various geometric properties of the stochastic porous lattice are discussed and we show how remaining area and number of blocks varies with time $t$. The conservation laws are also discussed in the section. In section V, multifractal analysis of the lattice is performed and we show that it is a infinitely many multifractal. In section VI, the area distribution function is introduced and we show that the lattice possesses temporal self-similarity. Finally, results are discussed and conclusions drawn in section VII.

§ FROM DYADIC CANTOR SET TO WEIGHTED PLANAR STOCHASTIC LATTICE

Weighted planar stochastic lattice was first proposed in 2010 by Hassan et al. which provides many non-trivial topological and geometric properties [13]. The substrate is
typically chosen as a square of unit area and is called an initiator. In step one,
a point on the substrate is chosen at random and two orthogonal lines parallel to its sides are grown to
divide the substrate randomly into four smaller blocks. In the next step and all
the steps thereafter one of the blocks is first picked with probability of their respective area and the generator
divides that randomly into four smaller blocks like it is done in step one to the initiator.
The resulting lattice is called weighted planar stochastic lattice or in short WPSL2
where the factor $2$ refers to the fact that two mutually perpendicular partitioning lines
are used. It has much richer properties than the square lattice and kinetic square lattice [13, 14]. One of the emergent behavior of the lattice is that the coordination number distribution function of WPSL2 follows inverse power-law which immediately implies that the degree distribution of the network
corresponding to its dual follow the same inverse power-law $P(k)\sim k^{-\gamma}$ with the same exponent $\gamma=5.66$ [13].
They have shown that it is governed by infinitely many conservation laws and one of conserved
quantity can be used as multifractal measure and geometrically the lattice is a multifractal. Later it was shown that each of the infinitely many conservation laws is actually a multifractal measure and hence WPSL2 is a multi-multifractal [19].

Recently, we have studied yet another variant of WPSL where a block chosen preferentially with respect to area is divided into two blocks
by dividing it either horizontally or vertically with equal a priori probability. The resulting lattice is called WPSL1, where the factor 1 represents the fact that at each time step only one partitioning line is grown rather than two which was the case for WPSL2. We show WPSL1 too provides awe-inspiring topological and geometric properties. For instance, the degree distribution again follow a power law
$P(k) \sim k^{-\gamma}$ with $\gamma = 4.13$ which is much closer to real life scale-free network.
It has also been found that it has small-world and nested hierarchical network properties. Geometrically WPSL1 also has infinitely many a multi-multifractal spectrum for each of the infinitely many non-trivial conserved quantities although trivial area conservation do not give rise to multifractal spectrum [20].

We are particularly interested in WPSL1. If we compare it with the stochastic dyadic Cantor set(SDCS), we can see that WPSL1 is exactly the two dimensional version of the stochastic dyadic Cantor set, except for the fact that no cell (interval for the case of SDCS) is removed from the initiator in each time step. So, the support of WPSL1 remains an Euclidean two dimensional space, although WPSL1 itself is regarded as a multi-multifractal. We now want to extend WPSL1 to a new model where there is probabilistic mass loss involved, which may introduce a fractal support and can mimic the SDCS in higher dimension.

§ CONSTRUCTION OF STOCHASTIC POROUS LATTICE

A snapshot of the stochastic for p=0.95 and t=5000. The shaded cells indicate that the cells were deleted.

There exist a great extent of demand for a structure that can mimic the porous media. That is,
structure that contain voids or gaps. To this end we now modify the construction process of WPSL1
so that at each time step of its construction process one of the smaller blocks will be
removed with a probability $p$ of our choice.
In step one, the generator divides the initiator randomly into two smaller blocks either horizontally or vertically.
We then label the left block by $a_1$ if it is divided horizontally and the top block
by $a_1$ if it is divided vertically. The other block is always labelled
as $a_2$. The second block $a_2$ is then kept with a probability $p$ and removed with probability $1-p$. In each step thereafter
only one block is picked preferentially with respect to their respective area
and then it is divided in the same way as has been done to the initiator.
In general, the $j$th step of the algorithm can be described as follows.

(i) Subdivide the interval $[0,\sum_{i=1}^{l} a_i]$ into $l$ sub-intervals of size $[0,a_1]$, $[a_1, a_1+a_2]$,$\  ...$,
$[\sum_{i=1}^{l-1} a_i,\sum_{i=1}^{l} a_i]$ each of which represents the blocks labeled by their areas $a_1,a_2,...,a_{l}$ respectively. Here we have assumed that there are $l$ number of remaining blocks in the $j^{th}$ step. Then, the area sub-intervals of the deleted cells or voids fall within the interval $[\sum_{i=1}^{l} a_i,1]$.

(ii) Generate a random number $R$ from the interval $[0,1]$ and find which of the $l$ sub-interval
contains this $R$. The corresponding block it represents, say the $k$th block of area $a_k$, is picked. If $R$ falls within any of the voids left by removal, then increase time by one unit and go back to step (i).

(iii) Calculate the length $x_k$ and the width $y_k$ of this block and keep note of the coordinate of the
lower-left corner of the $k$th block, say it is $(x_{low}, y_{low})$.

(iv) Generate two random numbers $x_R$ and $y_R$ from $[0,x_k]$ and $[0,y_k]$ respectively and hence
the point $(x_{R}+x_{low},y_{R}+y_{low})$ mimics a random point chosen in the block $k$.

(v) Generate a random number $p$ within $[0,1]$.

(vi) If $p<0.5$ then draw a vertical line else a horizontal line through the point $(x_{R}+x_{low},y_{R}+y_{low})$
to divide it into two smaller rectangular blocks. The label $a_k$ is now redundant and hence
it can be reused.

(vii) Label the left or top of the two newly created blocks as $a_{l}$ depending
whether a vertical or horizontal line is drawn respectively and the remaining block is then labeled as $a_{l+1}$.

(viii) Generate another random number $R_1$ within $[0,1]$. If $R_1<p$ then the block $a_{l+1}$ is kept and if $R_1\geq p$, the block is removed.

(ix) Increase time by one unit and repeat the steps (i) - (viii) ad infinitum.

We name the resulting lattice as stochastic porous lattice. A snapshot of the lattice is shown Fig. (<ref>), which reveals the intriguing pattern of surviving blocks and voids. The number of blocks $N(t)$ for the stochastic porous lattice grows with time $t$ but $N(t)$ now becomes probabilistic rather deterministic. The trivial conservation of total area also is no longer applicable for the porous lattice and it too, depends on time $t$.

§ GEOMETRIC PROPERTIES OF STOCHASTIC POROUS LATTICE

In this section we shall try to solve some aspects of the model analytically. To this
end, we  characterize the remaining blocks of the system by their length $x$ and width $y$. We
then define the block size distribution function $f(x,y,t)$ such that
$f(x,y,t)dxdy$ is the number of blocks whose length and width fall within the range
$[x,x+dx]$ and $[y,y+dy]$. Then
the evolution of their distribution function $f(x,y,t)$ can be described by the
master equation

\begin{eqnarray}
\label{eq:WPSL}
{{\partial f(x,y,t)}\over{\partial t}} & = & -xy f(x,y,t)+\dfrac{1+p}{2}x\int_y^\infty f(x,y_1,t)dy_1  \nonumber \\ & +
&  \dfrac{1+p}{2}y\int_x^\infty f(x_1,y,t)dx_1.
\end{eqnarray}

The first term on the right accounts for the loss of blocks of sides $x$ and $y$ due to nucleation of seed
of crack on to it. The pre-factor $xy$ here implies that the blocks are picked
preferentially according to their area.
The second term on the right account for the gain of block of sides $x$ and $y$ upon breaking
the block of sides $x$ and $y_1>y$ horizontally with probability 1/2. The factor $(1+p)$ stands for the fact that on average $1+p$ number of blocks are created with each breaking event. The last term similarly accounts for the contribution due to breaking  a block of sides $x_1>x$ and $y$ vertically.

Solving  Eq. (<ref>) for $f(x,y,t)$ with any given initial size distribution $f(x,y,0)$ is a
formidable task. Therefore, instead of solving for $f(x,y,t)$ we find
it instructive to solve it for $2$-tuple Mellin transform of $f(x,y,t)$ given by

\begin{equation}
\label{momenteq_1}
{M(m,n;t)  =  \int_0^\infty \int_0^\infty x^{m-1} y^{n-1} f(x,y,t) dx dy.
\end{equation}

Incorporating it in Eq. (<ref>) leads us to,

\begin{equation}
\label{momenteq_2}
{{dM(m,n;t)}\over{dt}} = - {{2mn-b(m+n)}\over{2mn}}M(m+1,n+1;t),
\end{equation}

where $b=(1+p)$. It can be seen that if we choose $m=m^\dagger$ and

\begin{equation}
n={{bm^\dagger}\over{2m^\dagger -b}},
\end{equation}

then the corresponding moment is a conserved quantity i.e.

\begin{equation}
\label{momenteq_const}
M\Big(m^\dagger, {{m^\dagger(1+p)}\over{2m^\dagger-(1+p)}};t\Big)=const.,
\end{equation}

$\ \forall \ m^\dagger$ and $\forall\  0<p<1$. Thus, the system is governed by infinitely
many non-trivial conserved quantities.
The credit for such striking feature goes to analytical approach to the problem without which would not
have been possible to even guess that such infinitely many conservation law exist. Unlike
WPSL1 and WPSL2, the present model does not obey conservation of mass law
which is obvious since mass is removed during the process with a non-zero probability.

p     M(1,1;t)     M(2,2;t)

$0.25$     $1.10326 t^{0.25}$     $1.10326 t^{-0.75}$

$0.50$     $1.12838 t^{0.50}$     $1.12838 t^{-0.50}$

$0.75$     $1.08807 t^{0.75}$     $1.08807 t^{-0.25}$

Exact values of the moments as a function of time.

To find the solution for $M(m,n;t)$ we find it convenient to re-write  Eq. (<ref>)

\begin{equation}
\label{momenteq_3}
{{dM(m,n;t)}\over{dt}} = - \Big ( {\alpha_{+}\alpha_{-}\over{m n}}\Big )M(m+1,n+1;t).
\end{equation}

\begin{equation}
\label{momenteq_4}
\alpha_{\pm}={{m+n-b}\over{2}}\pm\sqrt{ {{(m-n)^2}\over{4}}+{{b^2}\over{4}}}.
\end{equation}

On the other hand, the Taylor series expansion of M(m,n;t) about $t = 0$ is

\begin{equation}
\label{eq:taylor}
M(m,n;t)=M(m,n,0)+\sum_{j=1}^\infty {{t^j}\over{j!}}{{d^jM(m,n;t)}\over{dt^j}}\Big |_{t=0}.
\end{equation}

We now iterate Eq. (<ref>) to find various derivatives of $M(m,n;t)$ which are given by,

\begin{eqnarray}
\label{momenteq_5}
\dfrac {d^{j}M(m,n;t)} {dt^{j}} = (-1)^{j} \prod_{k = 1}^{j} \dfrac {(\alpha_{+}+k-1)(\alpha_{-}+k-1)} {(m+k-1) (n+k-1)} \nonumber \\ M(m+j,n+j;t). \nonumber \\
\end{eqnarray}

Using mono-disperse initial condition

\begin{equation}
\end{equation}

gives $M(m,n;0)$ and all its derivatives equal to one. Substituting all these in Eq. (<ref>)
gives a solution for $M(m,n;t)$ as

\begin{equation}
\label{momenteq_6}
M(m,n;t) =  ~_{2}F_{2} (\alpha_{+},\alpha_{-};m,n;-t),
\end{equation}

where $~_{2}F_{2}(a,b;c,d;z)$ is known as generalized hypergeometric function [21].

Using the asymptotic behavior of generalized hypergeometric function we find that the moment $M(m,n;t)$ decays algebraically with time

\begin{equation}
\label{momenteq_7}
M(m,n;t) =  \dfrac {\Gamma(m) \Gamma(n) \Gamma(\alpha_{+}-\alpha_{-})} {\Gamma(\alpha_{+})\Gamma(m-\alpha_{-})\Gamma(n-\alpha_{-})} t^{-\alpha_{-}},
\end{equation}

Using it we immediately find that the first moment that represents the number of intervals increases with
time as

\begin{equation}
\label{momenteq_nc}
M(1,1;t) \sim t^{p} ,
\end{equation}

and the second moment that represents total mass decays with time as

\begin{equation}
\label{momenteq_ac}
M(2,2;t) \sim t^{-(1-p)},
\end{equation}

as expected since the model describes WPSL1 with mass loss.
It can be seen from the above equations the total number of blocks grows following a power law with time with exponent $p$ and total area is not a conserved quantity, rather it decays inversely with time with an exponent $1-p$. The exact values of $M(1,1;t)$ and $M(2,2;t)$ for different $p$ are shown in Table (<ref>) as a function of time $t$ calculated using Eq. (<ref>).

Plot of $ln[(M(1,1;t)]$ vs $ln(t)$ for $p=0.25$, $p=0.5$ and $p=0.75$. The resulting plots are clearly straight line with slopes exactly equal to $p$ for all three cases.

We have performed extensive numerical simulations to verify our results.
First of all, the results based on numerical simulation
for $M(1,1;t)$ and $M(2,2;t)$ are shown in Fig. (<ref>) and Fig. (<ref>) respectively
for $p = 0.25, 0.5$ and $0.75$ in both cases. All the results are extracted from a large number of independent realizations. It can be seen from the plots that the results matches perfectly with the analytical prediction given by Eqs. (<ref>) and (<ref>). Moreover, the fitted straight line in Fig. (<ref>) follows $M(1,1;t)=1.10428t^{0.24961}$ for $p=0.25$, $M(1,1;t)=1.13006t^{0.49987}$ for $p=0.5$, and $M(1,1;t)=1.08655t^{0.75011}$ for $p=0.75$ which matches perfectly with the calculated values given in  Table (<ref>). On the other hand, the fitted straight line in Fig. (<ref>) follows $M(2,2;t)=1.10643t^{-0.74979}$ for $p=0.25$, $M(2,2;t)=1.12731t^{-0.49961}$ for $p=0.5$, and $M(2,2;t)=1.08546t^{-0.24979}$ for $p=0.75$ which also matches perfectly.  This shows the high accuracy and precision between the theory and numerical simulation.

Plot of $ln[(M(2,2;t)]$ vs $ln(t)$ for $p=0.25$, $p=0.5$ and $p=0.75$. The resulting plots are clearly straight lines with slope exactly equal to $-(1-p)$ for all three cases.

Plots of the moment $M(m^\dagger,{{m^\dagger b}\over{2m^\dagger-b}};t)$ vs $t$ for $p=0.5$ and $m^\dagger = 1, 2$ and $4$ for an average of $20000$ independent realizations. It is evident from the plot that the moments are constant with time and equal to $1$ for all values of $m^\dagger$. In the inset the same plot for a single realization where a considerable number of remaining blocks are present in the lattice is shown. It reveals that although $M(m^\dagger,{{m^\dagger b}\over{2m^\dagger-b}};t)$ remains constant with time $t$ for all $m^\dagger$, the constant values are different from one another for a single realization.

To verify the conservation laws given by Eq. (<ref>) we first write its discrete counterpart

\begin{equation}
\label{consteq}
\sum_i^N x_i^{m^\dagger -1} y_i^{ {{m^\dagger b}\over{2m^\dagger-b}}-1}=const.
\end{equation}

To this end we plot $M(m^\dagger,{{m^\dagger b}\over{2m^\dagger-b}};t)$ using Eq.(<ref>) as a function of $t$ in the inset of Fig. (<ref>) for $p=0.5$ $m^\dagger = 1, 2$ and $4$ for a single realization where there are a considerable number of remaining blocks. The moment for each $m^\dagger$ value remains constant with time albeit the value of the moments are different for different values of $m^\dagger$. In fact, the value of the moment for a fixed value of $m^\dagger$ is different in every realization, although it remains constant with time in every case.
In addition, for $m=m^\dagger$ and $n={{bm^\dagger}\over{2m^\dagger -b}}$ the coefficient of $t^{-\alpha_{-}}$ and $\alpha_{-}$ in Eq.(<ref>)  after some calculation turns out to be equal to $1$. This implies that,

\begin{equation}
\label{moment_final}
M\Big(m^\dagger,{{m^\dagger b}\over{2m^\dagger-b}};t\Big) = 1  \ \forall \ m^\dagger > 0, t > 0, p > 0
\end{equation}

The resulting data for $M(m^\dagger,{{m^\dagger b}\over{2m^\dagger-b}};t)$ as a function of $t$ is plotted in Fig. (<ref>) that is averaged over $20000$ independent realizations for $p=0.5$ and $m^\dagger=1, 2$ and $4$. The plot shows that the conserved quantity is equal to 1 independent of the value of $m^\dagger$ and $t$, which confirms our theoretical prediction. The moment for $p=0.25$ and $p=0.75$ were also checked and found to behave same as it is for $p=0.5$. So, the conserved quantity is always equal to 1 for every $p$, $m^\dagger$ and $t$.

§ MULTIFRACTAL ANALYSIS OF STOCHASTIC POROUS LATTICE

In this section, we show that each of the non-trivial conserved quantity can be used as a multifractal
We assume that $\mu_i=x_i^{m^\dagger-1} y_i^{ {{m^\dagger b}\over{2m^\dagger-b}}-1}$ is the
fraction of the total measure $\sum_i^N x_i^{m^\dagger-1} y_i^{ {{m^\dagger b}\over{2m^\dagger-b}}-1}$
which is found in the $i$th block of the remaining blocks. We can now
construct the partition function

\begin{equation}
\end{equation}

of the probability $\mu_i$. Comparing it with the definition of the two-tuple Mellin transform
of $f(x,y,t)$ we can write it as

\begin{equation}
\label{partitioneq}
Z_q=M \Bigg( (m^\dagger-1)q+1, \Big({{m^\dagger b}\over{2m^\dagger-b}}-1\Big)q+1;t\Bigg).
\end{equation}

According to Eq. (<ref>) we, therefore, have the asymptotic solution for the partition function

\begin{eqnarray}
\label{partitioneq1_2}
Z_q(t)\sim t^{\scaleobj{0.8}{\dfrac{\sqrt{(m^\dagger-{{{m^\dagger b}\over{2m^\dagger-b}})^2q^2+b^2}}-\Big((m^\dagger+{{{m^\dagger b}\over{2m^\dagger-b}}-2)q+2-b \Big)}}{2}}}
\end{eqnarray}

We now measure $Z_q$ using the square root of the mean area,

\begin{equation}
\label{delta}
\delta(t) = \sqrt{\dfrac{M(2,2;t)}{M(1,1;t)}} \sim t^{-1/2},
\end{equation}

as the yard-stick and find that it decays following a power-law

\begin{equation}
\label{weightednumber}
Z_q(\delta)\sim \delta^{-\tau(q)},
\end{equation}

where the mass exponent

\begin{eqnarray}
\label{massexponent}
\tau(q) & = & \sqrt{(m^\dagger-{{{m^\dagger b}\over{2m^\dagger-b}})^2q^2+b^2}} - \nonumber\\ &
& \Big((m^\dagger+{{{m^\dagger b}\over{2m^\dagger-b}}-2)q+2-b\Big)}.
\end{eqnarray}

The mass exponent must possess the properties such that $\tau(0)$ is the dimension of the support and $\tau(1)=0$ required by the normalization of the probabilities $p_i$. The later one is clearly true in our case since $\tau(1)=0$. On the other hand, $\tau(0)=2p$, which reveals that the dimension of the
support is not an Euclidean as long as $p<1$. Rather it is a fractal with fractal dimension $2p$
which is always less than the dimension of the space $d=2$ where the lattice is embedded for $p<1$.
This is shown in Fig. (<ref>), where $\tau(q)$ is plotted as a function of $q$ for $p=0.25$
and $m^\dagger=5, 6$ and $7$. It is clear from the plot that $\tau(0)=2p$, in this particular case
$\tau(0)=0.5$ for any value of $m^\dagger$. Also, $\tau(1)=0$ for any value of $m^\dagger$ and $p$.

The plot of $q$ vs $\tau(q)$ for $p=0.75$ and $m^\dagger = 3, 4$ and $5$.
All the plots intersects in two common points one of which represents $\tau(0)=0.5$ and
another one represents $\tau(1)=0$ for any $m^\dagger$.

Often we find that Legendre transform of a function can provide a better physical interpretation than
the function itself. The Legendre transformation of $\tau(q)$ means we use its slope as an independent
parameter instead of $q$ itself. The slope of $\tau(q)$ is known as
Lipschitz-Hölder exponent $\alpha(q)$ and it is defined as

\begin{equation}
\label{eq:alpha_q}
\alpha(q)=-{{d\tau(q)}\over{dq}},
\end{equation}

and hence we find

\begin{eqnarray}
\label{alpha_value}
\alpha(q) & = & \Big( m^\dagger+ {{m^\dagger b}\over{2m^\dagger-b}}-2\Big) -  \nonumber \\ & & \dfrac{q(m^\dagger-{{m^\dagger b}\over{2m^\dagger-b}})^2}{\sqrt{q^{2}(m^\dagger-{{m^\dagger b}\over{2m^\dagger-b}})^2+b^2}} .
\end{eqnarray}

Using the definition of $\alpha(q)$ we find the Legendre transform of $\tau(q)$ which is

\begin{equation}
\label{Legendre}
f(\alpha)=\tau(q) + \alpha q,
\end{equation}

and hence the $f(\alpha)$ spectrum of the porous lattice is

\begin{equation}
f(\alpha(q))=(p-1)+\dfrac{b^2}{\sqrt{q^{2}(m^\dagger-{{m^\dagger b}\over{2m^\dagger-b}})^2+b^2}}.
\end{equation}

It implies that each and every conservation laws obtained by tuning the $m^\dagger$ value results
in a spectrum of spatially intertwined fractal dimensions revealing the fact that the present
system has infinitely many multifractal $f(\alpha)$ spectrum which share the same support. Note that $f(\alpha(q))$ is always concave in character (see figure 6) with a single maximum $2p$ at $q=0$ which corresponds to the dimension of the support. Since the dimension of the support $2p$ is always less than the dimension of the support $d = 2$ for $p<1$, so the support is always a fractal.

The $f(\alpha)$ spectrum for $p=0.75$ and $m^\dagger = 3, 4$ and $5$.
The maximum of each plot occurs at $q=0$ which is in fact the dimension $d_f$ of the support
$2p=1.5$. It confirms the theoretical prediction.

§ AREA SIZE DISTRIBUTION FUNCTION AND DYNAMIC SCALING IN THE STOCHASTIC POROUS LATTICE

The snapshot of the lattice shown in Fig. (<ref>) reveals the complex structure of the lattice in long time limit with a pattern of removed and surviving cells. Although the lattice is highly disordered, we have seen some order and non-trivial properties of the lattice in previous sections. It is then reasonable to ask whether the lattice grows with time in a self-similar fashion, that is, whether two or more snapshots of the lattice taken at different times are similar to each other or not. In Physics,
there exist different kind of self-similarity. An object is called spatial self-similar
if a suitably chosen part of the object is similar to the whole. The litmus test of this property is the
power-law behaviour of an observable. On the other hand, when a system evolves probabilistically
with time and the snapshots taken at different times are look-alike at least in the statistical
sense, then the system is said to exhibit temporal self-similarity. It means that the
numerical values of the dimensional quantities are different in different snapshots but the corresponding dimensionless quantities remain the same. The litmus test of this property is the dynamic scaling [29, 30, 31]

We can characterize each block by their area of the blocks in the lattice and investigate
the nature of their distribution. To this end, we define an observable quantity $C(a,t)da$ as
the number of blocks whose area lies in the range $a$ and $a+da$ at time $t$. To calculate the distribution function we collect data for different times by finding the frequency of number of blocks
using $\delta a$ as an interval size. The distribution function is then normalized by dividing it with the interval size $\delta a$. The total area of the remaining blocks is the given by $\int_0^\infty a C(a,t) da $ which according to Eq.(<ref>) gives,

\begin{equation}
\label{dc_1}
\int_0^\infty a C(a,t) da = \dfrac{1}{\Gamma(1+p)} t^{-(1-p)},
\end{equation}

and the number of remaining blocks is

\begin{equation}
\int_0^\infty C(a,t) da =t^p,
\end{equation}

which are essentially $M(2,2;4)$ and $M(1,1;t)$ respectively.

We shall now apply the Buckingham Pi theorem and invoke the idea of dimensional analysis to
check whether $C(a,t)$ exhibits dynamic scaling or not. The mean area of the remaining blocks $<a(t)>$ varies with time $t$ as $<a(t)> \sim t^{-1}$, as shown in Eq.(<ref>). Now, the governed
parameter $C$ depends on two governing parameters $a$ and $t$. But since $<a(t)> \sim t^{-1}$, it immediately implies that using $t^{-1}$ as a the yard-stick the areas of all the remaining blocks
can be measured. In other words, time $t$ be taken as an independent parameter as area $a$ can
be expressed in terms of $t$ and hence we can define a dimensionless governing parameter,

\begin{equation}
\label{xi}
\xi = \dfrac{a} {t^{-1}} = at.
\end{equation}

It also means $C$ too can be expressed in terms of $t$ alone so that $C\sim t^\theta$
and therefore we can define a dimensionless governed parameter as

\begin{equation}
\label{Pi}
\Pi = \dfrac{C(a,t)}{t^{\theta}}.
\end{equation}

Now, according to the definition of dimensionless quantity, the numerical value of $\Pi$ will remain same even if the time $t$ is scaled by some factor. However $\Pi$ may still depend on the dimensionless governing parameter $\xi$ and hence we can write,

\begin{equation}
\label{Pi_sim}
\Pi \sim \phi(\xi).
\end{equation}

Using Eq. (<ref>) in the above equation
we find that the solution of $C(a,t)$ should have the following dynamic scaling form

\begin{equation}
\label{disfunc}
C(a,t) \sim t^{\theta} \phi (at),
\end{equation}

where $\phi(\xi)$ is called the scaling function. The exponent $\theta$ can be calculated if we substitute Eq. (<ref>) into Eq.(<ref>), which gives,

\begin{equation}
\label{area_eq1}
\int_0^\infty a t^{\theta} \phi(at) da = \dfrac{1}{\Gamma(1+p)} t^{-(1-p)} ,
\end{equation}

which after change of variable becomes,

\begin{equation}
\label{area_eq2}
t^{\theta-2} \int_0^\infty  \xi \phi(\xi) d\xi = \dfrac{1}{\Gamma(1+p)} t^{-(1-p)}.
\end{equation}

The integration over the dimensionless quantity will result into another dimensionless quantity which immediately gives $\theta = 1+p$. So,

\begin{equation}
\label{Ct}
C(a,t) \sim t^{1+p} \phi(at).
\end{equation}

It implies that the plots of $C(a,t)$ versus $a$ for different time $t$ will be distinct. However, if the
same data is plotted for $C(a,t)t^{-(1+p)}$ as a function of $at$ they would collapse into a universal
scaling function since both $C(a,t)t^{-(1+p)}$ and $at$ are dimensionless quantities.

The natural logarithm of the area distribution function $C(a,t)$ as a function of $a$ for $t=5000, 10000$ and $15000$ and $p=0.25$. The resulting plots are straight lines with slopes equal to $t$ and intercepts equal to $t^{1+p}$. The same graph is plotted in the inset of the figure where we have plotted $ln(C(a,t)t^{1+p})$ vs $at$ and it can be seen that all the distinct plots of the main plot collapses into a single universal curves.

We now examine the behavior of the area distribution function $C(a,t)$ to see how it varies with $a$. For this we have taken data for $p=0.25$ for $t=5000, 10000$ and $15000$ which are averaged over a large number of independent realizations. We have plotted $\ln[C(a,t)]$ as a function of $a$ for fixed times $t = 5000, 10000$ and $15000$ as shown in Fig. (<ref>). It is clear from the figure that the corresponding plots are distinct straight lines. The slopes of these straight lines were found to be equal to the time
$t=5000, 10000$ and $15000$ when the snapshots are taken and it confirms that,

\begin{equation}
\label{scaling}
C(a,t) \sim e^{-at}.
\end{equation}

The intercepts were found to be equal $10.682,11.545$ and $12.047$ for $t=5000,10000$ and $15000$ respectively. These intercepts are equal to $ln(t^{1+p})$ for all the three values of $t$ which is in accordance with Eq. (<ref>). So, we can write that,

\begin{equation}
\label{area_final}
C(a,t) \sim t^{1+p} e^{-at},
\end{equation}

and hence the solution for the scaling function is $\phi(\xi) \sim e^{-\xi}$.

Now, we examine whether the area distribution function possesses self-similarity using the idea of data collapse. To do this $C(a,t)$ is measured using $t^{(1+p)}$ as an yardstick and a is measured using $t^{-1}$ as an yardstick which makes them dimensionless governed quantity and dimensionless governing quantity respectively. As we have discussed earlier, the dimensionless quantities should remain the same for all $t$. We have plotted $\ln(C(a,t)t^{-(1+p)})$ vs at in the inset of Fig. (<ref>) and it can be seen that all the distinct plots of the main plot of $ln(C(a,t)t^{1+p})$ vs $at$ collapses into a single universal plot ignoring the fat tails. This universal plot now contains all the information about $C(a,t)$ in any time including the infinite time limit. This proves that although the numerical values of the dimensional quantities are different for each different time $t$, the corresponding dimensionless quantities remain invariant. So C(a,t) obeys dynamic scaling and the snapshots at different time $t$
are self-similar from the perspective of area size distribution.

§ SUMMARY AND CONCLUSIONS

In this article, we have formulated a bridge between stochastic dyadic Cantor set and weighted
planar stochastic lattice with one crack (WPSL1). We have shown that if probabilistic mass loss
is incorporated in WPSL1 such that a cell is kept with probability $p$ and removed with
probability $1-p$, it can mimic a lattice with porous character. To the best of our
knowledge there do not exist a
random lattice with coordination number disorder and have porosity of different sizes
too. We have given a generalized master equation and solved it to find
the various moments of the block size distribution functions. The trivial conserved
quantity of WPSL1, that is the total mass is no longer conserved rather it decreases
as $t^{-(1-p)}$. On the other hand, the total number of
blocks increases with time as $t^{p}$ revealing that the mean area of the block size decreases
as $t^{-1}$. The numerical results matched perfectly and are consistent with the analytical
results in all cases.

The stochastic porous lattice also possesses infinitely many hidden
conservation laws for each value of $p$ which are highly non-trivial
since we can only know about them thanks to analytical solution. Moreover the ensemble
average value of each of
the conserved quantities was found to be equal to $1$ and it was
also proven analytically albeit they are different from $1$ in every independent relization.
However, each of the infinitely many conservation laws is a
multifractal measure, and so WPSL1 with mass loss possesses multiple multifractality.
We have obtained an exact analytical generalized for the $f(\alpha)$ spectrum. It shows that
the peaks of $f(\alpha)$ spectrum is equal to $2p$ which is the support of the
multifractal measure and it is fractal provided $0<p<1$. Finding an exact analytical expression
for $f(\alpha)$ is often formidable task.

We have checked whether the snapshot of the complex structure of the stochastic porous lattice
taken at different times are similar or not. To this end, we have used the area distribution function and incorporated
dynamic scaling and Buckingham Pi theorem. The scaling solution of the area distribution
function $C(a,t)$ is found to exhibit dynamic scaling $C(a,t) \sim t^{1+p} e^{-at}$. To show
this we have used the idea of
data-collapse. We have shown that the distinct plots of $C(a,t)$ as a function
of $a$ for the snapshots at different time collapses superbly into on universal
curve if we plot $t^{-(1+p)}C(a,t)$ versus $at$ instead.
It clearly suggests that the porous structure of the stochastic lattice
possesses temporal self-similarity.

A two-dimensional counterpart of stochastic dyadic Cantor set leading to a porous structure
has the potential of great use studying spreading of virus (computer or biological), fluids, opinion,
rumors etc. Besides, it can also be used as a skeleton for percolation process to see the impact of
$p$ that determines the degree of porosity. We can also see if percolation on it belongs to a
new universality class or to the same class where WPSL belong to regardless of the $p$ value.
We intend to address some of these issues in our future endeavors.

[1] B. B. Mandelbrot, Fractals: Form, Chance, and Dimen-
sion, Freeman, San Francisco, 1977.
[2] B. B. Mandelbrot, The Fractal Geometry of Nature, Free-
man, San Francisco, 1982.
[3] Edgar E. Peters, Fractal Market Analysis: Applying
Chaos Theory to Investment and Economics, Wiley, New
York, 1994.
[4]  G. Korvin, Fractal models in the earth sciences, Amster-
dam, Elsevier, 1992.
[5] T. Vicsek, Fractal Growth Phenomena, Second ed.,
World Scientific, Singapore, 1992.
[6] S. N. Kempkes, M. R. Slot, S. E. Freeney, S. J. M. Zevenhuizen, D. Vanmaekelbergh, I. Swart and C. Morais Smith, Nat. Phys. 15, 127–131 (2019).
[7] J. Feder, Fractals, Plenum Press, New York, 1988.
[8] S. Sears, M. Soljacic, M. Segev, D. Krylov, and D.
Bergman, Cantor set fractals from solitons, Phys. Rev. Lett. 84 1902 (2000).
[9] P. L. Krapivsky and S. Redner, Random walk with shrinking steps, Am. J. Phys. 72 591 (2004).
[10] Hassan MK, Pavel NI, Pandit RK, Kurths J, Chaos Solit Fract 60, 31 (2014).
[11] M. K. Hassan, Eur Phys J Spec Top 228 209 (2019).
[12] Rakibur Rahman, Fahima Nowrin, M. S. Rahman, J. A. Wattis, M. K. Hassan, Phys Rev. E 103 022106 (2021).
[13] M. K. Hassan, M. Z. Hassan and N. I. Pavel, New J. Physics, 12 093045 (2010).
[14] M. K. Hassan, M. Z. Hassan and N. I. Pavel, J. Phyis. Conf. Ser. 297 012010 (2011).
[15] H. David, Multifractals Theory and Applications, London, Chapman and Hall, 2001.
[16] S. Faez, A. Strybulevych, J. H. Page, A. Lagendijk, B. A. van Tiggelen, Phys. Rev. Lett. 103, 155703 (2009).
[17] F. Evers , AD. Mirlin, Rev Mod Phys 80, 1355 (2008).
[18] P. C. Ivanov, L. A. N. Amaral, A. L. Goldberger, S. Havlin, M. G. Rosenblum, Z. Struzik, H. E. Stanley, Nature, 399 , 461-465 (1999).
[19] F. R. Dayeen, M. K. Hassan, Chaos, Solitons and Fractals, 91 228 (2016)
[20] T. Mitra, M. K. Hassan, arXiv:2107.00936 [cond-mat.stat-mech]
[21] Luke Y L 1969 The special functions and their approximations I
(New York: Academic Press)
[22] M. K. Hassan and G. J. Rodgers, Phys. Lett. A 218 207 (1996).
[23] G. J. Rodgers and M. K. Hassan, Phys. Rev. E 50 3458 (1994).
[24] P. L. Krapivsky, S. Redner, and E. Ben-Naim, 2010
A kinetic view of Statistical Physics (Cambridge University Press, New York)
[25] K. Klemm and V. M. Eguiluz, Phys. Rev. E 65, 036123 (2002).
[26] P. Holme and B. J. Kim, Phys. Rev. E 65, 026107 (2002).
[27] H. Zhou, Phys. Rev. E 66, 016125 (2002).
[28] M. A. Serrano and M. Boguna, Phys. Rev. E 72, 063133 (2005).
[29] Hassan MK, Hassan MZ, Islam N, Phys. Rev. E 88, 042137 (2013)
[30] Hassan MK, Hassan MZ, Pavel NI, J Phys A: Math Theor 44, 175101 (2011)
[31] T. Mitra, T. Hossain, S. Banerjee, M. K. Hassan, Chaos Solit Fract 145, 110790 (2021)


* [22] M. Doi and S. F. Edwards.  The Theory of Polymer Dynamics, volume 73.  Oxford University Press, 1986.
  * [23] I. Frankel and H. Brenner.  On the foundations of generalized Taylor dispersion theory.  Journal of Fluid Mechanics, 204:97–119, 1989.
  * [24] T. Gallay and M. C. Zelati.  Enhanced dissipation and Taylor dispersion in higher-dimensional parallel shear flows.  arXiv preprint arXiv:2108.11192, 2021.
  * [25] E. Grenier, T. T. Nguyen, and I. Rodnianski.  Landau damping for analytic and Gevrey data.  arXiv preprint arXiv:2004.05979, 2020.
  * [26] Y. Guo.  The Landau equation in a periodic box.  Comm. Math. Phys., 231(3):391–434, 2002.
  * [27] C. Hohenegger and M. J. Shelley.  Stability of active suspensions.  Physical Review E, 81(4):046311, 2010.
  * [28] G. B. Jeffery.  The motion of ellipsoidal particles immersed in a viscous fluid.  Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character, 102(715):161–179, 1922.
  * [29] N. Jiang, Y.-L. Luo, and T.-F. Zhang.  Coupled self-organized hydrodynamics and Navier–Stokes models: local well-posedness and the limit from the self-organized kinetic-fluid models.  Archive for Rational Mechanics and Analysis, 236(1):329–387, 2020\.
  * [30] L. Kanzler and C. Schmeiser.  Kinetic model for myxobacteria with directional diffusion.  arXiv preprint arXiv:2109.13184, 2021.
  * [31] D. L. Koch and G. Subramanian.  Collective hydrodynamics of swimming microorganisms: living fluids.  Annual Review of Fluid Mechanics, 43:637–659, 2011.
  * [32] J. La.  Global well-posedness of strong solutions of Doi model with large viscous stress.  Journal of Nonlinear Science, 29(5):1891–1917, 2019.
  * [33] J. La.  On diffusive 2d Fokker–Planck–Navier–Stokes systems.  Archive for Rational Mechanics and Analysis, 235(3):1531–1588, 2020\.
  * [34] E. Lauga.  The fluid dynamics of cell motility, volume 62.  Cambridge University Press, 2020.
  * [35] P.-L. Lions and N. Masmoudi.  Global existence of weak solutions to some micro-macro models.  Comptes Rendus Mathematique, 345(1):15–20, 2007.
  * [36] A. Manela and I. Frankel.  Generalized Taylor dispersion in suspensions of gyrotactic swimming micro-organisms.  Journal of Fluid Mechanics, 490:99–127, 2003.
  * [37] N. Masmoudi.  Equations for polymeric materials.  In Handbook of Mathematical Analysis in Mechanics of Viscous Fluids, pages 973–1005. Springer International Publishing, 2018.
  * [38] N. Masmoudi, P. Zhang, and Z. Zhang.  Global well-posedness for 2d polymeric fluid models and growth estimate.  Physica D: Nonlinear Phenomena, 237(10-12):1663–1675, 2008.
  * [39] C. Mouhot and C. Villani.  On Landau damping.  Acta Math., 207(1):29–201, 2011.
  * [40] L. Ohm and M. J. Shelley.  Weakly nonlinear analysis of pattern formation in active suspensions.  Journal of Fluid Mechanics, to appear, 2022.
  * [41] F. Otto and A. E. Tzavaras.  Continuity of velocity gradients in suspensions of rod–like molecules.  Communications in Mathematical Physics, 277(3):729–758, 2008.
  * [42] W. Rudin.  Real and complex analysis.  McGraw-Hill Book Co., New York, third edition, 1987.
  * [43] D. Saintillan and M. J. Shelley.  Orientational order and instabilities in suspensions of self-locomoting rods.  Physical review letters, 99(5):058102, 2007.
  * [44] D. Saintillan and M. J. Shelley.  Instabilities and pattern formation in active particle suspensions: kinetic theory and continuum simulations.  Physical Review Letters, 100(17):178103, 2008.
  * [45] D. Saintillan and M. J. Shelley.  Instabilities, pattern formation, and mixing in active suspensions.  Physics of Fluids, 20(12):123304, 2008.
  * [46] D. Saintillan and M. J. Shelley.  Emergence of coherent structures and large-scale flows in motile suspensions.  Journal of the Royal Society Interface, 9(68):571–585, 2012.
  * [47] D. Saintillan and M. J. Shelley.  Active suspensions and their nonlinear models.  Comptes Rendus Physique, 14(6):497–517, 2013.
  * [48] D. Saintillan and M. J. Shelley.  Theory of active suspensions.  In Complex Fluids in biological systems, pages 319–355. Springer, 2015.
  * [49] O. Sieber.  Existence of global weak solutions to an inhomogeneous Doi model for active liquid crystals.  arXiv preprint arXiv:2006.16832, 2020.
  * [50] E. M. Stein and T. S. Murphy.  Harmonic analysis: real-variable methods, orthogonality, and oscillatory integrals, volume 3.  Princeton University Press, 1993.
  * [51] G. Subramanian and D. L. Koch.  Critical bacterial concentration for the onset of collective swimming.  Journal of Fluid Mechanics, 632:359–400, 2009.
  * [52] G. Subramanian, D. L. Koch, and S. R. Fitzgibbon.  The stability of a homogeneous suspension of chemotactic bacteria.  Physics of Fluids, 23(4):041901, 2011.
  * [53] Y. Sun and Z. Zhang.  Global well-posedness for the 2d micro-macro models in the bounded domain.  Communications in mathematical physics, 303(2):361–383, 2011.
  * [54] G. I. Taylor.  The dispersion of matter in turbulent flow through a pipe.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 223(1155):446–468, 1954.
  * [55] C. Villani et al.  Landau damping.  Notes de cours, CEMRACS, 2010.
  * [56] D. Wei.  Diffusion and mixing in fluid flow via the resolvent estimate.  Sci. China Math., 64(3):507–518, 2021.
  * [57] D. Wei and Z. Zhang.  Enhanced dissipation for the Kolmogorov flow via the hypocoercivity method.  Sci. China Math., 62(6):1219–1232, 2019.
  * [58] H. Zhang and P. Zhang.  On the new multiscale rodlike model of polymeric fluids.  SIAM journal on mathematical analysis, 40(3):1246–1271, 2008.