# Quantum self-bound droplets in Bose-Bose mixtures: Effects of higher-order
quantum and thermal fluctuations

Nadia Guebli 1 and Abdelâali Boudjemâa 2,3 1 Laboratory for Theoretical
Phyiscs and Material Physics, Hassiba Benbouali University of Chlef, P.O. Box
78, 02000, Ouled-Fares, Chlef, Algeria.
2 Department of Physics, Faculty of Exact Sciences and Informatics, Hassiba
Benbouali University of Chlef, P.O. Box 78, 02000, Ouled-Fares, Chlef,
Algeria.
3Laboratory of Mechanics and Energy, Hassiba Benbouali University of Chlef,
P.O. Box 78, 02000, Ouled-Fares, Chlef, Algeria<EMAIL_ADDRESS>
###### Abstract

We systematically study the effects of higher-order quantum and thermal
fluctuations on the stabilization of self-bound droplets in Bose mixtures
employing the time-dependent Hartree-Fock-Bogoliubov theory. We calculate the
ground-state energy, the droplet equilibrium density, the depletion and
anomalous density of the droplet as well as the critical temperature as a
function of the relevant parameters. Our findings are compared with previous
analytical predictions and diffusion Monte Carlo simulations. We employ our
theory together with the local density approximation for quantum and thermal
fluctuations to obtain an extended finite-temperature Gross-Pitaevskii
equation. The density profiles and breathing modes of the droplet are deeply
examined in terms of the interaction strength and the temperature by
numerically solving the developed generalized Gross-Pitaevskii equation.

##  I Introduction

Recently, quantum self-bound droplets in Bose-Bose mixtures have attracted
significant theoretical and experimental interest Petrov ; Capp ; Cik ; Cab ;
Sem ; Err . This exotic state of matter could pave the way to investigate
beyond-mean-field intriguing quantum phenomena. It arises from the competition
of an attractive mean-field energy and repulsive beyond-mean-field effects
provided by the Lee-Huang-Yang (LHY) term. This latter stems from the quantum
fluctuations induced by interactions. Self-bound droplets in Bose mixtures
have been the object of several recent works. The role played by dipole-dipole
interactions in the stability of the droplets in binary Bose-Einstein
condensates (BECs) was first highlighted by one of us Boudj . The effects of
correlations in self-bound Bose mixtures have been explored using the
variational hypernetted-chain Euler-Lagrange method Stau . The dynamical
properties of low-dimensional self-bound quantum droplets of bosonic mixtures
have been studied in Refs.Astr ; Chiq ; Astr1 ; Astr2 ; Astr3  by solving the
corresponding amended Gross-Pitaevskii equation (GPE). Most recently, the
self-bound supersolid stripe phase in binary BEC has been predicted in
Ref.Sach .

The aim of this paper is to investigate the effects of quantum and thermal
fluctuations on the occurrence and stabilization of quantum droplets in
symmetric Bose mixtures using the beyond LHY description based on our time
dependent Hartree-Fock-Bogoliubov (TDHFB) theory Boudj ; Boudj1 ; Boudj2 ;
Boudj3 ; Boudj4 ; Boudj5 ; Boudj6 ; Boudj7 ; Boudj8 ; Boudj9 ; Boudj10 ;
Boudj11 ; Boudj12 ; Boudj13 . This latter includes automatically higher-order
corrections originating from normal and anomalous correlations to the equation
of state generalizing all existing models in the literature. The coupling
between the order parameter, and the normal and anomalous fluctuations makes
the TDHFB a promising approach for exploring quantum self-bound droplets in a
dipolar single BEC Boudj12  and Bose-Bose mixtures Boudj13  at both zero and
finite temperatures.

In this paper we extract useful analytical results for the ground-state energy
and the equilibrium density that extend naturally the seminal equations of
Petrov Petrov . Explicit expressions connecting normal and anomalous
correlations to the droplet equilibrium density are also obtained for the
first time to the best of our knowledge. We show that at zero temperature, our
theory not only provides fascinating results but also captures genuine higher-
order quantum effects predicted from recent diffusion Monte Carlo (DMC)
simulations Cik1 . At finite temperature, we calculate the free energy and
precisely determine the thermal equilibrium density of the droplet. We find
that the droplet destabilizes when the temperature becomes slightly larger
than the ground-state energy of the droplet, in good agreement with recent
results inferred from the macroscopic approach Ota  and the pairing theory Hu
. Our predictions point out that higher-order corrections and interspecies
interactions may shift the critical temperature.

On the other hand, we deal with the ground-state properties of the droplet
such as the density profiles and the collective modes. These latter play a
crucial role for probing many-body effects in ultracold quantum gases. Our
calculations are based on the extended finite-temperature GPE which is derived
in a self-consistent manner from our formalism taking into account higher-
order quantum and thermal corrections under the local density approximation.
Our goal here is to understand how the density distribution and the collective
oscillations behave when varying the interaction strength and temperature. To
this end, we solve numerically and variationally the developed GPE.

The rest of the paper is organized as follows. In Sec.II, we recapitulate the
TDHFB theory. We briefly review the expressions of the free energy, and the
normal and anomalous fluctuations. In Sec.III, we present our calculations up
to second-order terms for the energy per particle and equilibrium density. Our
results are compared with available theoretical treatments and DMC
simulations. The behavior of the noncondensed and anomalous densities of the
droplet is also highlighted. Furthermore, we look at how the thermal
fluctuations influence the nucleation and stability of the self-bound droplet.
We compute also the critical temperature and thermal equilibrium density of
the self-bound droplet. Section IV is devoted to the numerical analysis of the
ground-state properties of the self-bound droplet. We calculate in particular
diverse density profiles and the collective modes. The paper ends with
conclusions and an outlook in Sec.V.

##  II Model

We consider a weakly interacting two-component BEC with an atomic mass
$m_{j}$. In the mean field, the physics of this system is given by the
following energy functional which includes the LHY quantum corrections and
normal and anomalous densities Boudj ; Boudj5 ; Boudj13 ,

$\displaystyle{\cal E}$ $\displaystyle=\sum_{j=1}^{2}\bigg{[}\int d{\bf
r}\,\left(\Phi_{j}^{*}h_{j}^{sp}\Phi_{j}+\hat{\bar{\psi}}_{j}^{\dagger}h_{j}^{sp}\hat{\bar{\psi}}_{j}+\frac{g_{j}}{2}n_{j}^{2}\right)\bigg{]}$
(1) $\displaystyle+g_{12}\int d{\bf r}n_{1}n_{2}+{\cal E}_{\text{LHY}},$

where $h_{j}^{sp}=-(\hbar^{2}/2m_{j})\Delta-\mu_{j}$ is the single particle
Hamiltonian, $\mu_{j}$ represent the chemical potentials related to each
component, $n_{cj}=|\Phi_{j}|^{2}$ is the condensed density, and
$\hat{\bar{\psi}}_{j}({\bf r})=\hat{\psi}_{j}({\bf r})-\Phi_{j}({\bf r})$ is
the noncondensed part of the field operator with $\Phi_{j}({\bf
r})=\langle\hat{\psi}_{j}({\bf r})\rangle$. The coefficients
$g_{j}=(4\pi\hbar^{2}/m_{j})a_{j}$ and
$g_{12}=g_{21}=2\pi\hbar^{2}(m_{1}^{-1}+m_{2}^{-1})a_{12}$ with $a_{j}$ and
$a_{12}$ being the intraspecies and the interspecies scattering lengths,
respectively. In the spirit of our HFB theory, the LHY corrections to the
energy are defined as:

$\displaystyle{\cal E}_{\text{LHY}}$
$\displaystyle=\frac{1}{2}\sum_{j=1}^{2}g_{j}\int d{\bf
r}\big{(}2\tilde{n}_{j}n_{j}-\tilde{n}_{j}^{2}+|\tilde{m}_{j}|^{2}$ (2)
$\displaystyle+\tilde{m}_{j}^{*}\Phi_{j}^{2}+\tilde{m}_{j}{\Phi_{j}^{*}}^{2}\big{)},$

where
$\tilde{n}_{j}=\langle\hat{\bar{\psi}}_{j}^{\dagger}\hat{\bar{\psi}}_{j}\rangle$
is the noncondensed density, and
$\tilde{m}_{j}=\langle\hat{\bar{\psi}}_{j}\hat{\bar{\psi}}_{j}\rangle$ is the
anomalous density. The presence of these quantities enables us to derive the
LHY term without any ad-hoc assumptions in contrast to the standard
generalized GPE.

Let us start by briefly reviewing the main results predicted by our theory
Boudj13  that permit us to study the properties of the self-bound droplets in
Bose mixtures.

The free energy can be computed by resorting to the dimensional
regularization. This yields the famous $T^{4}$-law Boudj13

$\displaystyle\frac{F}{V}=$
$\displaystyle\frac{1}{2}\sum_{j=1}^{2}g_{j}n_{j}^{2}+g_{12}n_{1}n_{2}+\frac{16g_{1}\sqrt{a_{1}^{3}/\pi}}{15\sqrt{2}}n_{1}^{5/2}$
$\displaystyle\times\bigg{(}1+\frac{\tilde{m}_{1}-\tilde{n}_{1}}{n_{1}}\bigg{)}^{5/2}\Bigg{\\{}\sum_{\pm}\bigg{[}f_{\pm}^{5/2}(\Delta,\alpha)$
$\displaystyle-\frac{1}{2\sqrt{2}}\frac{f_{\pm}^{-3/2}(\Delta,\alpha)}{[1+(\tilde{m}_{1}-\tilde{n}_{1})/n_{1}]^{4}}\left(\frac{\pi
T}{g_{1}n_{1}}\right)^{4}\bigg{]}\Bigg{\\}},$ (3)

where $V$ is the volume of the system,
$f_{\pm}(\Delta,\alpha)=1+\alpha\pm\sqrt{(1-\alpha)^{2}+4\Delta^{-1}\alpha}$,
$\alpha=\bar{g}_{2}n_{2}/\bar{g}_{1}n_{1}$, and
$\Delta=\bar{g}_{1}\bar{g}_{2}/g_{12}^{2}$ is the miscibility parameter. The
last term which accounts for the effect of the thermal fluctuations, is
calculated at temperatures $T\ll g_{1}n_{1}$, where the main contribution to
Eq.(II) comes from the phonon region. In the limit of lower density
$n\rightarrow 0$, the thermal corrections to the free energy diverge as
$n^{-3/2}$ and vanish at zero temperature. A similar formula has been derived
in Ref.Ota  using the beyond-LHY theory, based on the calculation of second-
order terms in the Bogoliubov phonon modes. At $T=0$ and for
$\tilde{m}_{j}=\tilde{n}_{j}=0$, the free energy (II) reproduces the
analytical expression for the ground-state energy Petrov .

Our formalism provides useful expressions for the noncondensed and anomalous
densities at finite temperature Boudj13 ,

$\displaystyle\tilde{n}$
$\displaystyle=\frac{2\sqrt{2}}{3}\,n_{1}\sqrt{\frac{n_{1}a_{1}^{3}}{\pi}}\bigg{(}1+\frac{\tilde{m}_{1}-\tilde{n}_{1}}{n_{1}}\bigg{)}^{3/2}$
(4)
$\displaystyle\times\sum_{\pm}\bigg{[}f_{\pm}^{3/2}(\Delta,\alpha)+\frac{1}{3}\frac{f_{\pm}^{-1/2}(\Delta,\alpha)}{[1+(\tilde{m}_{1}-\tilde{n}_{1})/n_{1}]^{2}}\left(\frac{\pi
T}{g_{1}n_{1}}\right)^{2}\bigg{]},$

and

$\displaystyle\tilde{m}$
$\displaystyle=2\sqrt{2}\,n_{1}\sqrt{\frac{n_{1}a_{1}^{3}}{\pi}}\bigg{(}1+\frac{\tilde{m}_{1}-\tilde{n}_{1}}{n_{1}}\bigg{)}^{3/2}$
(5)
$\displaystyle\times\sum_{\pm}\bigg{[}f_{\pm}^{3/2}(\Delta,\alpha)+\frac{1}{3}\frac{f_{\pm}^{-1/2}(\Delta,\alpha)}{[1+(\tilde{m}_{1}-\tilde{n}_{1})/n_{1}]^{2}}\left(\frac{\pi
T}{g_{1}n_{1}}\right)^{2}\bigg{]}.$

For $g_{12}=0$, Eqs.(4) and (5) reduce to those obtained in our Ref.Boudj15
for a single component BEC. The expression of the anomalous density (5) has
been indeed obtained using the dimensional regularization Boudj3 ; Boudj12 ;
Anders  to overcome the ultraviolet divergence. It is worth stressing that the
involvement of the anomalous density in Bose systems leads to a double
counting of the interaction effects Boudj2 ; Boudj10 ; Griffin ; DStoof . The
solution of the self-consistent Eqs.(4) and (5) requires the use of some form
of an iterative scheme.

Figure 1: (Color online) (a) Ground state energy as a function of the density
for $a_{12}/a=-1.1$. (b) The same as (a) but for $a_{12}/a=-1.05$. Red dashed
lines correspond to the pairing theory Hu . Black dotted-dashed lines
correspond to our results to first-order corrections of quantum fluctuations.
Black dotted lines correspond to our results up to second-order corrections of
quantum fluctuations. Black solid lines correspond to Petrov’s predictions
Petrov . Circles represent the DMC results of Ref.Cik1 .

##  III Self-bound droplet

Here, we apply the general method presented in the previous section to
investigate the properties of the self-bound droplet at both zero and finite
temperatures in the presence of higher-order quantum effects.

From now on, we consider symmetric Bose-Bose mixtures with equal intraspecies
coupling constants $g_{1}=g_{2}=g$, and equal densities in each component
$n_{1}=n_{2}=n/2$, $\tilde{m}_{1}=\tilde{m}_{2}=\tilde{m}/2$, and
$\tilde{n}_{1}=\tilde{n}_{2}=\tilde{n}/2$.

###  III.1 Zero-temperature case

Figure 2:  (a) Equilibrium density of the droplet as a function of $a_{12}/a$.
(b) Equilibrium density in units of $n_{0}$, as a function of $a_{12}/a$.
Black solid lines correspond to Petrov’s equilibrium density Petrov . Black
dotted-dashed lines correspond to our results from Eq. (7). Magenta dashed
lines correspond to the pairing theory of Ref.Hu . Blue circles correspond to
DMC data Cik1 .

The ground state energy can be rewritten in the following dimensionless form:

$\displaystyle\frac{E}{NE_{0}}=$ $\displaystyle\pi\left(\frac{\delta
a}{a}\right)_{+}(na^{3})+\frac{32\sqrt{\pi}}{15}(na^{3})^{3/2}$
$\displaystyle\times\bigg{(}1+\frac{\tilde{m}-\tilde{n}}{n}\bigg{)}^{5/2}\sum_{\pm}\left(\frac{\delta
a}{a}\right)_{\pm}^{5/2},$ (6)

where $E_{0}=\hbar^{2}/ma^{2}$, and $(\delta a/a)_{\pm}=1\pm(a_{12}/a)$.
Equation (III.1) extends and unifies various equations of state in the
existing literature. For $\left(\delta a/a\right)_{+}<0$, the mean-field
theory energy provides a term $\propto n^{2}$ leading to a collapse of a
homogeneous state towards bright soliton formation. The beyond-mean-field LHY
term adds an extra repulsive term $\propto n^{5/2}$ compensating the
attractive mean-field term, enabling quantum and thermal fluctuations to
stabilize mixture droplets against collapse.

We solve iteratively the energy (III.1) up to second-order in $\tilde{n}$ and
$\tilde{m}$. Our results are compared with recent DMC data Cik1  and
theoretical analyses Petrov ; Hu . Figure 1 shows that our findings are in
excellent agreement with the DMC simulations Cik1 , and improve upon the
conventional Bogoliubov and the pairing predictions previously reported in
Petrov ; Hu  regardless o the value of $a_{12}$. This indicates the importance
of the beyond-LHY quantum corrections.

Figure 3:  Anomalous and noncondensed equilibrium density of the self-bound
droplet state from Eqs.(8) and (9) as a function of $a_{12}/a$.

The equilibrium density of the droplet can be obtained by minimizing the
energy with respect to the density Petrov . Using the identity
$(1+a_{12}/a)^{j/2}+(1-a_{12}/a)^{j/2}=2^{j/2}$, we immediately get

$\displaystyle\frac{n_{eq}}{n_{0}}$ $\displaystyle=\frac{0.0576}{\left(\delta
a/a\right)_{+}^{2}}\Big{[}18-25\left(\delta a/a\right)_{+}$ (7)
$\displaystyle-6\sqrt{9-25\left(\delta a/a\right)_{+}}\Big{]}.$

This equation differs by the term $[0.0576/\left(\delta
a_{+}/a\right)^{2}][18-25\left(\delta a_{+}/a\right)-6\sqrt{9-25\left(\delta
a_{+}/a\right)}]$ from the equilibrium density
$n_{0}a^{3}=\left(25\pi/16384\right)\left(\delta a_{+}/a\right)^{2}$ predicted
by Petrov Petrov . In Fig.2 we plot our predictions for the equilibrium
density and compare the results with the DMC simulation points Cik1  and the
existing previous theoretical models Petrov ; Hu . Figure 2.(a) shows that our
theoretical results based on the HFB theory are in good agreement with the DMC
simulation practically in the whole range of interactions. However, for large
values of $|a_{12}/a|$, the ratio $n_{eq}/n_{0}$ assessed by our theory
diverges from that of the Bogoliubov approach Petrov  and the pairing theory
Hu  [see Fig.2.(b)]. One might argue that quantum fluctuations are responsible
for such an increase in our equilibrium density.

Figure 4:  (a) Free energy as a function of the density for $a_{12}=-1.1a$ at
different values of temperature. (b) Critical temperature $T_{c}$ of a droplet
as a function of the inter-species interaction strength $a_{12}/a$. Dotted
lines are our results up to first-order in $\tilde{m}$ and $\tilde{n}$. Dashed
lines are our results up to the second-order in $\tilde{m}$ and $\tilde{n}$.

We now focus ourselves to calculate the noncondensate and anomalous densities
in the equilibrium phase. Inserting Eq.(7) into Eqs.(4) and (5), one can
rewrite the noncondensate and anomalous densities of the droplet in terms of
the equilibrium density up to first-order :

$\displaystyle\tilde{n}_{eq}=\frac{2\sqrt{2}}{3}n_{eq}\sqrt{\frac{a^{3}n_{eq}}{\pi}}\sum_{\pm}\left(\frac{\delta
a}{a}\right)_{\pm}^{3/2},$ (8)

and

$\displaystyle\tilde{m}_{eq}=2\sqrt{2}\,n_{eq}\sqrt{\frac{a^{3}n_{eq}}{\pi}}\sum_{\pm}\left(\frac{\delta
a}{a}\right)_{\pm}^{3/2}.$ (9)

To the best of our knowledge, these quantities have never been determined in
the literature. Remarkably, Eqs.(8) and (9) show that the anomalous density of
the droplet is larger than the noncondensed density specifically for large
$|a_{12}/a|$ as in the case of an ordinary BEC. We see from Fig.3 that
$\tilde{n}_{eq}$ and $\tilde{m}_{eq}$ increase with the ratio $|a_{12}/a|$
which is indeed natural since both quantities arise from interactions.

###  III.2 Self-bound droplet at finite temperature

At finite temperature, the free energy (II) can be rewritten as a function of
the small parameter of the theory $na^{3}\ll 1$ as:

$\displaystyle\frac{F}{NE_{0}}=$
$\displaystyle\frac{E}{NE_{0}}-\frac{8\sqrt{\pi}}{675}(na^{3})^{1/6}\bigg{(}1+\frac{\tilde{m}-\tilde{n}}{n}\bigg{)}^{-3/2}$
$\displaystyle\times\sum_{\pm}\left(\frac{\delta
a}{a}\right)_{\pm}^{-3/2}\left(\frac{T}{E_{0}}\right)^{4},$ (10)

where $E$ is given in Eq.(III.1).

In Fig.4.(a), we present the selfconsistent solution of the free energy
(III.2) up to second order for several values of temperature. We see that the
local minimum disappears at the critical temperature $T_{c}\simeq 1.3\times
10^{-4}E_{0}$ at which the thermal fluctuations compensate the repulsive
quantum fluctuations. This value is greater by $\sim 23\%$ than the threshold
temperature obtained in Refs.Hu ; Ota  owing to the higher-order beyond-LHY
corrections. At higher temperatures ($T\simeq 1.5\times 10^{-4}E_{0}$), the
self-bound droplet is completely evaporated due to the strong thermal
fluctuations. Another important remark is that the higher-order correlations
may lead to a shift in the critical temperature.

A useful analytical expression for the critical temperature above which the
droplet destabilizes can be extracted by minimizing the free energy (III.2)
with respect to the density. This yields

$\displaystyle\frac{T_{c}}{E_{0}}=0.0114\left(\frac{\delta
a}{a}\right)_{+}^{2}\bigg{[}72+25\left(\frac{\delta
a}{a}\right)_{+}\bigg{]}^{1/4}.$ (11)

In Fig.4.(b) we plot the critical temperature from Eq.(11) and compare it with
the results of Ref.Hu . We observe that $T_{c}$ decreases with decreasing the
ratio $|a_{12}/a|$. Our findings agree with the pairing predictions only for
small $|a_{12}/a|$ and then both theories diverge from each other on account
of the higher-order contribution of the anomalous correlations.

Figure 5:  Thermal equilibrium density as a function of the temperature for
different values of interspecies interaction strength.

Treating the temperature term as a perturbation Ayb , one gets for the
equilibrium density:

$\displaystyle n_{eq}^{T}\approx n_{eq}\left[1+\frac{5(\delta
a/a)_{+}^{-4}}{128\sqrt{3}}\Big{(}\frac{T/E_{0}}{a^{3}n_{eq}}\Big{)}^{4}\right],$
(12)

and its behavior is displayed in Fig.5. We see that $n_{eq}^{T}$ is increasing
with temperature and interspecies interactions.

##  IV Ground-state properties

In this section we examine the ground-state properties of the droplet by
discussing the effects of higher-order quantum and thermal corrections in the
different ground-state solutions such as the density profiles and the
collective modes.

###  IV.1 Density profiles

Following the method outlined by Petrov Petrov , we rescale the wavefunction
of each species as: $\Phi(r,t)=\sqrt{n_{0}}\phi(r,t)$, where $\phi(r,t)$ is a
scalar wavefunction common to both species. It is convenient to introduce the
dimensionless parameters : $\tilde{N}=N/(n_{0}\xi^{3})$, $\tilde{r}=r/\xi$,
and $\tilde{t}=t/\tau$, where $\xi=\sqrt{6\hbar^{2}/\left(m|\delta
g|n_{0}\right)}$, and $\tau=6\hbar/\left(|\delta g|n_{0}\right)$.

Figure 6:  (a) Density profile of the self-bound droplet obtained from
Eq.(IV.1) at zero temperature for different values of interaction strength
$\delta g/g$, and $\tilde{N}=3000$. Solid lines: Our results. Dashed lines:
Petrov’s results Petrov . Dotted lines: Findings of Ref.Ota . (b) Density
profile of the self-bound droplet at different values of temperature for
$g_{12}/g=-1.1$.

Upon inserting expressions of $\tilde{m}$ and $\tilde{n}$ up to the second-
order into Eq.(1), and using our variational TDHFB equation $id\phi/dt=d{\cal
E}/d\phi$ Boudj1 ; Boudj2 , we get the generalized finite-temperature GPE

$\displaystyle
i\frac{d\phi}{d\tilde{t}}=\bigg{\\{}-\frac{1}{2}\Delta_{\bf\tilde{r}}-3|\phi|^{2}+\frac{5}{2}[1+\alpha|\phi|(1+\alpha|\phi|)^{3/2}]^{5/2}|\phi|^{3}$
$\displaystyle+\frac{5}{2}\alpha^{\text{T}}|\phi|^{-1}\bigg{\\}}\phi,$ (13)

where $\alpha=(4\sqrt{2}/3)(\delta g/g)_{+}^{3/2}\sqrt{a^{3}n_{0}/\pi}$, and
$\alpha^{\text{T}}=-(4/9)\sqrt{a^{3}n_{0}/\pi}(\delta
g/g)_{+}^{-1/2}\big{(}T/4a^{3}n_{0}E_{0}\big{)}^{2}$. For
$\alpha^{\text{T}}=0$, Eq.(IV.1) reduces to that used in Ref.Ota . For
$\alpha=\alpha^{\text{T}}=0$, one recovers the Petrov’s generalized GPE Petrov
. Equation (IV.1) is appealing since it obviously contains higher-order
quantum corrections and temperature effects. It is valid when the spatial
variations of the condensate density are small over the extended healing
length, $\xi$, validating a description of a locally homogeneous system.

We solve numerically Eq.(IV.1) and the corresponding results are captured in
Fig. 6.

It is clearly seen in Fig. 6.(a) that as the interaction strength $|\delta
g/g|$ increases, the central density decreases and the size of the droplet
increases in agreement with the findings of Ota . This may provoke a shift in
the critical number compared to the universal results of Ref.Petrov .

Figure 6 (b) shows that the droplet exhibits a flat-top profile at the
equilibrium (i.e. at $T<T_{c}$). When the temperature augments, a dramatic
change in the shape of the spatial density profile is observed where the
condensate is a Gaussian-like shape at $T\geq T_{c}$. This signals the
evaporation of the droplet to the gas state as is foreseen above.

###  IV.2 Collective modes

The collective modes of the self-bound droplet can be analytically explored
employing a Gaussian ansatz

$\displaystyle\phi(r)=A\exp\bigg{[}-\left(\frac{1}{2\sigma^{2}}+i\beta\right)r^{2}\bigg{]},$
(14)

where $A=\sqrt{\tilde{N}/\left(\pi^{3/2}\sigma^{3}\right)}$ is the
normalization factor, $\sigma$, and $\beta$ are the variational parameters
standing for the width and the phase of the condensate, respectively.

The dynamics of the droplet width can be obtained from the Euler-Lagrange
equations

$\ddot{\sigma}=-\frac{dV_{eff}(\sigma)}{d\sigma},$ (15)

where the effective potential is given by:

$\displaystyle V_{eff}(\sigma)$
$\displaystyle=-\frac{\tilde{N}}{2\sigma^{2}}+\frac{\sqrt{2}}{3\pi^{3/2}}\frac{\tilde{N}}{\sigma^{3}}+\frac{8\sqrt{2/5}}{27\pi^{9/4}}\frac{\tilde{N}^{3/2}}{\sigma^{9/2}}$
(16)
$\displaystyle+\frac{4\sqrt{2/5}}{18\pi^{3}}\frac{\tilde{N}}{\sigma^{6}}\alpha+\frac{16\sqrt{2/\tilde{N}}\pi^{3/4}\sigma}{9}\alpha^{T}.$

Frequencies of the low-lying collective modes are determined by linearizing
Eq.(16) around the equilibrium solutions. The breathing modes can be found by
inserting the decomposition $\sigma(t)=\sigma+\delta\sigma(t)$ (with
$\delta\sigma(t)=\delta\sigma\,e^{i\omega t}$ and $\delta\sigma\ll\sigma$)
into Eq.(16) and expanding the effective potential into a Taylor series. As
shown in Fig.7, both temperature and interactions $|\delta g/g|$ tend to
increase the frequency of the breathing modes. Another important remark is
that higher-order effects in the breathing modes are substantial.

Figure 7: Frequencies of low-lying collective excitation modes as a function
of the total atom number in units of $\tilde{N}$. Dotted line: $\delta
g/g=-0.05$, dashed line: $\delta g/g=-0.1$, solid line: $\delta g/g=-0.5$.
Here $N_{c}=18.65$ is the critical number found in Ref.Petrov .

##  V Conclusion and outlook

In this work we investigated the effects of higher-order quantum and thermal
fluctuations on the properties of a symmetric mixture self-bound droplet using
the beyond-LHY theory. We derived useful formulas for the energy, the
equilibrium, the normal and anomalous densities, and the critical temperature
of the self-bound droplet. Our predictions have been tested against DMC data
and the existing theoretical results, and excellent agreement is found.
Finite-temperature effects on the ground-state properties of the self-bound
droplet have been also addressed by numerically solving the extended finite-
temperature GPE. We showed that the inclusion of higher-order quantum and
thermal corrections may shift the density profiles and frequencies of the
breathing oscillations. We conjecture that our model can be extended to the
asymmetric droplet for a direct comparison with recent experiments. This is
one of our future works.

##  VI Acknowledgments

We are grateful to Jordi Boronat and Viktor Cikojević for giving us the DMC
data. We acknowledge support from the Algerian ministry of higher education
and scientific research under Research Grant No. PRFU-B00L02UN020120190001.

## References

  * (1) D. S. Petrov, Phys. Rev. Lett. 115, 155302 (2015).
  * (2) A. Cappellaro, T. Macrí, G. F. Bertacco and L. Salasnich, Sci. Rep. 7, 13358 (2017).
  * (3) V. Cikojević, K. Dzelalija, P. Stipanović, L. Vranješ Markić and J. Boronat, Phys. Rev. B 97, 140502 (R) (2018).
  * (4) C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, L. Tarruell, Science 359, 301 (2018); P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor, L. Tanzi, and L. Tarruell, Phys. Rev. Lett. 120, 135301 (2018).
  * (5) G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F.Minardi, M. Modugno, G. Modugno, M. Inguscio, and M. Fattori, Phys. Rev. Lett. 120, 235301 (2018).
  * (6) C. D’Errico, A. Burchianti, M. Prevedelli, L. Salasnich, F. Ancilotto, M. Modugno, F. Minardi, and C. Fort, Phys. Rev. Research 1, 033155 (2019).
  * (7) A. Boudjemâa, Phys. Rev. A 98, 033612 (2018).
  * (8) C. Staudinger, F. Mazzanti, and R. E. Zillich, Phys. Rev. A 98, 023633 (2018).
  * (9) D. S. Petrov and G. E. Astrakharchik, Phys. Rev. Lett. 117, 100401 (2016).
  * (10) E. Chiquillo, Phys. Rev. A 99, 051601(R) (2019).
  * (11) G. E. Astrakharchik, B. A. Malomed, Phys. Rev. A 98, 013631 (2018).
  * (12) M. Tylutki, G. E. Astrakharchik, . A. Malomed, and D. S. Petrov, Phys. Rev. A 101, 051601(R) (2020).
  * (13) I. Morera, G. E. Astrakharchik, A. Polls, B. Juliá-Díaz, Phys. Rev. Research 2, 022008 (R) (2020).
  * (14) R. Sachdeva, M. N. Tengstrand, S. M. Reimann, 2001.04751Phys. Rev. A 102, 043304 (2020).
  * (15) A. Boudjemâa and M. Benarous, Eur. Phys. J. D 59, 427 (2010) .
  * (16) A. Boudjemâa and M. Benarous, Phys. Rev. A 84, 043633 (2011).
  * (17) A. Boudjemâa, Phys. Rev. A 86, 043608 (2012).
  * (18) A. Boudjemâa, Phys. Rev. A 88, 023619 (2013).
  * (19) A. Boudjemâa, Phys. Rev. A 90, 013628 (2014).
  * (20) A. Boudjemâa, J. Phys. A: Math. Theor. 48 045002 (2015).
  * (21) A. Boudjemâa, Phys. Rev. A 91, 063633 (2015).
  * (22) A. Boudjemâa, Commun. Nonlinear Sci. Numer. Simul. 33, 85 (2016).
  * (23) A. Boudjemâa, Commun. Nonlinear Sci. Numer. Simul. 48, 376 (2017).
  * (24) A. Boudjemâa, Phys. Rev. A 94, 053629 (2016).
  * (25) A. Boudjemâa, and N. Guebli, J. Phys. A: Math. Theor. 50, 425004 (2017).
  * (26) A. Boudjemâa, and N. Guebli, Phys. Rev. A 102, 023302 (2020).
  * (27) A. Boudjemâa, Phys. Rev. A 97, 033627 (2018).
  * (28) V. Cikojevíc, L. Vranjes Markic, G. E. Astrakharchik, and J. Boronat, Phys. Rev. A 99, 023618 (2019).
  * (29) M. Ota, and G. E. Astrakharchik, SciPost Phys. 9, 020 (2020).
  * (30) H. Hu and X-J. Liu, Phys. Rev. Lett. 125, 195302 (2020).
  * (31) A. Boudjemâa, Degenerate Bose Gas at Finite Temperatures, (Lambert Academic Publishing, Saarbrücken, Germany, 2017).
  * (32) J. O. Andersen, Rev. Mod. Phys 76, 599 (2004).
  * (33) A. Griffin and H. Shi, Phys. Rep. 304, 1 (1998).
  * (34) R.A. Duine and H.T.C. Stoof, J. Opt. B: Quantum Semiclass. Opt. 5, S212 (2003).
  * (35) E. Aybar and M. Ö. Oktel, Phys. Rev. A 99, 013620 (2019).
  * (36) Q. Gu and L. Yin, Phys. Rev. B 102, 220503 (R) (2020).
  * (37) A. Boudjemâa, Ann. Phys. 381, 68 (2017).


# Billion-years old proteins show the importance of N-lobe orientation in
Imatinib-kinase selectivity.

Zahra Shamsi
Department of Chemical and Biomolecular Engineering
University of Illinois at Urbana-Champaign
Urbana
IL
Diwakar Shukla
Department of Chemical and Biomolecular Engineering
University of Illinois at Urbana-Champaign
Urbana
IL
<EMAIL_ADDRESS>

###### Abstract

The molecular origins of proteins’ functions are a combinatorial search
problem in the proteins’ sequence space, which requires enormous resources to
solve. However, evolution has already solved this optimization problem for us,
leaving behind suboptimal solutions along the way. Comparing suboptimal
proteins along the evolutionary pathway, or ancestors, with more optimal
modern proteins can lead us to the exact molecular origins of a particular
function. In this paper, we study the long-standing question of the
selectivity of Imatinib, an anti-cancer kinase inhibitor drug. We study two
related kinases, Src and Abl, and four of their common ancestors, to which
Imatinib has significantly different affinities. Our results show that the
orientation of the N-lobe with respect to the C-lobe varies between the
kinases along their evolutionary pathway and is consistent with Imatinib’s
inhibition constants as measured experimentally. The conformation of the DFG-
motif (Asp-Phe-Gly) and the structure of the P-loop also seem to have
different stable conformations along the evolutionary pathway, which is
aligned with Imatinib’s affinity.

_K_ eywords Evolutionary pathway $\cdot$ Kinase inhibitor $\cdot$ Simulations

##  1 Introduction

Protein kinases are a class of enzyme that transfer phosphate groups from ATP
to other proteins, thereby signaling growth and cell proliferation. Mutations
in kinases can lead to uncontrolled cell growth and eventually cancer, making
kinases a prime target for drug design, typically of small molecule inhibitors
[1, 2]. Imatinib is one of the clinically successful drugs for the treatment
of multiple cancers like chronic myelogenous leukemia [3]. It selectively
inhibits Abl and not other structurally similar kinases like Src. Since the
overactive Abl mutant only exists in cancer cells, Imatinib has a limited
effect on healthy cells [4]. Why does the drug inhibit Abl and not Src despite
the high protein sequence identity between Src and Abl ($\sim$46%) [4]? After
two decades of study on the basis of Imatinib selectivity between Abl and Src,
the question still remains unanswered. During these years, extensive work has
been done to elucidate this problem using different approaches from NMR and
fast kinetics [5], sequence swapping [4], and ancestral gene reconstruction
[6] experiments to free energy calculations and long time-scale molecular
dynamics (MD) simulations [7, 8, 9]. Each of these studies gave insight into
different aspects of the problem, but a full answer is still missing.

DFG-motif (Asp-Phe-Gly) is a highly conserved segment of the activation loop
in kinase domains that is proposed to play a major role in the selection
mechanism since Imatinib only binds a specific configuration of DFG. Two
conformations of the DFG motif are the inactive “DFG-out” and the active “DFG-
in”. Imatinib only binds to DFG-out conformation of kinases. Multiple groups
have argued that the kinetic basis of Imatinib selectivity for Abl, compared
to Src kinase is because of a pre-existing equilibrium between DFG-out and
DFG-in, or so-called a conformational selection mechanism [10, 7]. More
recently despite the general belief on the critical role of conformational
selection, Agafonov $et$ $al.$ claimed that the Imatinib selectivity is rooted
in conformational changes after drug binding, not before [5]. They support
their hypothesis using NMR studies and showed the presence of induced-fit
mechanism in Abl kinase. The next step was to find a sequence-function
relationship, and clarify which set of residues is responsible for the
accessibility of induced-fit mechanism in Abl, but not in Src. Sequence
swapping had been performed to make Src similar to Abl but all studies had
failed to illuminate the atomistic determinants of selectivity [4, 7]. On the
other hand, phylogenetic studies have been used as powerful tools to study
protein sequence-function relationships in different problems [11, 12, 13].
Therefore, following the NMR study, Wilson $et$ $al.$ recreated the
evolutionary pathway between Src and Abl by resurrecting the common ancestors
between them [6]. Using x-ray structures of an ancestral kinase and binding
kinetics data from NMR, they showed the evolution of their proposed induced-
fit mechanism and its effect on the selectivity of Imatinib [6].

However, proteins are not just sequences of amino acids; they can adopt
thousands of different 3D conformations, among which only specific sets of
conformations are functional. In the aforementioned studies with the plausible
sequence-function relationship, still, the sequence-structure-function
relationship is missing. Even though these studies have shown the presence of
induced-fit mechanism, they do not give any details on what are these
mechanisms in the protein’s structure. In this study, we investigate the
sequence-structure relationship by reconstructing the common ancestors of Src
and Abl computationally.

Study of evolutionary pathways is a natural way to identify the key amino acid
changes that differentiate one family member from another [14]. Evolution
generates functional proteins in every stage despite the large sequence
differences between them. It diversifies functions by altering their structure
and the associated free-energy landscapes. The differences between Abl and Src
also have evolved over a billion years from their common ancestor. Here, we
borrow the sequences of the kinase ancestors from literature [6] and
reconstruct them computationally (Figure 1(a) which is adapted from ref.
[15]). To elucidate the sequence-structure-function relationship, we simulated
four ancestors, ANC-A1, ANC-A2, ANC-AS, and ANC-S1, plus ANC-AS with 15
suggested mutations which changed () the affinity significantly [6]. We refer
to ANC-AS+15 as ANC-AS(+15) in this paper to prevent any confusion. We compare
the results of these simulations with the experimental values of free energies
and inhibition constants of the ancestors. We also denote the evolution of
conformations that play role in the conformational selection mechanism and
induced-fit mechanism.

Figure 1: (a) Phylogenetic tree of Abl and Src families showing the
reconstructed nodes.(b) Crystal structure of Abl kinase (PDB ID: 2HYY [16]) is
shown with labels for important regions.

##  2 Results

###  2.1 Evolution of conformational differences: N-lobe rotation

Comparison of crystal structures shows a significant difference in the
orientation of N-lobe with respect to the C-lobe between Src and Abl. This
rotational angle is known to play a role in kinase activation, as the
catalytic cleft between the lobes is relatively closed in inactive kinases
[6]. The ancestral study by Wilson $et$ $al.$ suggested that altering only 15
residues in ANC-AS’s N-lobe to the corresponding residues in Abl kinase
drastically increased the Imatinib affinity to a level similar to Abl. This
indicates the importance of N-lobe residues in the drug binding mechanism [6].
Therefore, it is likely that lobe rotation can be the dynamics effect of this
15 residue mutations in N-lobe. Here, we want to calculate the orientation of
N-lobe along the evolutionary pathway and see if it is correlated with the
Imatinib binding affinities.

We studied the dynamics of four common ancestors of Src and Abl, an ancestor
with the 15 mutations, and the two modern kinases. In these long timescale
unbiased MD simulations no Imatinib molecule was present. The simulations were
performed using adaptive sampling technique and Markov State Model (MSM)
analysis [17, 18, 19]. In total, we performed 0.545 ms of aggregated unbiased
MD simulation for the seven apo kinases. In order to quantify the lobe
rotation, we define vectors in C-lobe (V1) and N-lobe (V2), and calculate the
angle between them (named as $\theta$) as shown in Figure 2(a). Our results
show a gradual shift in the density distributions of $\theta$ as we move from
Abl to Src kinase on the phylogenetic tree (Figure 2 (b)). The peak values of
$\theta$ distributions are well correlated with the experimental values of
Imatinib inhibition constants ($K_{i}$) at 25°C as reported in the literature
[6]. This result suggests that the available area at the N-lobe and C-lobe is
one of the major factors affecting the Imatinib binding affinity.

Figure 2:  N-lobe rotational angle of different kinases. (a) The vectors
defining the rotational angle are shown on the kinase structure. (b) The
distribution of angles for different kinases reveals how position of N-lobe
with respect to C-lobe evolved along the evolutionary pathway and the angle
shifted from Src to Abl. (c) The correlation between the most likely value of
$\theta$ calculated from the simulation and the experimental values of
Imatinib inhibition constants ($K_{i}$) at 25°C as reported in literature [6].

###  2.2 Evolution of induced-fit mechanism: secondary structure of P-loop.

Changes in the secondary structure of P-loop are critical reaction coordinates
for the Imatinib binding process. To observe the evolution of P-loop
conformation, we calculated the root mean square deviation (RMSD) of the
P-loop from the crystal structures of Abl and Src in the kinases (Figure 3).
Based on our observations, Src-like P-loop is stable in all seven kinases,
whereas the ability to form Abl-like P-loop has been lost in ANC-S1 and Src.
This suggests Imatinib binding from the P-loop side and the following induced-
fit mechanism is not possible in Src and ANC-S1, but it could be feasible in
the other five kinases. The exact residues responsible for Abl-like kinked
P-loop should be a subset of the residue differences between ANC-AS and ANC-S1
(ANC-AS and ANC-S1 are 82% identical as shown in Figure S1).

Looking at the simulations of Abl kinase and its crystal structure, we can
observe that two inter residue interactions responsible for Abl P-loop kinked
conformation are Y253-N322 and Q252-N322. Both of these interactions are
conserved in ANC-A1, ANC-A2, and ANC-AS(+15), which makes P-loop more likely
to form a helical structure as shown in Figure 3(a). However, only one of the
interactions is conserved in ANC-AS and ANC-S1, and none in Src. This
partially explains the lower helical content of P-loop in ANC-AS, ANC-S1, and
Src. The corresponding residue pairs are F-S and Q-S in ANC-AS/ANC-S1 and F-S
and C-S in Src kinase.

Figure 3: Evolutionary pathway of P-loop’s conformation. Secondary structure
of P-loop in (a) Abl [16] and (b) Src’s [4] crystal structure is shown.
Configuration of P-loop is the most distinct difference between Abl and Src’s
crystal structures. (c) Free energy landscapes of P-loop conformation show how
ancestors lose their Abl-like P-loop structure as they get closer to Src. (d)
Src and ancestors closer to Src and less likely to forms helical P-loop.
Residue numbering is based on the sequence of ANC-AS(+15) as presented in the
Figure S2. Colorbar shows free energies in kcal/mol.

###  2.3 Evolution of conformational selection: DFG flip’s mechanism.

To investigate the dynamics of the DFG-motif, distances between two pairs of
residues were measured as shown in Figure 4. Our results indicate that all
ancestors are able to adopt both DFG-in and DFG-out conformations. In Src
kinase, DFG-out conformation is relatively less stable, whereas it is stable
in Abl and all the ancestors. The difference in the stability of DFG-out
conformation between Src and other kinases is $\sim$1 kcal/mol, which is not
sufficient to justify the $\sim$3000 fold difference in their inhibition
constants [4]. Stable DFG-out conformation in the ancestors and Abl, indicates
that the residues responsible for less stable DFG-out conformation in Src are
among the set of residue differences between ANC-S1 and Src (ANC-S1 and Src
are 75% identical as shown in Figure S1).

Figure 4:  Stability of DFG-in and DFG-out conformations. (a) Free energy
landscapes of DFG flip for Abl, ancestors and Src are shown. When both the
distances are greater than 1 nm, DFG motif is in DFG-out conformation, and
when both of them are less than 1 nm it is in DFG-in conformation. Imatinib
only binds DFG-out conformation. Colorbar shows free energies in kcal/mol.
DFG-in and DFG-out conformations is shown in (b). Residue numbering is based
on the sequence of ANC-AS(+15) as presented in the SI Figure S2

###  2.4 Evolution of conformational differences: secondary structure of
A-loop.

The activation-loop (A-loop) in a kinase is usually includes the highly
conserved DFG motif and 17 of its following residues. This region adopts a
closed or inactive conformation and an open or active conformation. In open
conformations of A-loop, substrate proteins are able to bind the kinases for
phosphate transfer, whereas in closed configurations, ATP is not exposed to
the substrate proteins and phosphate transfer is not feasible. Imatinib binds
the closed conformation of Abl and is incapable of binding to its open
configuration [20].

Analyzing the simulations of the seven apo kinases showed that the Abl-like
inactive configuration of A-loop is accessible in Abl’s closer ancestors and
disappears in ANC-S1. Therefore residues responsible for inactive Abl-like
A-loop conformation should be among the sequence differences between ANC-S1
and ANC-AS (Figure 5). On the other hand, Src-like inactive structure of
A-loop appears in ANC-AS, which shows residues coding for Src-like helical
conformation of A-loop are among the sequence differences between ANC-AS and
ANC-AS(+15) (Figure 5). We also see in Figure 5(d) that A-loop in Src and
ancestors closer to it tend to be less structured as compared to Abl.

Figure 5: Evolutionary pathway of A-loop’s conformation. Secondary structure
of A-loop in (a) Abl [16] and (b) Src [4] crystal structure is shown. (c) Free
energy landscapes of A-loop show how ancestors gain the ability to form Src-
like A-loop structure as it gets closer to Src. Colorbar shows free energies
in kcal/mol. (d) Src and ancestors closer to Src and less likely to
$\beta$sheet A-loop. Residue numbering is based on the sequence of ANC-
AS(+15).

##  3 Discussion

Protein-ligand selectivity remains a mysterious phenomenon in biology due to
the lack of knowledge in biophysical mechanism of ligand binding. Here, we
studied the long-lasting problem of Imatinib selectivity toward Abl kinase,
using long MD simulations of Src, Abl, and their common ancestors. We compared
the simulation results with experimental activity and binding free energy
values to better understand the mechanism. We evaluated two of the best-known
hypotheses that try to explain the selectivity mechanism by following their
corresponding conformational changes along the evolutionary pathway
connections Abl and Src kinases. The first hypothesis assumes the high
selectivity of Imatinib towards Abl is due to its particular conformation of
DFG motif. We compared the dynamics of DFG motif in Src, Alb and their
ancestors and observed that all of them can adopt DFG-out conformation, which
is crucial for Imatinib binding. However, DFG-out conformation seemed to be
less stable in Src kinase as compared to the other six. This suggests that
Imatinib can potentially bind to all of these kinases, but it is less likely
to find the desirable conformation in Src kinase and bind to it. The second
hypothesis that we studied was induced-fit mechanism. Even though we still do
not know what exactly the induced-fit mechanism is, we believe P-loop
conformation plays a role in the process. Therefore we studied the evolution
of 3D structure of P-loop and observed that the helical P-loop adapts a
$\beta$-sheet-rich conformation as we go from Src to Abl kinase on their
phylogenetic tree. This change seemed to be very gradual along the
evolutionary pathway and align with the Imatinib affinity towards the kinases
measured by experiments [6].

Although Imatinib has achieved remarkable success in treating chronic myeloid
leukemia, the emergence of resistance to this agent weakens the prospect of a
cure for this leukemia. Imatinib resistance in most of the patients coincides
with single point mutations in Abl kinase. In this study, we also rationalized
the origin of some of these mutations.

A set of mutations that reported to be among the most frequent mutation in
Imatinib resistance patients are Y253F, Y253H, and Q252H in Abl kinase [21].
We showed that these two residues are crucial to forming the helical
conformation in the P-loop by making Y253-N322 and Q252-N322 interactions.

In conclusion, our data suggest that future drug design efforts should focus
more on understanding the binding pathways and their corresponding induced-fit
steps. To this end, the full binding process of more drugs and drug targets
need to be studied at atomic resolution with the hope that new insights will
translate into improved next-generation compounds. As we learn more about
Imatinib selectivity mechanism, these insights could be transferred to
understand other selectivity challenges.

## References

  * [1] Zahra Shamsi and Diwakar Shukla.  How does evolution design functional free energy landscapes of proteins? a case study on the emergence of regulation in the cyclin dependent kinase family.  Molecular Systems Design & Engineering, 5(1):392–400, 2020.
  * [2] Charles L Sawyers.  Disabling abl—perspectives on abl kinase regulation and cancer therapeutics.  Cancer Cell, 1(1):13–15, February 2002.
  * [3] Renaud Capdeville, Elisabeth Buchdunger, Juerg Zimmermann, and Alex Matter.  Glivec (STI571, imatinib), a rationally developed, targeted anticancer drug.  Nature Reviews Drug Discovery, 1(7):493–502, July 2002.
  * [4] Markus A. Seeliger, Bhushan Nagar, Filipp Frank, Xiaoxian Cao, M. Nidanie Henderson, and John Kuriyan.  c-src binds to the cancer drug imatinib with an inactive abl/c-kit conformation and a distributed thermodynamic penalty.  Structure, 15(3):299–311, mar 2007.
  * [5] Roman V Agafonov, Christopher Wilson, Renee Otten, Vanessa Buosi, and Dorothee Kern.  Energetic dissection of gleevec’s selectivity toward human tyrosine kinases.  Nature Structural & Molecular Biology, 21(10):848–853, September 2014.
  * [6] C. Wilson, R. V. Agafonov, M. Hoemberger, S. Kutter, A. Zorba, J. Halpin, V. Buosi, R. Otten, D. Waterman, D. L. Theobald, and D. Kern.  Using ancient protein kinases to unravel a modern cancer drug's mechanism.  Science, 347(6224):882–886, February 2015.
  * [7] Y.-L. Lin, Y. Meng, W. Jiang, and B. Roux.  Explaining why gleevec is a specific and potent inhibitor of abl kinase.  Proceedings of the National Academy of Sciences, 110(5):1664–1669, January 2013.
  * [8] Yen-Lin Lin and Benoît Roux.  Computational analysis of the binding specificity of gleevec to abl, c-kit, lck, and c-src tyrosine kinases.  Journal of the American Chemical Society, 135(39):14741–14753, September 2013.
  * [9] Lulu Wang, Guodong Zheng, Xianxian Liu, Duan Ni, Xinheng He, Jinying Cheng, and Shaoyong Lu.  Molecular dynamics simulations provide insights into the origin of gleevec’s selectivity toward human tyrosine kinases.  Journal of Biomolecular Structure and Dynamics, 37(10):2733–2744, November 2018.
  * [10] Y. Shan, M. A. Seeliger, M. P. Eastwood, F. Frank, H. Xu, M. O Jensen, R. O. Dror, J. Kuriyan, and D. E. Shaw.  A conserved protonation-dependent switch controls drug binding in the abl kinase.  Proceedings of the National Academy of Sciences, 106(1):139–144, December 2008.
  * [11] Michael J. Harms and Joseph W. Thornton.  Evolutionary biochemistry: revealing the historical and physical causes of protein properties.  Nature Reviews Genetics, 14(8):559–571, July 2013.
  * [12] Yosephin Gumulya, Jong-Min Baek, Shun-Jie Wun, Raine E. S. Thomson, Kurt L. Harris, Dominic J. B. Hunter, James B. Y. H. Behrendorff, Justyna Kulig, Shan Zheng, Xueming Wu, Bin Wu, Jeanette E. Stok, James J. De Voss, Gerhard Schenk, Ulrik Jurva, Shalini Andersson, Emre M. Isin, Mikael Bodén, Luke Guddat, and Elizabeth M. J. Gillam.  Engineering highly functional thermostable proteins using ancestral sequence reconstruction.  Nature Catalysis, 1(11):878–888, October 2018.
  * [13] Jordi Paps and Peter W. H. Holland.  Reconstruction of the ancestral metazoan genome reveals an increase in genomic novelty.  Nature Communications, 9(1), April 2018.
  * [14] Michael J Harms and Joseph W Thornton.  Analyzing protein structure and function using ancestral gene reconstruction.  Current Opinion in Structural Biology, 20(3):360–366, June 2010\.
  * [15] Roman V. Agafonov, Christopher Wilson, and Dorothee Kern.  Evolution and intelligent design in drug development.  Frontiers in Molecular Biosciences, 2, May 2015.
  * [16] Sandra W. Cowan-Jacob, Gabriele Fendrich, Andreas Floersheimer, Pascal Furet, Janis Liebetanz, Gabriele Rummel, Paul Rheinberger, Mario Centeleghe, Doriano Fabbro, and Paul W. Manley.  Structural biology contributions to the discovery of drugs to treat chronic myelogenous leukaemia.  Acta Crystallographica Section D Biological Crystallography, 63(1):80–93, dec 2006.
  * [17] Zahra Shamsi, Alexander Moffett, and Diwakar Shukla.  Enhanced unbiased sampling of protein dynamics using evolutionary coupling information.  Sci. Rep., 5, sept 2017.  in-press.
  * [18] Saurabh Shukla, Zahra Shamsi, Alexander S. Moffett, Balaji Selvam, and Diwakar Shukla.  Application of hidden markov models in biomolecular simulations.  In Hidden Markov Models, pages 29–41. Springer New York, 2017.
  * [19] Zahra Shamsi, Kevin J Cheng, and Diwakar Shukla.  Reinforcement learning based adaptive sampling: Reaping rewards by exploring protein conformational landscapes.  The Journal of Physical Chemistry B, 122(35):8386–8395, 2018.
  * [20] Paul La Rosée, Amie S Corbin, Eric P Stoffregen, Michael W Deininger, and Brian J Druker.  Activity of the bcr-abl kinase inhibitor pd180970 against clinically relevant bcr-abl isoforms that cause resistance to imatinib mesylate (gleevec, sti571).  Cancer research, 62(24):7149–7153, 2002.
  * [21] Shundong Cang and Delong Liu.  P-loop mutations and novel therapeutic approaches for imatinib failures in chronic myeloid leukemia.  Journal of Hematology & Oncology, 1(1), October 2008.
  * [22] Szymon Chojnacki, Andrew Cowley, Joon Lee, Anna Foix, and Rodrigo Lopez.  Programmatic access to bioinformatics tools from EMBL-EBI update: 2017\.  Nucleic Acids Research, 45(W1):W550–W553, apr 2017.
  * [23] Xavier Robert and Patrice Gouet.  Deciphering key features in protein structures with the new ENDscript server.  Nucleic Acids Research, 42(W1):W320–W324, apr 2014.

Supplemental Materials

|  $K_{i}$ (nM) |  $\Delta G$ (kcal/mol)
---|---|---
Abl | 13 $\pm$ 10 | -4.5
ANC-A2 | 10 $\pm$ 4 | -5
ANC-AS(+15) | 21 $\pm$ 12 | -
ANC-A1 | 190 $\pm$ 65 | -3
ANC-AS | 910 $\pm$ 120 | -2.5
ANC-S1 | 8000 $\pm$ 1000 | -0.5
Src | 2000 $\pm$ 3000 | -0.5
Table S1:  Experimental values of Imatinib inhibition constants ($K_{i}$) at
25°C and $\Delta G$ associated with the induced-fit step from literature [6].

Figure S1: Multiple sequence alignment of Src, Abl, and reconstructed
ancestral sequences. Identical and similar residues are boxed in red and
yellow, respectively. Sequence alignments were created in T-Coffee [22] and
plotted by ESPRIPT [23]. (a) ANC-AS and ANC-AS(+15) are 96% identical. (b)
ANC-AS and ANC-S1 are 82% identical. (c) ANC-S1 and Src are 75% identical.

Figure S2: Multiple sequence alignment of Src, Abl, and reconstructed
ancestral sequences. Identical and similar residues are boxed in red and
yellow, respectively. Sequence alignments were created in T-Coffee [22] and
plotted by ESPRIPT [23].

| N-lobe pocket interactions | A-loop pocket interactions | P-loop kinked Conf.
---|---|---|---
| Y326L | F278Y | S346N
Abl | L | Y | N
ANC-A1 | L | Y | N
ANC-A2 | L | Y | N
AS(+15) | L | Y | N
ANC-AS | Y | F | S
ANC-S1 | Y | F | S
Src | Y | F | S
Table S2:  Conservation of some functionally important residues.

Figure S3: The angle between two lobs as a function of Imatinib binding.
Imatinib binds the kinases when $\theta$ is equal or greater than 120°.
Interaction with Imatinib increases the average value of the angles between
the lobes.


# High-Quality Supersampling via Mask-reinforced Deep Learning for Real-time
Rendering

Hongliang Yuan1, Boyu Zhang1,2, Mingyan Zhu1,3, Ligang Liu4, Jue Wang1
1Tencent AI Lab, 2Southeast University, 3Tsinghua University,
4University of Science and Technology of China
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(2018)

###### Abstract.

To generate high quality rendering images for real time applications, it is
often to trace only a few samples-per-pixel (spp) at a lower resolution and
then supersample to the high resolution. Based on the observation that the
rendered pixels at a low resolution are typically highly aliased, we present a
novel method for neural supersampling based on ray tracing 1/4-spp samples at
the high resolution. Our key insight is that the ray-traced samples at the
target resolution are accurate and reliable, which makes the supersampling an
interpolation problem. We present a mask-reinforced neural network to
reconstruct and interpolate high-quality image sequences. First, a novel
temporal accumulation network is introduced to compute the correlation between
current and previous features to significantly improve their temporal
stability. Then a reconstruct network based on a multi-scale U-Net with skip
connections is adopted for reconstruction and generation of the desired high-
resolution image. Experimental results and comparisons have shown that our
proposed method can generate higher quality results of supersampling, without
increasing the total number of ray-tracing samples, over current state-of-the-
art methods.

Monte Carlo denoising, neural networks, path tracing

††submissionid: 467††copyright: acmcopyright††journalyear: 2018††doi:
XXXXXXX.XXXXXXX††conference: Make sure to enter the correct conference title
from your rights confirmation emai; June 03–05, 2018; Woodstock,
NY††booktitle: Woodstock ’18: ACM Symposium on Neural Gaze Detection, June
03–05, 2018, Woodstock, NY††price: 15.00††isbn: 978-1-4503-XXXX-X/18/06

(a) 0.25-spp input

(b) NSRR

(c) RAE

(d) Ours

(e) Ground truth

Figure 1.  Left to right: (a) noisy image generated using hybrid path-tracer
at 0.25 sample per pixel; (b) Neural supersampling network (Xiao et al., 2020)
(10.3ms at $1024\times 2048$, SSIM: 0.7737); (c) RAE (Chaitanya et al., 2017)
(6.5ms, SSIM: 0.7556); (d) our sparse sampling reconstruction (7.8ms, SSIM:
0.9036); (e) reference path-traced image with 32768 samples per pixel.

##  1\. Introduction

Rendering noise-free Monte Carlo (MC) ray-traced images at real-time frame
rates is still challenging. Despite the widely used of modern RTX GPU
accelerators, only a few rays per pixel can be traced at target resolution for
real-time applications, resulting in severe noise in renderings. The most
efficient strategy is to denoise and reconstruct the rendering results in
image-space, usually as a post-process pass of a physically-based renderer.

Until recently, most MC denoisers were proposed based on convolutional neural
networks (CNN). Chaitanya et al. (Chaitanya et al., 2017) proposed a recurrent
model for interactive applications that are targeted at images rendered with
low sample per pixel (1~4 spp). In addition, the NVIDIA OptiX ray-tracing
engine introduces an AI-accelerated denoiser based on this work. We also
developed a hybrid ray tracer based on Vulkan and we use it to export training
datasets. The source code of our ray tracer and paper will be available soon.
Mustafa et al. (Işık et al., 2021) adopt dilated spatial kernels to filter the
noisy image-guided by pairwise affinity over the features and target in the
low-sample count regime (2~8 spp). Meng et al. (Meng et al., 2020) also
denoise 1-spp noisy input images with a neural bilateral grid at real-time
frame rates. Hasselgren et al. (Hasselgren et al., 2020) proposed a neural
temporal adaptive sampling method for denoising image sequences rendered at
4-spp. Fan et al. (Fan et al., 2021) expands the kernel-prediction method to
remove noise at low spp (more than one) in a strict time budget. All the
state-of-the-art denoising and reconstruction methods aim at removing noise of
images rendered with more than 1-spp.

In this paper, we propose a novel approach to reconstruct less than 1-spp
renderings at real-time frame rates. Following traditional temporal anti-
aliasing (Karis, 2014) (TAA), our method uses renderer generated motion vector
to warp previous frames and accumulate sparse samples from the previous frame
based on the temporal accumulation factor computed according to the
correlation of current and previous frame, effectively increasing the number
of samples per pixel. The module can also detect ghosting artifacts at
disocclusion regions and remove mismatched pixels at inconsistent shading
regions. Mustafa et al. (Işık et al., 2021) also compute temporal accumulation
factor for a pixel using neural network, but they concatenate features of
current and previous frames and feed them into network together. Compared to
this method, our method can produce better temporal stable and high-quality
results.

After accumulating sparse samples, we use a residual block (He et al., 2015)
to fusion the accumulated features. Then we implement a multi-scale U-Net
(Ronneberger et al., 2015a) with skip connections for the reconstruction
subnetwork. The multi-scale predicting network is similar to the method
suggested by Vogels et al. (Vogels et al., 2018) which uses kernel prediction.
We directly predict denoised images for the current frame and two additional
channels as blending factors. We also predict a 2 $\times$ downscaled image
from the layer of the last but one. We composite the final denoised image from
the current denoised image, the previous warped image, and 2 $\times$
downscaled images. Comprehensive experiment results show that our approach is
good at reconstructing 0.25-spp images at a real-time frame rate. To
summarize, our contributions are the following:

  * •

We introduce a temporally-stable neural network to reconstruct image sequences
rendered at 0.25-spp at real-time frame rates. To the best of our knowledge,
we are the first that utilize 0.25-spp images as input for the neural network.

  * •

A novel temporal accumulation network which computes the correlation between
current and previous features to significantly improve the temporal stability
of Monte Carlo denoising.

  * •

Extensive experiments demonstrate that our method outperforms state-of-the-art
methods both quantitatively and qualitatively.

##  2\. RELATED WORK

Traditional best-performing MC denoisers were mainly based on local
neighborhood regression models (Zwicker et al., 2015). With the advent of
powerful modern GPUs, lots of researchers utilize CNN to build their MC
denoisers. In this section, we will mainly discuss CNN-based real-time
denoising techniques, which are most related to our approach. For a
comprehensive study of deep learning-based MC denoising and reconstruction
techniques, please refer to the recent survey of Huo et al. (Huo and Yoon,
2021).

###  2.1. Image-space Methods

Traditional MC denoisers are based on zero-order regression (Rousselle et al.,
2012; Delbracio et al., 2014; Li et al., 2012; Kalantari et al., 2015;
Rousselle et al., 2013; Moon et al., 2013), first-order regression (Bauszat et
al., 2011; Bitterli et al., 2016; Moon et al., 2014) and even higher-order
regression models (Moon et al., 2016).The filtering-based methods are based on
using the auxiliary feature buffers to guide the construction of image-space
filters. Most of the above methods run in offline rendering. To increase the
effective sample count, real-time denoisers leverage temporal accumulation
between frames over time to amortize supersampling (Yang et al., 2009), i.e.
temporal anti-aliasing (TAA). The previous frame is reprojected according to
the motion vector and blended with the current frame using a temporal
accumulation factor $\alpha$. The $\alpha$ can be constant (Schied et al.,
2017; Mara et al., 2017; Meng et al., 2020) and changed (Schied et al., 2018)
per frame and per pixel. The fixed temporal accumulation factor inevitably
leads to ghosting and temporal lag. By setting the parameter adaptively, the
temporal filter can fastly respond to times in case of sudden changes between
frames. Yang et al. (Yang et al., 2020) survey recent TAA techniques and
provide an in-depth analysis of the image quality trade-offs with these
heuristics. Koskela et al. (Koskela et al., 2019) propose a blockwise
regression for real-time path tracing reconstruction and also do accumulation
to improve temporal stability.

###  2.2. CNN-based Monte Carlo Denoising

Recent deep learning denoisers (Bako et al., 2017; Vogels et al., 2018) use
deep CNN to estimate the local per-pixel filtering kernels used to compute
each denoised pixel from its neighbors. Dahlberg et al. (Dahlberg et al.,
2019) implement the approach of (Vogels et al., 2018) as a practical
production tool used on the animated feature film. Layer-based denoiser
(Munkberg and Hasselgren, 2020) designs a hierarchical kernel prediction for
multi-resolution denoising and reconstruction. Since the high computational
cost of predicting large filtering kernels, these methods mostly target
offline renderings. There are also other methods (Kuznetsov et al., 2018; Xu
et al., 2019; Gharbi et al., 2019; Yu et al., 2021) that target denoising
rendering results at more than 4 spp.

To reduce the kernel prediction methods’ overhead, Fan et al. (Fan et al.,
2021) predict an encoding of the kernel map, followed by a high-efficiency
decoder to construct the complete kernel map. Chaitanya et al. (Chaitanya et
al., 2017) proposed a recurrent connection based on U-Net (Ronneberger et al.,
2015b) to improve temporal stability for sequences of sparsely sampled input
images. Hasselgren et al. (Hasselgren et al., 2020) proposed a neural spatio-
temporal joint optimization of adaptive sampling and denoising with a
recurrent feedback loop. Hofmann et al. (Hofmann et al., 2021) also utilized
the neural temporal adaptive sampling architecture to denoise rendering
results with participating media. Xiao et al. (Xiao et al., 2020) presented a
neural supersampling method for TAA, which is similar to deep-learned
supersampling (DLSS) (Edelsten et al., 2019). Meng et al. (Meng et al., 2020)
denoised 1-spp noisy input images with a neural bilateral grid at real-time
frame rates. Mustafa et al. (Işık et al., 2021) adopted dilated spatial
kernels to filter the noisy image guiding by pairwise affinity over the
features. Compare with these real-time denoising framework targeting for more
than 1-spp renderings, our method is designed to work with 0.25-spp.

##  3\. SPARSE SAMPLING DENOISING

###  3.1. Problem Statement

Our goal is to reconstruct temporally stable video from 0.25-spp hybrid path
traced image sequences in real-time frame rates, and we achieve this with a
supervised deep learning method. We use our hybrid path traced renderer to
generate a set of data
$\mathscr{D}$={($\textbf{c}_{1}$,$\textbf{f}_{1}$,$\textbf{r}_{1}$), …
,($\textbf{c}_{N}$,$\textbf{f}_{N}$,$\textbf{r}_{N}$)} where c stands for
noisy image rendered by sparse sampling, f is the auxiliary features (e.g.
albedo, normal, depth, metallic, roughness, shadow and transparent) obtained
in the rendering process, r is the reference image with high spp. We train a
deep neural function $\Phi$ with parameters $\Theta$ to reconstruct the noisy-
free image. The loss function $\ell$ is measured as the difference between the
denoised image and its reference image. We then minimize the loss function
with gradient descent algorithm across the dataset $\mathscr{D}$ with $N$
samples to get the optimal parameters $\hat{\theta}$:

(1)
$\hat{\theta}=\underset{\theta}{\arg\min}\sum_{i=1}^{N}\ell(\Phi(c_{i},\textbf{f}_{i}),r_{i})$

The loss function combines four-loss items, including spatial, temporal,
relative edge, and albedo loss, see section 3.4 for details.

###  3.2. Sparse Sampling

We developed a hybrid ray tracer to generate our dataset. To accelerate ray
tracing, we leverage a rasterization pipeline to get the first hit position
from the camera and store its associated shading attributes including albedo,
normal, depth, motion vector, metallic, and roughness. After this
rasterization pass, we trace a shadow ray to record the soft shadow attribute.
If there are transparent materials in the scene, we also save the transparent
attribute at the first hit position. We divide the full resolution into non-
overlapping blocks with spatial size $4\times 4$. We use the MC method to
solve the rendering equation (Kajiya, 1986) for one pixel in the block at each
frame and other pixels remain zero, see Figure 2. If the camera is static, the
radiance of all pixels will be computed once at every four frames. For image
sequences, we use two-layer CNNs to accumulate history frames, see section
3.3.1.

Figure 2. Sampling pattern. In _t_ and _t+1_ frame, we compute radiance for
the top left and top right pixel, respectively. In _t+2_ and _t+3_ frame,
bottom left and right pixel is estimated, respectively.

The input for our network is 18-channels features, including a 3D vector
(noised image, albedo, normal, shadow, and transparent) and a 1D vector
(depth, metallic, and roughness). Following prior method (Chaitanya et al.,
2017), we demodulate the noisy RGB image by the albedo of the directly visible
material, and the untextured irradiance $x$ is transformed to log space,
$\ln(1+x)$. Different from the prior method (Chaitanya et al., 2017), after
the untextured irradiance has been reconstructed, we re-modulate by the
accumulated albedo predicted by our temporal accumulator network which is our
key module for producing temporally stable results.

Figure 3. Network pipeline of our sparse sampling reconstruction (SSR) method.
The pipeline includes feature accumulator, feature fusion, and reconstruction
networks. The numbers under each network layer represent the output channels
at corresponding layers. The kernel size is $3\times 3$ at all layers. The
operator $\odot$ denotes dot product between features. ⓒ indicates
concatenation operation. $\oplus$ and $\otimes$ represent element-wise
addition and multiplication, respectively.

###  3.3. Network Pipeline

In this section, we describe our method in details with Figure 3.

####  3.3.1. Temporal Accumulator

The temporal accumulator module contains two neural networks each with 2-layer
CNNs. One network accepts normal and depth of current frame as input and
outputs reference embedding. Another network computes embeddings for the
current frame and warped the previous frame. These two embeddings are then
multiplied in a pixel-wise manner to the reference embedding and then call
softmax($\cdot$) to get $\mathbf{\alpha}$ and $\mathbf{\beta}$
($\mathbf{\alpha}+\mathbf{\beta}=1$) blending factors for current features and
previous features, respectively. We only accumulate noisy images, shadow and
albedo (see Figure 4). Take shadow as an example,we use the following equation
to accumulate shadow over the frame:

(2)
$\textbf{f}^{s}_{t}=\mathbf{\alpha}\mathcal{W}(\textbf{f}^{s}_{t-1})+\mathbf{\beta}\textbf{f}^{s}$

where $\textbf{f}^{s}_{t}$ is accumulated shadow until $t$ frame,
$\textbf{f}^{s}$ is shadow buffer for $t$ frame. $\mathcal{W}(\cdot)$ is a
warping operator that reprojects previous frame to current one using motion
vector. For the first frame, we set $\textbf{f}^{s}_{t-1}$ to
$\textbf{f}^{s}$.

(a) warped albedo (b) current albedo (c) accumulated albedo

Figure 4. The history albedo (a) is first wared and then is blended with the
current frame albedo (b). Our temporal accumulator not only fills missing
pixels but also smooths artifacts at the edge.

####  3.3.2. Feature Fusion

After accumulating images, shadow, and albedo, we concatenate accumulated
features, normal, depth, transparent, metallic, and roughness. Then we feed
them into a feature fusion network. Since our image is sparse, we use this
network to fusion the features and spread signals across spatial space.

####  3.3.3. Reconstruction Network

Finally, fused features and warped denoised images of the previous frame are
concatenated and fed into a reconstruction network, which outputs the high-
quality image for the current frame. The reconstruction network details are
given in Figure 3. Our network directly predict denoised fine image
$\textbf{d}^{f}$ for current frame and two additional channels as blending
factor, i.e., $\mathbf{\alpha}^{s}$ and $\mathbf{\alpha}^{t}$. We also
directly predict a 2 $\times$ downscaled coarse image $\textbf{d}^{c}$ from
the layer of the last but one. We use scale composition suggested by Vogels et
al (Vogels et al., 2018) to combine fine and coarse images:

(3)
$\textbf{O}_{p}=\textbf{d}^{f}_{p}-\mathbf{\alpha}^{s}_{p}[\textbf{UD}\textbf{d}^{f}]_{p}+\mathbf{\alpha}^{s}_{p}[\textbf{U}\textbf{d}^{c}]_{p}$

where D and U are $2\times 2$-downsampling and nearest-neighbor upsampling
operators. The filtered history $\textbf{O}_{t-1}$ is linearly blended with
the result of the scale composition O using $\mathbf{\alpha}^{t}$:

(4)
$\textbf{O}_{t}=\mathbf{\alpha}^{t}\textbf{O}+(1.0-\mathbf{\alpha}^{t})\textbf{O}_{t-1}$

###  3.4. Losses

We use the symmetric mean absolute percentage error (SMAPE):

(5)
$\ell(\textbf{r},\textbf{d})=\frac{1}{3N}\sum_{p=1}^{p=N}\sum_{c=1}^{c=3}\frac{\left|\textbf{d}_{p,c}-\textbf{r}_{p,c}\right|}{\left|\textbf{d}_{p,c}\right|+\left|\textbf{r}_{p,c}\right|+\varepsilon}$

Here, $N$ is the number of pixels in image and $\varepsilon$ is $10^{-2}$. d
and r are the denoised frame and the corresponding reference frame.

Our loss combines two parts, the first one is computed on a sequence of 5
images, including spatial loss $\ell_{s}=\ell(\textbf{r},\textbf{d})$,
temporal loss $\ell_{t}=\ell(\Delta\textbf{r},\Delta\textbf{d})$ where
$\Delta$ is temporal gradient computed between two consecutive frames,
relative edge loss
$\ell_{e}=L_{1}(\frac{\nabla\textbf{d}}{\textbf{r}+\varepsilon},\frac{\nabla\textbf{r}}{\textbf{r}+\varepsilon})$,
where gradient $\nabla$ is computed using a High Frequency Error Norm (HFEN),
an image comparison metric from medical imaging (Ravishankar and Bresler,
2011). As suggested by Chaitanya et al. (Chaitanya et al., 2017), we assign
higher weight to three loss functions ($\ell_{s}$, $\ell_{t}$ and $\ell_{e}$)
of frames later in the sequence to amplify temporal gradients. For our
training sequence of 5 images, we use (0.05, 0.25, 0.5, 0.75, 1).

The second part is warped temporal loss
$\ell_{wt}=\ell(\omega\textbf{r},\omega\textbf{d})$ where
$\omega\textbf{r}=r_{4}-\mathcal{W}(r_{3})$, $\mathcal{W}(\cdot)$ is a warping
operator that reprojects previous frame to current one. We also include albedo
loss $\ell_{a}=\ell(\textbf{a}_{acc},\textbf{a}_{r})$ where $\textbf{a}_{acc}$
is accumulated albedo computed by our feature accumulator network. We only
compute albedo loss on last frame and warped temporal loss on last two frames.

We use a weighted combination of these losses as the final training loss:

(6) $\ell=0.7\ell_{s}+0.1\ell_{t}+0.2\ell_{e}+0.4\ell_{wt}+5.0\ell_{a}$

##  4\. Dateset and Training Procedure

Since our method is designed for 3A game and virtual character rendering, we
train a separate network for each 3D scene same as (Xiao et al., 2020). Due to
the input image being generated at 0.25-spp, training robust denoiser requires
a large number of images. We train our method on 6 scenes, see Figure 5).
BistroInterior and BistroExterior (Lumberyard, 2017) have more than one
million triangles and support transparency, diffuse, specular, and soft shadow
features. Sponza, Diningroom, Angel, and Warmroom are simple scenes. Each
scene in the training set contains 100 to 1000 frames with resolution
$1024\times 2048$ depending on its complexity. We also rendered a validation
set with 10 frames and a test set with $50$ frames for each scene. For each
frame, we rendered the reference image at 32768 spp which is the target of our
denoiser.

(a) BistroInterior

(b) BistroExterior

(c) Sponza

(d) Diningroom (e) Angel (f) Warmroom

Figure 5. An overview of reference images in our generated dataset

When we train the denoiser, we randomly select $5$ consecutive frames for
training in consecutive clips of each scene. The inputs, including the noisy
image and the auxiliary features, of each frame, are randomly cropped with
resolution $256\times 256$ to make full use of GPU.

We optimize our denoiser network with ADAM optimizer (Kingma and Ba, 2015). We
set the initial learning rate to $1\times 10^{-4}$ and half it at one-third
and two-thirds of the total number of iterations. The batch size is 7 and the
epoch is 200 for each scene. Our denoiser is implemented by PyTorch (Paszke et
al., 2019) and all the models we presented were trained and tested parallel on
four GPUs of NVIDIA Tesla A100. Each network takes around 9 hours.

##  5\. Results

In this section, we evaluate the performance of our method. We describe the
implementation of compared baseline and metrics in Section 5.1, analyze the
algorithm with various ablation experiments in Section 5.2, and describe its
limitations and future work in Section 5.3.

###  5.1. Baseline and metrics

We compare our method with several state-of-the-art denoising and
reconstruction work, including real-time methods RAE (Chaitanya et al., 2017),
ANF (Işık et al., 2021), offline method MCD (Yu et al., 2021), and super-
resolution model NSRR (Xiao et al., 2020). Although NSRR is a method for the
super-resolution task, it can also reconstruct images from zero-padding inputs
which means it fits the sparse sampling task well. So we choose it as one of
our competitors. We removed the zero-sampling modules so that it can apply to
our dataset. We follow all these papers and use PyTorch to re-implement them.
We train all the methods on the same datasets as in our method with the same
training procedure.

To evaluate quality, we use three quality metrics: peak signal to noise ratio
(PSNR), structural similarity index (SSIM) (Wang et al., 2004), and root mean
squared error (RMSE). The higher the better in both PSNR and SSIM, while the
lower the better in RMSE.

Quantitative comparison results are shown in Table1. Average results are
reported on the 50 test videos of six scenes. We only show the results of SSIM
due to space limit, and please refer to our supplemental material for more
comparison results. As shown in Table 1, our method achieves the best
performance with all six scenes.

Scene | MCD | ANF | NSRR | RAE | SSR
---|---|---|---|---|---
BistroInterior | 0.7650 | 0.7583 | 0.7405 | 0.7751 | 0.8921
BistroExterior | 0.8071 | 0.7201 | 0.8538 | 0.8006 | 0.8962
Sponza | 0.8119 | 0.8219 | 0.8113 | 0.8898 | 0.9410
Diningroom | 0.8637 | 0.7226 | 0.8843 | 0.9007 | 0.9375
Warmroom | 0.8021 | 0.8774 | 0.9740 | 0.9675 | 0.9758
Angel | 0.8601 | 0.8813 | 0.9804 | 0.9161 | 0.9763
Table 1. Quantitative comparison results on six scenes. We choose four
baseline methods to compare with our SSR method.

At inference time, all methods are applied to a single frame at a time, Table
2 shows inference time at 1024 $\times$ 2048 resolution. We tested all the
models using a GPU, NVIDIA Tesla A100. All network models are not optimized
with Nvidia TensorRT at 16-bit precision, so inference time still has room for
improvement.

Method | MCD | ANF | NSRR | RAE | SSR
---|---|---|---|---|---
Time(ms) | 13.5 | 32 | 33.5 | 6.5 | 7.8
Table 2. Comparison results of inference time.

In Figure 6, we compare reconstructed images visually. Our method outperforms
all other methods on all scenes by a large margin. Previous state-of-the-art
methods are not good at denoising renderings at 0.25-spp. MCD originally
targets offline rendering and transformer needs large memory to train and
inference. RAE, NSRR, and ANF feed previous and current features into the
network directly. The difference between our approach and the previous ones is
that we compute the correlation for each pixel between normal and depth
features of current and previous frame. Please refer to the supplementary
material and videos, our method produces significantly more temporally stable
video results than existing methods.

(a) Ours

(b) Input

(c) MCD

(d) ANF

(e) NSRR

(f) RAE

(g) Ours

(h) Reference

Figure 6. Visual results on BistroInterior, BistroExterior, Sponza,
Diningroom, Warmroom, and Angel scenes.

###  5.2. Analysis

####  5.2.1. Rendering Efficiency

We test rendering time of each stage in NVIDIA RTX 3060 GPU at resolution
$1024\times 2048$, see Table 3. With our sparse sampling, the total rendering
time of scene BistroInterior is 8.75 ms. Without sparse sampling, the total
rendering time is 18.19 ms. This leads to an about $3\times$ rendering
performance improvement. After applying our SSR model, high-fidelity results
are produced.

Rasterization | Transparent and Shadow | W-SS | Wo-SS
---|---|---|---
1.08 ms | 2.32 ms | 5.35 ms | 14.79 ms
Table 3. Rendering time of scene BistroInterior. W-SS means rendering with our
sparse sampling, Wo-SS means without sparse sampling.

####  5.2.2. Quality Gain with Shadow and Transparent

Our training images are produced by MC path tracer at 0.25-spp average. Due to
light occlusion, more than three-fourths of pixels remain zero, so we need
more features to train our model. We add direct noisy shadow as input of our
model. Our feature accumulator will accumulate noisy shadows between the
current frame and the history shadow buffer. The accumulated shadow can help
our model to detect the continuous edge of shadow and improve temporal
stability. The synthesis video of test sequences can show that the edge of
Figure 7 without shadow feature will jitter over frames. If noisy features
feed into regression-based method (Rousselle et al., 2013), the quality of the
denoised image will decrease. These methods need another filter to prefilter
noisy features, but CNN-based methods can accept more than one noisy buffer
except noisy images.

(a) Wo-shadow (b) Ours (c) Noisy shadow (d) Ground truth

Figure 7.  The result (a) is generated by training model without shadow
feature, our result (b) is trained with shadow feature. (c) and (d) is noisy
shadow feature and ground truth, respectively.

We also add the transparent feature into our model for training, but we did
not accumulate it before feeding it into the feature fusion module. The reason
is that the transparent feature includes less noise than the shadow, see
Figure 8. If the scene didn’t have a transparent object, such as
BistroExterior, we also feed transparent features with zero.

Our model without shadow and transparent only gets 27.98 dB on testing
BistroInterior. With shadow and transparent, our model not only gets higher
PSNR (28.78 dB) but also generates high-quality image.

(a) Wo-transparent (b) W-transparent (c) Transparent (d) Ground truth

Figure 8.  The result (a) is generated by training our model without the
transparent feature, our result (b) is trained with the transparent feature.
(c) and (d) is the transparent feature and ground truth, respectively.

####  5.2.3. Quality Gain with Feature Accumulator

We demodulate the image with the albedo at a primary hit position. After our
network reconstructs the untextured illumination, we re-modulate by the albedo
to include the texture detail in the final rendering. If the albedo in the
corresponding frame has an artifact, the artifact will transfer to the final
rendering. Chaitanya et. al (Chaitanya et al., 2017) apply TAA as a
supplemental post-process pass to fix the artifact. We re-modulate by the
accumulated albedo generated by the temporal accumulator module to achieve
some efficiency as multisampling antialiasing (MSAA). (Akeley, 1993).

In summary, from Figure 4, Figure 7 and Figure 9, we can see that our feature
accumulator plays a key role in reconstructing sparse sampling renderings at
less than 1-spp.

(a) SSIM: 0.8626 (b) SSIM: 0.8972 (c) Ground truth

Figure 9.  (a) Artifact is transferred to the final result, SSIM is 0.8626 (b)
Re-modulating by the accumulated albedo leads to high-quality image, SSIM is
0.8927.

####  5.2.4. Network Modules

In Table 4, we report the ablation experiments for analyzing the quality
improvements from the temporal accumulator (Section 3.3.1) and feature fusion
(Section 3.3.2) modules. Average results are reported on the 50 test image
sequences of the BistroInterior scene. If without the temporal accumulator and
feature fusion, the PSNR decreases about 0.58dB, but temporal stability will
decrease dramatically in the video results. See our supplemental materials for
more detailed information.

Feature Accumulator | Feature Fusion | SSIM | PSNR (dB)
---|---|---|---
✕ | ✕ | 0.8600 | 28.20
✕ | ✓ | 0.8617 | 28.15
✓ | ✕ | 0.8866 | 28.56
✓ | ✓ | 0.8911 | 28.78
Table 4. Ablation experiment for the feature accumulator and feature fusion
modules. The network is trained with each (and both) of these subnetworks
removed, and results on the BistroInterior scene are reported.

###  5.3. Limitations and Future Work

While our method provides a significant improvement for neural sparse sampling
reconstruction, the inference time still has room for improvement. We will
adopt TensorRT for acceleration and deploy our model on our game engine and
virtual character rendering platform in the future. In addition, there is
still little jitter for small objects on the temporal domain. Modern game
engine all has a TAA pass, applying TAA post-processing can get more
temporally stable results. We also try to add a layer of Swin Transformer (Liu
et al., 2021) to the first layer of our reconstruction network. It can truely
improve quantitative number about 0.23 dB, but inference time will increase
1.1 ms at 1024 $\times$ 2048 resolution on NVIDIA Tesla A100.

##  6\. CONCLUSION

We have presented the first CNN-based method for reconstructing Monte Carlo
renderings at 0.25-spp and experiments show that our method reconstructs high-
quality results compared with current state-of-the-art methods. We propose an
efficient feature accumulator network to compute the blending factor for each
pixel between current and previous frames. Then the accumulated features are
fused and fed into a multi-scale U-Net to reconstruct final results. We
evaluated our method by comparing its performance to previous works
demonstrating better results across all test scenes.

##  7\. Acknowledgments

We thank Open Research Content Archive (ORCA) of NVIDIA for providing
BistroInterior and BistroExterior scenes for training and testing. We also
thank all students who participate in the development of our hybrid renderer.

## References

  * (1)
  * Akeley (1993) Kurt Akeley. 1993\.  Reality Engine Graphics. In _Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques_ (Anaheim, CA) _(SIGGRAPH ’93)_. Association for Computing Machinery, New York, NY, USA, 109–116.  https://doi.org/10.1145/166117.166131
  * Bako et al. (2017) Steve Bako, Thijs Vogels, Brian Mcwilliams, Mark Meyer, Jan NováK, Alex Harvill, Pradeep Sen, Tony Derose, and Fabrice Rousselle. 2017.  Kernel-Predicting Convolutional Networks for Denoising Monte Carlo Renderings.  _ACM Trans. Graph._ 36, 4, Article 97 (jul 2017), 14 pages.  https://doi.org/10.1145/3072959.3073708
  * Bauszat et al. (2011) Pablo Bauszat, Martin Eisemann, and Marcus Magnor. 2011\.  Guided Image Filtering for Interactive High-Quality Global Illumination. In _Proceedings of the Twenty-Second Eurographics Conference on Rendering_ (Prague, Czech Republic) _(EGSR ’11)_. Eurographics Association, Goslar, DEU, 1361–1368.  https://doi.org/10.1111/j.1467-8659.2011.01996.x
  * Bitterli et al. (2016) Benedikt Bitterli, Fabrice Rousselle, Bochang Moon, José A. Iglesias-Guitián, David Adler, Kenny Mitchell, Wojciech Jarosz, and Jan Novák. 2016.  Nonlinearly Weighted First-Order Regression for Denoising Monte Carlo Renderings.  _Comput. Graph. Forum_ 35, 4 (jul 2016), 107–117.
  * Chaitanya et al. (2017) Chakravarty Chaitanya, Anton Kaplanyan, Christoph Schied, Marco Salvi, Aaron Lefohn, Derek Nowrouzezahrai, and Timo Aila. 2017\.  Interactive reconstruction of Monte Carlo image sequences using a recurrent denoising autoencoder.  _ACM Transactions on Graphics_ 36 (07 2017), 1–12.  https://doi.org/10.1145/3072959.3073601
  * Dahlberg et al. (2019) Henrik Dahlberg, David Adler, and Jeremy Newlin. 2019\.  Machine-Learning Denoising in Feature Film Production. In _ACM SIGGRAPH 2019 Talks_ (Los Angeles, California) _(SIGGRAPH ’19)_. Association for Computing Machinery, New York, NY, USA, Article 21, 2 pages.  https://doi.org/10.1145/3306307.3328150
  * Delbracio et al. (2014) Mauricio Delbracio, Pablo Musé, Julien Chauvier, Nicholas Phelps, and Jean-Michel Morel. 2014.  Boosting Monte Carlo Rendering by Ray Histogram Fusion.  _ACM Transactions on Graphics (TOG)_ 33 (02 2014).  https://doi.org/10.1145/2532708
  * Edelsten et al. (2019) Andrew Edelsten, Paula Jukarainen, and Anjul Patney. 2019\.  Truly next-gen: Adding deep learning to games and graphics.  _In NVIDIA Sponsored Sessions (Game Developers Conference)_ (2019).
  * Fan et al. (2021) Hangming Fan, Rui Wang, Yuchi Huo, and Hujun Bao. 2021\.  Real-time Monte Carlo Denoising with Weight Sharing Kernel Prediction Network.  _Computer Graphics Forum_ 40, 4 (2021), 15–27.  https://doi.org/10.1111/cgf.14338 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.14338
  * Gharbi et al. (2019) Michaël Gharbi, Tzu-Mao Li, Miika Aittala, Jaakko Lehtinen, and Frédo Durand. 2019\.  Sample-Based Monte Carlo Denoising Using a Kernel-Splatting Network.  _ACM Trans. Graph._ 38, 4, Article 125 (jul 2019), 12 pages.  https://doi.org/10.1145/3306346.3322954
  * Hasselgren et al. (2020) Jon Hasselgren, Jacob Munkberg, Marco Salvi, Anjul Patney, and Aaron Lefohn. 2020.  Neural temporal adaptive sampling and denoising. In _Computer Graphics Forum_ , Vol. 39. Wiley Online Library, 147–155.
  * He et al. (2015) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. 2015\.  Deep Residual Learning for Image Recognition.  _CoRR_ abs/1512.03385 (2015).  arXiv:1512.03385 http://arxiv.org/abs/1512.03385
  * Hofmann et al. (2021) Nikolai Hofmann, Jon Hasselgren, Petrik Clarberg, and Jacob Munkberg. 2021.  Interactive Path Tracing and Reconstruction of Sparse Volumes.  4, 1, Article 5 (apr 2021), 19 pages.  https://doi.org/10.1145/3451256
  * Huo and Yoon (2021) Yuchi Huo and Sung-eui Yoon. 2021.  A survey on deep learning-based Monte Carlo denoising.  _Computational Visual Media_ 7 (03 2021).  https://doi.org/10.1007/s41095-021-0209-9
  * Işık et al. (2021) Mustafa Işık, Krishna Mullia, Matthew Fisher, Jonathan Eisenmann, and Michaël Gharbi. 2021.  Interactive Monte Carlo Denoising Using Affinity of Neural Features.  _ACM Trans. Graph._ 40, 4, Article 37 (jul 2021), 13 pages.  https://doi.org/10.1145/3450626.3459793
  * Kajiya (1986) James T. Kajiya. 1986\.  The Rendering Equation. In _Proceedings of the 13th Annual Conference on Computer Graphics and Interactive Techniques_ _(SIGGRAPH ’86)_. Association for Computing Machinery, New York, NY, USA, 143–150.  https://doi.org/10.1145/15922.15902
  * Kalantari et al. (2015) Nima Khademi Kalantari, Steve Bako, and Pradeep Sen. 2015\.  A Machine Learning Approach for Filtering Monte Carlo Noise.  _ACM Trans. Graph._ 34, 4, Article 122 (jul 2015), 12 pages.  https://doi.org/10.1145/2766977
  * Karis (2014) Brian Karis. 2014\.  High-Quality Temporal Supersampling.  _In ACM SIGGRAPH 2014 Courses: Advances in Real-time Rendering in Games, Part I (SIGGRAPH ’14)_ (2014).  https://doi.org/10.1145/2614028.2615455
  * Kingma and Ba (2015) Diederik P. Kingma and Jimmy Ba. 2015.  Adam: A Method for Stochastic Optimization. In _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , Yoshua Bengio and Yann LeCun (Eds.).  http://arxiv.org/abs/1412.6980
  * Koskela et al. (2019) Matias Koskela, Kalle Immonen, Markku Mäkitalo, Alessandro Foi, Timo Viitanen, Pekka Jääskeläinen, Heikki Kultala, and Jarmo Takala. 2019.  Blockwise Multi-Order Feature Regression for Real-Time Path-Tracing Reconstruction.  38, 5, Article 138 (jun 2019), 14 pages.  https://doi.org/10.1145/3269978
  * Kuznetsov et al. (2018) Alexandr Kuznetsov, Nima Khademi Kalantari, and Ravi Ramamoorthi. 2018\.  Deep Adaptive Sampling for Low Sample Count Rendering.  _Computer Graphics Forum_ 37 (07 2018), 35–44.  https://doi.org/10.1111/cgf.13473
  * Li et al. (2012) Tzu-Mao Li, Yu-Ting Wu, and Yung-Yu Chuang. 2012.  SURE-Based Optimization for Adaptive Sampling and Reconstruction.  _ACM Trans. Graph._ 31, 6, Article 194 (nov 2012), 9 pages.  https://doi.org/10.1145/2366145.2366213
  * Liu et al. (2021) Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo. 2021.  Swin Transformer: Hierarchical Vision Transformer using Shifted Windows.  _CoRR_ abs/2103.14030 (2021).  arXiv:2103.14030 https://arxiv.org/abs/2103.14030
  * Lumberyard (2017) Amazon Lumberyard. 2017\.  Amazon Lumberyard Bistro, Open Research Content Archive (ORCA).  http://developer.nvidia.com/orca/amazon-lumberyard-bistro http://developer.nvidia.com/orca/amazon-lumberyard-bistro.
  * Mara et al. (2017) Michael Mara, Morgan McGuire, Benedikt Bitterli, and Wojciech Jarosz. 2017.  An Efficient Denoising Algorithm for Global Illumination. In _Proceedings of High Performance Graphics_ (Los Angeles, California) _(HPG ’17)_. Association for Computing Machinery, New York, NY, USA, Article 3, 7 pages.  https://doi.org/10.1145/3105762.3105774
  * Meng et al. (2020) Xiaoxu Meng, Quan Zheng, Amitabh Varshney, Gurprit Singh, and Matthias Zwicker. 2020.  Real-time Monte Carlo Denoising with the Neural Bilateral Grid. In _Eurographics Symposium on Rendering - DL-only Track_ , Carsten Dachsbacher and Matt Pharr (Eds.). The Eurographics Association.  https://doi.org/10.2312/sr.20201133
  * Moon et al. (2014) Bochang Moon, Nathan Carr, and Sung-Eui Yoon. 2014\.  Adaptive Rendering Based on Weighted Local Regression.  _ACM Trans. Graph._ 33, 5, Article 170 (sep 2014), 14 pages.  https://doi.org/10.1145/2641762
  * Moon et al. (2013) Bochang Moon, Jong Yun Jun, JongHyeob Lee, Kunho Kim, Toshiya Hachisuka, and Sung-Eui Yoon. 2013\.  Robust Image Denoising Using a Virtual Flash Image for Monte Carlo Ray Tracing.  _Computer Graphics Forum_ 32, 1 (2013), 139–151.  https://doi.org/10.1111/cgf.12004 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.12004
  * Moon et al. (2016) Bochang Moon, Steven McDonagh, Kenny Mitchell, and Markus Gross. 2016.  Adaptive Polynomial Rendering.  _ACM Trans. Graph._ 35, 4, Article 40 (jul 2016), 10 pages.  https://doi.org/10.1145/2897824.2925936
  * Munkberg and Hasselgren (2020) Jacob Munkberg and Jon Hasselgren. 2020.  Neural Denoising with Layer Embeddings.  _Computer Graphics Forum_ 39, 4 (2020), 1–12.  https://doi.org/10.1111/cgf.14049 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.14049
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Z. Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. 2019.  PyTorch: An Imperative Style, High-Performance Deep Learning Library.  _CoRR_ abs/1912.01703 (2019).  arXiv:1912.01703 http://arxiv.org/abs/1912.01703
  * Ravishankar and Bresler (2011) Saiprasad Ravishankar and Yoram Bresler. 2011.  MR Image Reconstruction From Highly Undersampled k-Space Data by Dictionary Learning.  _IEEE Transactions on Medical Imaging_ 30, 5 (2011), 1028–1041.  https://doi.org/10.1109/TMI.2010.2090538
  * Ronneberger et al. (2015a) Olaf Ronneberger, Philipp Fischer, and Thomas Brox. 2015a.  U-Net: Convolutional Networks for Biomedical Image Segmentation.  _CoRR_ abs/1505.04597 (2015).  arXiv:1505.04597 http://arxiv.org/abs/1505.04597
  * Ronneberger et al. (2015b) Olaf Ronneberger, Philipp Fischer, and Thomas Brox. 2015b.  U-Net: Convolutional Networks for Biomedical Image Segmentation.  _CoRR_ abs/1505.04597 (2015).  arXiv:1505.04597 http://arxiv.org/abs/1505.04597
  * Rousselle et al. (2012) Fabrice Rousselle, Claude Knaus, and Matthias Zwicker. 2012\.  Adaptive Rendering with Non-Local Means Filtering.  _ACM Trans. Graph._ 31, 6, Article 195 (nov 2012), 11 pages.  https://doi.org/10.1145/2366145.2366214
  * Rousselle et al. (2013) Fabrice Rousselle, Marco Manzi, and Matthias Zwicker. 2013\.  Robust Denoising using Feature and Color Information.  _Computer Graphics Forum_ 32, 7 (2013), 121–130.  https://doi.org/10.1111/cgf.12219 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.12219
  * Schied et al. (2017) Christoph Schied, Anton Kaplanyan, Chris Wyman, Anjul Patney, Chakravarty R. Alla Chaitanya, John Burgess, Shiqiu Liu, Carsten Dachsbacher, Aaron Lefohn, and Marco Salvi. 2017.  Spatiotemporal Variance-Guided Filtering: Real-Time Reconstruction for Path-Traced Global Illumination. In _Proceedings of High Performance Graphics_ (Los Angeles, California) _(HPG ’17)_. Association for Computing Machinery, New York, NY, USA, Article 2, 12 pages.  https://doi.org/10.1145/3105762.3105770
  * Schied et al. (2018) Christoph Schied, Christoph Peters, and Carsten Dachsbacher. 2018.  Gradient Estimation for Real-Time Adaptive Temporal Filtering.  _Proc. ACM Comput. Graph. Interact. Tech._ 1, 2, Article 24 (aug 2018), 16 pages.  https://doi.org/10.1145/3233301
  * Vogels et al. (2018) Thijs Vogels, Fabrice Rousselle, Brian Mcwilliams, Gerhard Röthlin, Alex Harvill, David Adler, Mark Meyer, and Jan Novák. 2018\.  Denoising with Kernel Prediction and Asymmetric Loss Functions.  _ACM Trans. Graph._ 37, 4, Article 124 (jul 2018), 15 pages.  https://doi.org/10.1145/3197517.3201388
  * Wang et al. (2004) Zhou Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli. 2004\.  Image quality assessment: from error visibility to structural similarity.  _IEEE Transactions on Image Processing_ 13, 4 (2004), 600–612.  https://doi.org/10.1109/TIP.2003.819861
  * Xiao et al. (2020) Lei Xiao, Salah Nouri, Matt Chapman, Alexander Fix, Douglas Lanman, and Anton Kaplanyan. 2020\.  Neural Supersampling for Real-Time Rendering.  _ACM Trans. Graph._ 39, 4, Article 142 (jul 2020), 12 pages.  https://doi.org/10.1145/3386569.3392376
  * Xu et al. (2019) Bing Xu, Junfei Zhang, Rui Wang, Kun Xu, Yong-Liang Yang, Chuan Li, and Rui Tang. 2019.  Adversarial Monte Carlo Denoising with Conditioned Auxiliary Feature Modulation.  _ACM Trans. Graph._ 38, 6, Article 224 (nov 2019), 12 pages.  https://doi.org/10.1145/3355089.3356547
  * Yang et al. (2020) Lei Yang, Shiqiu Liu, and Marco Salvi. 2020.  A Survey of Temporal Antialiasing Techniques.  _Computer Graphics Forum_ 39, 2 (2020), 607–621.  https://doi.org/10.1111/cgf.14018 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.14018
  * Yang et al. (2009) Lei Yang, Diego Nehab, Pedro V. Sander, Pitchaya Sitthi-amorn, Jason Lawrence, and Hugues Hoppe. 2009\.  Amortized Supersampling.  _ACM Trans. Graph._ 28, 5 (dec 2009), 1–12.  https://doi.org/10.1145/1618452.1618481
  * Yu et al. (2021) Jiaqi Yu, Yongwei Nie, Chengjiang Long, Wenju Xu, Qing Zhang, and Guiqing Li. 2021\.  Monte Carlo Denoising via Auxiliary Feature Guided Self-Attention.  _ACM Trans. Graph._ 40, 6, Article 273 (dec 2021), 13 pages.  https://doi.org/10.1145/3478513.3480565
  * Zwicker et al. (2015) M. Zwicker, W. Jarosz, J. Lehtinen, B. Moon, R. Ramamoorthi, F. Rousselle, P. Sen, C. Soler, and S.-E. Yoon. 2015\.  Recent Advances in Adaptive Sampling and Reconstruction for Monte Carlo Rendering.  _Comput. Graph. Forum_ 34, 2 (may 2015), 667–681.