aainstitutetext: Rudolf Peierls Centre for Theoretical Physics, Clarendon
Laboratory, University of Oxford, Oxford OX1 3PU, United
Kingdombbinstitutetext: Institut für Theoretische Teilchenphysik und
Kosmologie, RWTH Aachen University, D-52056 Aachen, Germanyccinstitutetext:
Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United
Kingdom111Preprint: Cavendish-HEP-21/07, OUTP-21-13P, P3H-21-032

# NNLO QCD corrections to diphoton production with an additional jet at the
LHC

Herschel A. Chawdhry b Michał Czakon c Alexander Mitov c Rene Poncelet

###### Abstract

We calculate the NNLO QCD corrections to diphoton production with an
additional jet at the LHC. Our calculation represents the first NNLO-accurate
prediction for the transverse momentum distribution of the diphoton system.
The improvement in the accuracy of the theoretical prediction is significant,
by a factor of up to four relative to NLO QCD. Our calculation is exact except
for the finite remainder of the two-loop amplitude which is included at
leading color. The numerical impact of this approximated contribution is
small. The results of this work are expected to further our understanding of
the Higgs boson sector and of the behavior of higher-order corrections to LHC
processes.

##  1 Introduction

The production of a pair of photons at the LHC is of special interest. On the
one hand this process represents the main background to the cleanest Higgs
boson decay channel $h\to{\gamma\gamma}$. On the other hand it is a process
where our ability to accurately predict LHC cross sections by including
higher-order QCD corrections can reliably be tested.

Inclusive diphoton production $pp\to{\gamma\gamma}+X$ has been studied
extensively at NNLO in QCD Catani:2011qz ; Campbell:2016yrh ; Catani:2018krb ;
Gehrmann:2020oec ; Alioli:2020qrd . Prior work also includes NLO QCD
corrections DelDuca:2003uz  or electroweak effects Chiesa:2017gqx , as well as
photon isolation Gehrmann:2013aga . Beyond NLO, resummation effects have been
included at least to NNLL Cieri:2015rqa ; Coradeschi:2017zzw ; Becher:2020ugp
; Ebert:2016gcn . Interference between ${\gamma\gamma+{\rm jet}}$ and
$h\to{\gamma\gamma}+j$ have also been investigated Cieri:2017kpq . A detailed
analysis of the limitations of existing results can be found in ref.
Amoroso:2020lgh .

A distinguishing feature of this process is the presence of very large higher-
order QCD corrections which has raised some questions about the reliability of
such higher-order predictions for this process. With time our understanding of
the behavior of higher-order corrections in this case has developed
significantly. It is presently believed Catani:2018krb  that starting with
N3LO, perturbative corrections will be much more mild and consistent with
perturbative convergence. A similar conclusion has also been reached recently
for the process $pp\to{\gamma\gamma\gamma}$ Chawdhry:2019bji ;
Kallweit:2020gcp . Clearly it is very desirable to have a full N3LO accurate
calculation of inclusive diphoton production where these ideas can be tested
and hopefully validated. The present work, together with the recently computed
3-loop amplitudes for diphoton production at the LHC Caola:2020dfu ,
represents a significant step in this direction.

The transverse momentum of the photon pair, ${p_{T}(\gamma\gamma)}$, plays a
special role in inclusive diphoton production $pp\to{\gamma\gamma}+X$. As is
well known, due to the fact that at leading order ${p_{T}(\gamma\gamma)}=0$,
an NNLO-accurate calculation of inclusive diphoton production is only NLO-
accurate for the ${p_{T}(\gamma\gamma)}$ distribution. To achieve NNLO
accuracy for ${p_{T}(\gamma\gamma)}$ at nonzero ${p_{T}(\gamma\gamma)}$ one
needs to compute the NNLO QCD corrections for the process
$pp\to{\gamma\gamma+{\rm jet}}+X$. The present work presents the first
calculation of $pp\to{\gamma\gamma+{\rm jet}}+X$ in NNLO QCD and makes public
NNLO QCD predictions for a number of diphoton observables at nonzero
${p_{T}(\gamma\gamma)}$.

The main reason the $p_{T}$ distribution of the diphoton system is of special
interest is that it represents the main background for Higgs production at
high $p_{T}$. High-$p_{T}$ Higgs production is relevant for Dark Matter
searches Sirunyan:2018fpy ; Aad:2021qks  and for disentangling the nature of
the Higgs boson’s local vertex Lindert:2018iug  which is not possible at low
$p_{T}$ where the effective $ggh$ vertex describes Higgs production well. The
${p_{T}(\gamma\gamma)}$ distribution, together with the angular distribution
of the two photons in the Collins-Soper frame Collins:1977iv , represents a
strong discriminator for the spin of a possible resonance decaying to two
photons Aad:2015mxa . For further details about high-$p_{T}$ Higgs production
we refer the reader to the recent review Becker:2020rjp .

This work is organized as follows: in sec. 2 we briefly describe our
calculation while in sec. 3 we present our predictions for a number of
differential distributions. Our conclusions are given in sec. 4.

##  2 Setup of the Calculation

The calculation is performed in the STRIPPER approach Czakon:2010td ;
Czakon:2011ve ; Czakon:2014oma . The approach has already been applied in the
calculation of NNLO QCD corrections to top-pair Czakon:2014xsa ;
Czakon:2015owf ; Czakon:2016ckf ; Czakon:2016dgf ; Behring:2019iiv , inclusive
jet Czakon:2019tmo , three-photon Chawdhry:2019bji , $W+c~{}{\rm jet}$
Czakon:2020coa , identified $B$-hadron Czakon:2021ohs  and polarized $W$-pair
Poncelet:2021jmj  production at the LHC. A detailed description of the
technical aspects of our implementation has been given in ref. Czakon:2019tmo
.

All tree-level diagrams are computed with the avhlib library avhlib ;
Bury:2015dla . The contributing one-loop amplitudes, including the loop-
induced contribution, are obtained from the library OpenLoops Cascioli:2011va
; Buccioni:2019sur .

The relevant two-loop contributions $q\bar{q}\to g{\gamma\gamma}$ and $qg\to
q{\gamma\gamma}$ are handled in the following way. We first separate the
finite remainders ${\cal H}^{(2)}(\mu_{R}^{2})$ of the two-loop amplitudes
from their infrared poles. The latter can be predicted exactly and we have
included them, including their finite contributions, without any
approximation. In terms of the scale dependence of the two-loop finite
remainder ${\cal H}^{(2)}$ (defined as in ref. Czakon:2014oma ):

$\displaystyle{\cal H}^{(2)}(\mu_{R}^{2})={\cal
H}^{(2)}(s_{12})+\sum_{i=1}^{4}c_{i}\ln^{i}\left({\mu_{R}^{2}\over
s_{12}}\right)\,,$ (1)

where $s_{12}$ is the squared partonic center-of-mass energy, we have included
without any approximation all two-loop terms $c_{i}$ corresponding to
$\ln^{i}(\mu_{R}^{2})$ with $i\geq 1$. The scale-independent part ${\cal
H}^{(2)}(s_{12})$ is included in the leading color approximation as derived in
ref. Chawdhry:2021mkw  with the help of refs. Chawdhry:2018awn ;
Chicherin:2020oor ; Chawdhry:2020for  (an equivalent expression for the spin-
averaged two-loop squared amplitude has also been derived in ref.
Agarwal:2021grm ). This is the only approximation made in this paper. Further
details about the implementation of the leading color approximation of the
two-loop finite remainder can be found in ref. Chawdhry:2019bji .

As a justification of the leading color approximation just described we have
verified that the numerical contribution of the scale-independent part of the
two-loop finite remainder ${\cal H}^{(2)}(s_{12})$ is small. For all
distributions computed here it is about 1-2% of the complete NNLO prediction
in all bins. This makes it smaller than the Monte Carlo integration error in
the differential distributions. More details can be found in sec. 3.

We have also included the so-called loop-induced (LI) contribution $gg\to
g{\gamma\gamma}$ which begins to contribute starting at NNLO. As we explain in
sec. 3 its effect is about five percent and is strongly dependent on the
distribution.

The QED coupling is taken to be $\alpha=1/137$. The strong coupling constant
and parton distributions are renormalized such that they evolve with $n_{f}=5$
active flavors. Diagrams with top-quark loops are included in all
contributions except the two-loop finite remainder and the one-loop squared
contribution. The current calculation uses the NNPDF3.1 pdf set Ball:2017nwa
of order that matches the order of the perturbative calculation. The value of
the strong coupling constant is taken from the LHAPDF interface
Buckley:2014ana . The central values of the factorization and renormalization
scales have been fixed to:

$\mu_{F}^{2}=\mu_{R}^{2}=\frac{1}{4}\left(m^{2}(\gamma\gamma)+p_{T}(\gamma\gamma)^{2}\right)\,.$
(2)

This scale choice has been motivated by the scales used in refs.
DelDuca:2003uz ; Gehrmann:2013aga : we have replaced the $p_{T}$ of the jet
with that of the photon pair since in our setup, see below, we have no
explicit jet requirements. At LO the two scales are equivalent.

Scale uncertainty has been estimated with the help of a 7-point restricted
independent variation by a factor of 2 of the scales $\mu_{F}$ and $\mu_{R}$.
Since in this work we are primarily concerned with perturbative convergence
and estimates of missing higher-order corrections, we have not included pdf
error estimates. We expect that those are not dominant over the scale
variation in the kinematic ranges considered here. We hope to include them in
a future update of the present work.

Our calculation has been performed for the LHC at 13 TeV and is subject to the
following set of selection cuts: 222This set of cuts is based on typical
selection requirements, see refs. Aad:2015mxa ; Sirunyan:2019iwo . we require
two photons satisfying the following criteria

  * •

$p_{T}(\gamma_{1})>30$ GeV, $p_{T}(\gamma_{2})>18$ GeV and
$|\eta(\gamma)|<2.4\,,$

  * •

Smooth photon isolation Frixione:1998jh  with $\Delta R_{0}=0.4$ and
$E_{T}^{\rm max}=10~{}{\rm GeV}$ (see ref. Chawdhry:2019bji  for details) ,

  * •

$m_{\gamma\gamma}\geq 90~{}\text{GeV}\,,$

  * •

$\Delta R(\gamma,\gamma)>0.4\,,$

  * •

$p_{T}(\gamma\gamma)>20$ GeV (for lower values resummation effects become
important Becher:2020ugp ) .

No additional jet requirements are imposed. In particular, infrared safety is
ensured by the $p_{T}(\gamma\gamma)$ cut specified above.

##  3 Phenomenological results

Figure 1: Absolute ${p_{T}(\gamma\gamma)}$ (left) and $m(\gamma\gamma)$
(right) differential distributions. Shown are the predictions in LO (green),
NLO (blue), NNLO (red) QCD. The colored bands around the central scales are
from 7-point scale variation. The grey band shows the estimated Monte Carlo
integration error in each bin. The lower panel shows the same distributions
but relative to the NLO central scale prediction.

Figure 2: As in fig. 1 but for the $m(\gamma\gamma)$ distribution subjected to
different ${p_{T}(\gamma\gamma)}$ cuts: ${p_{T}(\gamma\gamma)}>50$ GeV (left),
${p_{T}(\gamma\gamma)}>100$ GeV (center) and ${p_{T}(\gamma\gamma)}>200$ GeV
(right).

In this work we calculate the NNLO QCD corrections to a number of one-
dimensional distributions in the following variables: the transverse momentum
of the photon pair $p_{T}(\gamma\gamma)$, the invariant mass of the two
photons $m(\gamma\gamma)$, the angle between the two photons in the Collins-
Soper frame $\phi_{CS}$, the absolute difference in rapidities of the two
photons $\Delta y(\gamma\gamma)=|y(\gamma_{1})-y(\gamma_{2})|$, the azimuthal
angle between the two photons $\Delta\phi(\gamma\gamma)$ and the absolute
rapidity of the photon pair $|y(\gamma\gamma)|$. We also calculate the NNLO
QCD corrections to the following two-dimensional distributions:
$m(\gamma\gamma)\otimes p_{T}(\gamma\gamma)$ and $\phi_{CS}\otimes
m(\gamma\gamma)$.

We first discuss the ${p_{T}(\gamma\gamma)}$ differential distribution which
is of central interest to this work. The distribution is shown in fig. 1. As
can be seen from this figure, the NLO QCD correction is very significant
relative to the LO one. In particular, the scale uncertainty bands at LO and
NLO do not overlap anywhere. This behavior is easy to understand based on the
properties of inclusive diphoton production through NNLO. Clearly, a reliable
prediction of this observable requires the inclusion of, at least, the NNLO
QCD corrections.

As can be seen from fig. 1 the inclusion of the NNLO corrections has a major
stabilizing impact on the ${p_{T}(\gamma\gamma)}$ distribution. With the
exception of the very low ${p_{T}(\gamma\gamma)}$ region which we will discuss
shortly, the scale uncertainty decreases significantly, by a factor of about
four, relative to NLO. Moreover the NLO and NNLO scale bands now overlap
everywhere. Such a behavior is consistent with the expected stabilization of
the inclusive diphoton production cross-section starting at N3LO.

The low-${p_{T}(\gamma\gamma)}$ behavior of this distribution deserves special
attention. A fixed-order perturbative description would not be adequate for
${p_{T}(\gamma\gamma)}$ below about 20 GeV due to the importance of
resummation effects, but we expect it to be reliable for larger
${p_{T}(\gamma\gamma)}$ values. For this reason one may wonder why the
low-${p_{T}(\gamma\gamma)}$ part of the spectrum shows significantly increased
scale dependence and larger NNLO/NLO $K$-factor. This behavior may be
influenced by resummation effects however we do not expect them to be the
dominant ones. We suspect that the main factor behind it is the loop-induced
contribution $gg\to g{\gamma\gamma}$ which only starts to contribute at NNLO.

The impact of this contribution is shown in fig. 6 as a ratio of the full NNLO
to the NNLO excluding this contribution. As fig. 6 indicates the loop-induced
correction is concentrated at relatively low ${p_{T}(\gamma\gamma)}$ values
and becomes completely negligible for ${p_{T}(\gamma\gamma)}$ values about 200
GeV or larger. If the loop induced correction is excluded, the scale
dependence of the first bin becomes smaller by about a factor of two and the
NNLO/NLO $K$-factor also decreases by a factor of about two.

Our findings indicate that at the level of NNLO QCD corrections, the loop-
induced contribution becomes significant. This contribution can be tamed
further, by including the NLO QCD correction to the loop-induced contribution
(which is a partial N3LO contribution for this process). Such a calculation
requires the two-loop amplitude for the process $gg\to g{\gamma\gamma}$. This
result is not available in the literature but it is certainly within reach
given the number of other five-point two-loop QCD amplitudes that have been
computed.

Overall, the scale uncertainty of the ${p_{T}(\gamma\gamma)}$ spectrum at NNLO
is rather small - about couple of percent for diphoton $p_{T}$’s above 200 GeV
or so. This implies that this observable is well described within perturbative
QCD and can be used in precision analyses of Higgs physics and searches for
resonances decaying to diphoton final states. The dominant uncertainty at
large ${p_{T}(\gamma\gamma)}$ is due to the Monte Carlo integration error of
the calculation itself. It can be further improved albeit at a significant
computational cost. A future update may also include pdf uncertainties,
electroweak effects as well as the NLO correction to the loop-induced process.
Finally, the effects from the the leading color approximation used here may
also need to be improved upon. As can be seen from fig. 6 the approximated
contribution is a rather small 1-2% effect and is much smaller than the MC
error. It can also be improved upon once the complete two-loop amplitudes for
this process become available.

We next turn our attention to the ${m(\gamma\gamma)}$ distribution. It is
shown in figs. 1 and 2. This distribution is significant for any search of
resonances decaying to diphotons at non-zero $p_{T}$. To thoroughly understand
the interplay between ${m(\gamma\gamma)}$ and ${p_{T}(\gamma\gamma)}$ we have
shown the ${m(\gamma\gamma)}$ distribution in several ways. In fig. 1 we show
the ${m(\gamma\gamma)}$ distribution subjected only to our default selection
cuts. In fig. 2 we show the same distribution but with more stringent
${p_{T}(\gamma\gamma)}$ cuts of 50 GeV, 100 GeV and 200 GeV. A summary of the
same result (only the NLO and NNLO are displayed) is shown in fig. 5 (right).
Fig. 5 (center) shows the ${m(\gamma\gamma)}$ distributions for several slices
of ${p_{T}(\gamma\gamma)}$.

From these plots it is clear that the ${m(\gamma\gamma)}$ distribution has a
pattern of higher-order corrections that is similar to
${p_{T}(\gamma\gamma)}$: large NLO/LO corrections and much smaller NNLO/NLO
ones. The size of the NNLO corrections strongly depends on the
${p_{T}(\gamma\gamma)}$ cut and they decrease as the cut increases. For small
values of the ${p_{T}(\gamma\gamma)}$ cut the NNLO and NLO scale bands do not
overlap while they fully overlap for ${p_{T}(\gamma\gamma)}$ cuts above 100
GeV. One may wonder if such a behavior is related to the loop-induced
contribution. In fig. 6 we have shown its effect for all
${p_{T}(\gamma\gamma)}$ cuts considered in this work. From this one can
conclude that indeed the size of the loop-induced correction is consistent
with the non-overlap of the NLO and NNLO scale bands. This means that for
theoretical predictions to be reliable with full NNLO accuracy for
${p_{T}(\gamma\gamma)}$ cuts below 100 GeV or so, the NLO corrections to the
loop induced contributions might need to be included.

In general, the effect of the NNLO correction on the ${m(\gamma\gamma)}$
distribution is a rather flat shift with respect to the NLO one and leads to a
decrease of the scale uncertainty by a factor of about two at low
${p_{T}(\gamma\gamma)}$ and four or more at large ${p_{T}(\gamma\gamma)}$.
Another important source of error is the MC integration one. The effect from
the leading color approximation in the finite remainder is at the percent
level and therefore insignificant.

Figure 3: As in fig. 1 but for the angular distributions in $\phi_{CS}$ (left)
and $\Delta\phi(\gamma\gamma)$ (right).

In fig. 3 we show distributions in the angular variables $\phi_{CS}$ and
$\Delta\phi(\gamma\gamma)$ while in fig. 4 we show the $\Delta
y(\gamma\gamma)$ and $|y(\gamma\gamma)|$ rapidity distributions. The
$\phi_{CS}$ distribution in slices of ${m(\gamma\gamma)}$ is shown in fig. 5
(left). All these distributions have very large NLO/LO $K$-factors. Unlike the
${p_{T}(\gamma\gamma)}$ and ${m(\gamma\gamma)}$ distributions, however, they
also have sizable NNLO corrections which in most bins are outside the NLO
uncertainty bands. This pattern of higher order corrections indicates that for
the scale used in this work, the NLO approximation is inadequate for
describing these distributions.

Based on the above observations one may question the presence of perturbative
stability in these variables. As a first step towards analyzing this we
consider the behavior of the NNLO prediction without the loop-induced
contribution (in the following we refer to it as NNLO-minus-LI). The numerical
impact of the loop-induced contribution for each differential distribution can
be seen in fig. 6 and fig. 7. We observe the following. For the $\phi_{CS}$
distribution the NNLO-minus-LI scale uncertainty band is mostly within the NLO
one or the two bands overlap. This is not the case for the first and last bins
of this distribution, however, the behavior of the $\phi_{CS}$ distribution in
these two bins is strongly affected by the kinematic cuts. The NNLO-minus-LI
band for the $\Delta y(\gamma\gamma)$ distribution overlaps in all bins with
the NLO one. Same can be observed for the case of the $|y(\gamma\gamma)|$
distribution. In fact, the only distribution for which the NNLO-minus-LI and
NLO scale uncertainty bands do not mostly overlap is the
$\Delta\phi(\gamma\gamma)$ one. For this distribution we observe that the
NNLO-minus-LI and NLO scale uncertainty bands overlap for
$\Delta\phi(\gamma\gamma)/\pi>0.6$ while below this value they are not very
far apart, see fig. 7. Given that the NLO/LO $K$-factor in this region is more
than a factor of two it seems that such a non-overlap is not too concerning.

From the above discussion it seems reasonable to conclude that the non-overlap
between NNLO and NLO scale uncertainty bands observed in the angular and
rapidity diphoton distributions is somewhat affected by the loop-induced
contribution. It is therefore plausible to assume that the inclusion of this
contribution’s NLO correction may alleviate this non-overlap. Other factors
that may be affecting this behavior is the choice of scale as well as
resummation effects which are relevant at low ${p_{T}(\gamma\gamma)}$. A
detailed investigation of those is however outside the scope of this work. On
the other hand, as can also be seen from fig. 6, the two-loop finite remainder
has a rather small contribution and, therefore, we do not expect these
distributions to be significantly affected by two-loop subleading color
corrections.

Figure 4: As in fig. 1 but for the following rapidity distributions: $\Delta
y(\gamma\gamma)$ (left) and $|y(\gamma\gamma)|$ (right).

Figure 5: Two-dimensional differential distributions: in $\phi_{CS}\otimes
m(\gamma\gamma)$ (left) and in $m(\gamma\gamma)\otimes p_{T}(\gamma\gamma)$
but shown in two alternative forms and for different choice of bins: in
${p_{T}(\gamma\gamma)}$ bins (center) and with ${p_{T}(\gamma\gamma)}$ cuts
(right). Only the NLO and NNLO central scale predictions and scale variation
bands are shown. Note that the figure to the right shows the same results that
already appear in figs. 1,2.

Figure 6: Comparison of various approximations to the NNLO differential
predictions: the complete NNLO prediction (red), NNLO excluding the loop-
induced contribution (dark brown) and NNLO excluding the finite remainder
${\cal H}^{(2)}(s_{12})$ defined in eq. (1) (green). The grey band shows the
MC integration error of the complete NNLO prediction while the red band shows
its scale variation.

Figure 7: Comparison of NLO (blue) scale bands to NNLO QCD excluding the loop-
induced contribution (dark brown).

##  4 Conclusions

In this work we calculate the NNLO QCD corrections to the process
$pp\to{\gamma\gamma}+{\rm jet}$. This process is the main background to
high-$p_{T}$ Higgs boson production decaying to two photons. The main result
of this work is the calculation of the diphoton $p_{T}$ spectrum with NNLO
accuracy. The NNLO correction to this variable is important and it brings the
uncertainty from un-calculated higher-order corrections down to a couple of
percent at intermediate and large values of ${p_{T}(\gamma\gamma)}$. Overall,
the quality of the theoretical prediction for this distribution is very high
and it appears to be under good theoretical control. The same conclusion
applies for double differential distributions in ${p_{T}(\gamma\gamma)}$ and
${m(\gamma\gamma)}$.

We have suggested various possible avenues for further improving the quality
of the theoretical predictions in this process. They include the calculations
of the partial N3LO corrections due to loop induced processes which can be
calculated with the help of NLO technology. The only missing ingredient for
such a calculation is the two-loop amplitude for the process $gg\to
g{\gamma\gamma}$ whose calculation is within reach. A more extensive study of
possible scale choices for this process might also be beneficial given the
very high precision reached in the ${p_{T}(\gamma\gamma)}$ distribution.
Merging our fixed-order calculations with resummed calculations will allow for
a quality description of the ${p_{T}(\gamma\gamma)}$ spectrum from very high
down to very low values of ${p_{T}(\gamma\gamma)}$.

We conclude by stressing that the quality of the theoretical description
achieved for this process is high which makes it possible to use it in
background estimates for Higgs boson studies and related searches as well as
in dedicated measurements of diphoton production.

Note Added: After the completion of the current work, ref. Agarwal:2021vdh
appeared. It provides the subleading-colour expressions for the two-loop
amplitudes for this process. We will include them in a future update of this
work.

###### Acknowledgements.

The work of M.C. was supported by the Deutsche Forschungsgemeinschaft under
grant 396021762 - TRR 257. The research of A.M. and R.P. has received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 Research and Innovation Programme (grant agreement no. 683211). A.M. was
also supported by the UK STFC grants ST/L002760/1 and ST/K004883/1. The
research of H.C. is supported by the ERC Starting Grant 804394 HipQCD. A.M.
acknowledges the use of the DiRAC Cumulus HPC facility under Grant No.
PPSP226.

## References

  * (1) S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini, Phys. Rev. Lett. 108, 072001 (2012) [erratum: Phys. Rev. Lett. 117, no.8, 089901 (2016)] [arXiv:1110.2375 [hep-ph]].
  * (2) J. M. Campbell, R. K. Ellis, Y. Li and C. Williams, JHEP 07, 148 (2016) [arXiv:1603.02663 [hep-ph]].
  * (3) S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini, JHEP 04, 142 (2018) [arXiv:1802.02095 [hep-ph]].
  * (4) T. Gehrmann, N. Glover, A. Huss and J. Whitehead, JHEP 01, 108 (2021) [arXiv:2009.11310 [hep-ph]].
  * (5) S. Alioli, A. Broggio, A. Gavardi, S. Kallweit, M. A. Lim, R. Nagar, D. Napoletano and L. Rottoli, JHEP 04, 041 (2021) [arXiv:2010.10498 [hep-ph]].
  * (6) V. Del Duca, F. Maltoni, Z. Nagy and Z. Trocsanyi, JHEP 04, 059 (2003) [arXiv:hep-ph/0303012 [hep-ph]].
  * (7) M. Chiesa, N. Greiner, M. Schönherr and F. Tramontano, JHEP 10, 181 (2017) [arXiv:1706.09022 [hep-ph]].
  * (8) T. Gehrmann, N. Greiner and G. Heinrich, JHEP 06, 058 (2013) [erratum: JHEP 06, 076 (2014)] [arXiv:1303.0824 [hep-ph]].
  * (9) L. Cieri, F. Coradeschi and D. de Florian, JHEP 06, 185 (2015) [arXiv:1505.03162 [hep-ph]].
  * (10) F. Coradeschi and T. Cridge, Comput. Phys. Commun. 238, 262-294 (2019) [arXiv:1711.02083 [hep-ph]].
  * (11) M. A. Ebert and F. J. Tackmann, JHEP 02, 110 (2017) [arXiv:1611.08610 [hep-ph]].
  * (12) T. Becher and T. Neumann, JHEP 03, 199 (2021) [arXiv:2009.11437 [hep-ph]].
  * (13) L. Cieri, F. Coradeschi, D. de Florian and N. Fidanza, Phys. Rev. D 96, no.5, 054003 (2017) [arXiv:1706.07331 [hep-ph]].
  * (14) S. Amoroso, P. Azzurri, J. Bendavid, E. Bothmann, D. Britzger, H. Brooks, A. Buckley, M. Calvetti, X. Chen and M. Chiesa, et al. [arXiv:2003.01700 [hep-ph]].
  * (15) H. A. Chawdhry, M. L. Czakon, A. Mitov and R. Poncelet, JHEP 02, 057 (2020) [arXiv:1911.00479 [hep-ph]].
  * (16) S. Kallweit, V. Sotnikov and M. Wiesemann, Phys. Lett. B 812, 136013 (2021) [arXiv:2010.04681 [hep-ph]].
  * (17) F. Caola, A. Von Manteuffel and L. Tancredi, Phys. Rev. Lett. 126, no.11, 112004 (2021) [arXiv:2011.13946 [hep-ph]].
  * (18) A. M. Sirunyan et al. [CMS], JHEP 09, 046 (2018) [arXiv:1806.04771 [hep-ex]].
  * (19) G. Aad et al. [ATLAS], [arXiv:2104.13240 [hep-ex]].
  * (20) J. M. Lindert, K. Kudashkin, K. Melnikov and C. Wever, Phys. Lett. B 782, 210-214 (2018) [arXiv:1801.08226 [hep-ph]].
  * (21) J. C. Collins and D. E. Soper, Phys. Rev. D 16, 2219 (1977)
  * (22) G. Aad et al. [ATLAS], Eur. Phys. J. C 75, no.10, 476 (2015) [erratum: Eur. Phys. J. C 76, no.3, 152 (2016)] [arXiv:1506.05669 [hep-ex]].
  * (23) K. Becker, F. Caola, A. Massironi, B. Mistlberger, P. F. Monni, X. Chen, S. Frixione, T. Gehrmann, N. Glover and K. Hamilton, et al. [arXiv:2005.07762 [hep-ph]].
  * (24) M. Czakon, Phys. Lett. B 693, 259 (2010) [arXiv:1005.0274 [hep-ph]].
  * (25) M. Czakon, Nucl. Phys. B 849, 250 (2011) [arXiv:1101.0642 [hep-ph]].
  * (26) M. Czakon and D. Heymes, Nucl. Phys. B 890, 152 (2014) [arXiv:1408.2500 [hep-ph]].
  * (27) M. Czakon, P. Fiedler and A. Mitov, Phys. Rev. Lett. 115, no. 5, 052001 (2015) [arXiv:1411.3007 [hep-ph]].
  * (28) M. Czakon, D. Heymes and A. Mitov, Phys. Rev. Lett. 116, no. 8, 082003 (2016) [arXiv:1511.00549 [hep-ph]].
  * (29) M. Czakon, P. Fiedler, D. Heymes and A. Mitov, JHEP 1605, 034 (2016) [arXiv:1601.05375 [hep-ph]].
  * (30) M. Czakon, D. Heymes and A. Mitov, JHEP 1704, 071 (2017) [arXiv:1606.03350 [hep-ph]].
  * (31) A. Behring, M. Czakon, A. Mitov, A. S. Papanastasiou and R. Poncelet, Phys. Rev. Lett. 123, no. 8, 082001 (2019) [arXiv:1901.05407 [hep-ph]].
  * (32) M. Czakon, A. van Hameren, A. Mitov and R. Poncelet, JHEP 10, 262 (2019) [arXiv:1907.12911 [hep-ph]].
  * (33) M. Czakon, A. Mitov, M. Pellen and R. Poncelet, [arXiv:2011.01011 [hep-ph]].
  * (34) M. L. Czakon, T. Generet, A. Mitov and R. Poncelet, [arXiv:2102.08267 [hep-ph]].
  * (35) R. Poncelet and A. Popescu, [arXiv:2102.13583 [hep-ph]].
  * (36) https://bitbucket.org/hameren/avhlib.
  * (37) M. Bury and A. van Hameren, Comput. Phys. Commun. 196, 592 (2015) [arXiv:1503.08612 [hep-ph]].
  * (38) F. Cascioli, P. Maierhofer and S. Pozzorini, Phys. Rev. Lett. 108, 111601 (2012) [arXiv:1111.5206 [hep-ph]].
  * (39) F. Buccioni, J. N. Lang, J. M. Lindert, P. Maierhofer, S. Pozzorini, H. Zhang and M. F. Zoller, arXiv:1907.13071 [hep-ph].
  * (40) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, [arXiv:2103.04319 [hep-ph]].
  * (41) H. A. Chawdhry, M. A. Lim and A. Mitov, Phys. Rev. D 99, no.7, 076011 (2019) [arXiv:1805.09182 [hep-ph]].
  * (42) D. Chicherin and V. Sotnikov, JHEP 12, 167 (2020) [arXiv:2009.07803 [hep-ph]].
  * (43) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, [arXiv:2012.13553 [hep-ph]].
  * (44) B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi, JHEP 04, 201 (2021) [arXiv:2102.01820 [hep-ph]].
  * (45) R. D. Ball et al. [NNPDF], Eur. Phys. J. C 77, no.10, 663 (2017) [arXiv:1706.00428 [hep-ph]].
  * (46) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr and G. Watt, Eur. Phys. J. C 75, 132 (2015) [arXiv:1412.7420 [hep-ph]].
  * (47) A. M. Sirunyan et al. [CMS], JHEP 11, 109 (2019) [arXiv:1908.08500 [hep-ex]].
  * (48) S. Frixione, Phys. Lett. B 429, 369-374 (1998) [arXiv:hep-ph/9801442 [hep-ph]].
  * (49) B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi, [arXiv:2105.04585 [hep-ph]].


# Dimensional reduction of gradient-like stochastic systems with
multiplicative noise via Fokker-Planck diffusion maps

Andrew Baumgartner<EMAIL_ADDRESS>Sui Huang
<EMAIL_ADDRESS>Jennifer Hadlock<EMAIL_ADDRESS>Cory Funk<EMAIL_ADDRESS>Institute for Systems Biology
401 Terry Ave, Seattle, WA 98109

###### Abstract

Dimensional reduction techniques have long been used to visualize the
structure and geometry of high dimensional data. However, most widely used
techniques are difficult to interpret due to nonlinearities and opaque
optimization processes. Here we present a specific graph based construction
for dimensionally reducing continuous stochastic systems with multiplicative
noise moving under the influence of a potential. To achieve this, we present a
specific graph construction which generates the Fokker-Planck equation of the
stochastic system in the continuum limit. The eigenvectors and eigenvalues of
the normalized graph Laplacian are used as a basis for the dimensional
reduction and yield a low dimensional representation of the dynamics which can
be used for downstream analysis such as spectral clustering. We focus on the
use case of single cell RNA sequencing data and show how current diffusion map
implementations popular in the single cell literature fit into this framework.

††preprint: APS/123-QED

##  I Introduction

The problem of dimensional reduction of single cell RNA sequencing data (sc-
RNAseq) is a difficult one. The analysis of these high dimensional, discrete
systems usually involves a panoply of preprocessing tools that map a vector
$x\in\mathbb{N}_{0}^{m}$ (where $m$ is the total number of genes in the
sample) to some dense subset of a submanifold embedded in the positive orthant
of $\mathbb{R}^{m}$. A suite of tools is then used to “learn" this data
manifold and utilize geometric information to make claims about various
biological processes, typically using $k$ nearest neighbor (kNN) graphs and/or
clustering. The validity of this procedure has, and should be, questioned [1,
2] but nonetheless is motivated by the existence of a continuum limit of the
more fundamental chemical master equations. In this work, we introduce a new
twist on an old method that is meant to act as a bridge between this “top-
down" approach of manifold learning and the “bottom-up" approach of chemical
master equations via the spectra of the Fokker-Planck equation.

Specifically, we introduce a specific kNN graph whose normalized graph
Laplacian generates the backward Fokker-Planck equation of generic, continuous
model of single cell dynamics. The diffusion maps [3, 4] associated with this
graph are then used for spectral clustering to identify cell types and
subtypes. A general analysis of graph Laplacians was carried out in [5], where
a correspondence was established between general, isotropic stochastic
differential equations and graph Laplacians via drift-diffusion equations. By
finding the correct normalization of the graph Laplacian, we can model the
dynamics of a single cell moving through a potential landscape, a popular
model for cellular development initially due to Waddington [6, 7, 8].

Diffusion maps have been used by the sc-RNAseq community since at least 2015
[9]. They have primarily been used low dimensional representation of the
geometry of high dimensional data. This approach comes with many benefits,
including being insensitive to the sampling density of the data. A more recent
implementation introduced in [10], and integrated into popular single cell
packages such as scanpy [11] and Seurat [12], introduced local bandwidths and
weights to the graph normalization for better visualization, but no longer
carries a direct interpretation to the geometry of the data manifold, as we
will see below.

In section II we will summarize prior results and explain how we used them in
our specific graph construction. In section III we derive a general formula
for the potential of the system in terms of nearest neighbors, and show how
this simplifies when one approximates the strength of the noise as a linear
function of the distance to the $k^{th}$ nearest neighbor. We introduce our
choice of kernel to obtain the Fokker-Planck diffusion maps in section IV and
apply it to numerous simulated and real biological systems in section V. We
conclude with a discussion in section VI. For simplicity, we assume
$x,y\in\mathbb{R}^{m}$ and omit bold faced type for vectors. Additionally, we
assume the existence of an optimal $k$ such that the kNN graph yields an
appropriate and accurate discretizaiton of the data manifold. In reality, many
$k$’s should be used to track any results across many scales.

##  II The circle of diffusion processes

The following relies on a result concerning general diffusion processes. These
are defined by the solution to the following stochastic differential equation

$dx=\mu(x)dt+\sigma(x)dW$ (1)

where $dW$ is a Wiener process. The arbitrary functions $\mu(x)$ and
$\sigma(x)$ are the drift and diffusion terms and are given by vectors and
matrices, respectively, which depend on the state. The generator $\mathcal{G}$
of this diffusion processes is given by

$\mathcal{G}f=\frac{1}{2}\sum_{i,j,k}\sigma(x)_{ik}\sigma(x)_{kj}^{T}\frac{\partial^{2}f}{\partial
x_{i}\partial x_{j}}+\sum_{i}\mu_{i}(x)\frac{\partial f}{\partial x_{i}}.$ (2)

Thus, for every diffusion process there is a corresponding elliptic equation
which generates the process. What follows is a result from [5] that
establishes a connection between eq. (2) and normalized graph Laplacians. The
graph construction can be viewed as a discrete approximation to eq. (2).

We assume we have a translationally and rotationally invariant kernel
$K_{0}(|x-y|)$ with compact support. We then introduce a weight function
$w_{x}(y)\geq 0$ and a bandwidth function $r_{x}(y)>0$ which admit a Taylor
expansion about $y=x$ with bounded error, denoted by $w^{(n)}_{x}(y)$ and
$r^{(n)}_{x}(y)$. We then rescale the kernel:

$K(x,y)=w^{(n)}_{x}(y)K_{0}\left(\frac{|x-y|}{h_{n}r^{(n)}_{x}(y)}\right)$ (3)

where $h_{n}>0$ is a bandwidth scaling which goes to 0 in the continuum limit.
The following theorem can be found in [5]:

######  Theorem 1.

[5] Assume the standard assumptions for the diffusion approximation to hold
eventually with probability 1. If the bandwidth scaling $h_{n}$ satisfies
$h_{n}\to 0$ and $nh_{n}^{m+2}/log(n)\to\infty$ then for graphs constructed
using the kernel

$K(x,y)=w^{(n)}_{x}(y)K_{0}\left(\frac{|x-y|}{h_{n}r^{(n)}_{x}(y)}\right)$

there exists a constant $Z_{K_{0},m}>0$ depending only on the base kernel
$K_{0}$ and the dimension of the data manifold $m$ such that for
$c_{n}=Z_{K_{0},m}/h^{2}$,

$-c_{n}L^{(n)}_{rw}f\to Af$

where A is the infinitesimal generator of a diffusion process with the
following drift and diffusion terms given in normal coordinates

$\displaystyle\begin{split}\mu(x)&=r_{x}(x)^{2}\left(\nabla_{y}\ln
p(y)+\nabla_{y}\ln w_{x}(y)\right.\\\
&\left.\quad+(m+2)\nabla_{y}r_{x}(y)\vphantom{\nabla_{y}\ln
p(y)}\,\right)_{|y=x}\end{split}$ (4) $\displaystyle\sigma(x)\sigma(x)^{T}$
$\displaystyle=r_{x}(x)^{2}I$ (5)

where all derivatives are taken with respect to $y$ (in appropriate
coordinates) and evaluated at $y=x$, $p$ is the local sampling density at $x$
and $I$ is the $m\times m$ identity matrix.

This allows us to approximate any isotropic equation of the form (2) via a
graph with edge weights given by equation (3) with a clever choice of
bandwidth and weight functions.

The original authors of the diffusion maps in [4, 3, 13] introduced the
following normalization:

$w_{x}(y)=p(x)^{-\alpha}p(y)^{-\alpha},\,\,\,r_{x}(y)=1$ (6)

which leads to

$\mathcal{G}f=\frac{1}{2}\nabla^{2}f+(1-\alpha)\nabla(\ln p)\cdot\nabla f.$
(7)

When choosing $\alpha=1$ we recover the Laplacian of the manifold $\nabla^{2}$
and the $p$ dependence drops out. This normalization was used in [9] for
single cell data. A locally scaled version was introduced in [10] that induces
a spurious drift term.

##  III Noise induced drift and the equilbrium distribtuion

The work presented in [4, 3, 13] offers a way to represent gradient-like
stochastic systems with additive noise via diffusion maps. This is done by
taking $\alpha=1/2$ in equation (6) and assuming the local density is of
Boltzmann form, i.e. $p=e^{-U}$ for some potential $U$. This results in the
following drift-diffusion equation

$-\frac{\partial f}{\partial t}=\nabla^{2}f+\nabla U\cdot\nabla f$ (8)

which is the backwards Fokker-Planck equation of the system

$dx=-\nabla Udt+\sqrt{2}dW.$

Solutions to (8) can be approximated by an eigenvalue expansion using the
eigenvalues and eigenvectors of the graph Laplacian. The diffusion map gives a
low dimensional representation of the system allowing for the identification
of local minima and the paths connecting them. Moreover, this is the optimal
dimensional reduction technique with respect to the standard error between the
truncated eigenvalue expansion and the true solution [13]. This choice of
local density is equivalent to saying that the system is in local
“information-thermodynamic" equilibrium with respect to the potential $U$.

Extending this procedure to systems with multiplicative noise is not so
straightforward due to an ambiguity in the integration procedure of the
associated stochastic differential equation eq. (1). The ambiguity is most
easily seen through the corresponding Langevin equation

$\frac{dx}{dt}=\mu(x)+\sigma(x)\eta(t)$ (9)

where $\eta(t)\sim dW/dt$. This is the source of the Ito vs. Stratonovich
controversy. The choice does not matter in systems with additive noise as well
as systems with internal noise. The only situation in which the choice has
physical consequences is in systems with external noise, whereby the
Stratonovich interpretation is the correct one [14]. scRNAseq data is unique
in that it has both internal and external noise. Common preprocessing
pipelines are meant to eliminate technical (external) noise and only keep
biological (internal) noise, although their efficacy is widely debated (See
[1, 2] for alternative approaches).

In the Stratonovich picture, the Fokker-Planck equation with $\mu=-\nabla U$
becomes

$\frac{\partial p}{\partial t}=-\nabla\cdot\left(\left(-\nabla
U+\frac{1}{2}\nabla\sigma^{2}\right)p-\frac{1}{2}\nabla\left(\sigma^{2}p\right)\right)$
(10)

which leads to an equilibrium distribution of the form

$p=\frac{\mathcal{N}}{\sigma^{2}}\exp\left\\{\int^{x}\frac{2\left(-\nabla
U+\frac{1}{2}\nabla\sigma^{2}\right)}{\sigma^{2}}dz\right\\}$ (11)

where $\mathcal{N}$ is a normalization factor. The term
$\frac{1}{2}\nabla\sigma^{2}$ is a “noise induced drift" term. The
consequences of this term on the Waddington landscape was beautifully
illustrated in [15]. There, the authors show that this term flattens out
existing minima of the potential and changes their location [15]. We could
proceed by defining an effective potential as $U_{eff}=-\ln p$, but this
obscures the physics of the underlying stochastic process. Luckily, when we
take the noise term to be proportional to the $k^{th}$ nearest neighbor of
$x$, things simplify and we can approximate the potential using only knowledge
of the nearest neighbors.

To achieve this we utilize a trick from Ting et. al [5] and use the kNN
density estimate [16] as our local density. In $\mathbb{R}^{m}$ we have

$p_{kNN}(x)=\frac{k}{N}\frac{1}{\Omega_{m}\rho(x)^{m}}.$ (12)

where $\Omega_{m}$ is the volume of the unit $m$-ball, and $\rho(x)$ is the
distance to the $k^{th}$ nearest neighbor. If we equate this to the local
equilibrium density with the noise induced drift term, we can solve for the
potential in terms of the $k$-nearest neighbor field and the strength of the
noise. To do so, we start with

$\frac{\mathcal{N}}{\sigma^{2}}\exp\left\\{\int^{x}\frac{2\left(-\nabla
U+\frac{1}{2}\nabla\sigma^{2}\right)}{\sigma^{2}}\cdot
dz\right\\}\approx\frac{C}{\rho(x)^{m}}$ (13)

where $C=k/(N\Omega_{m})$ and the integral is shorthand for a line integral
measured from the zero point of the potential, i.e. from $x_{0}$ to $x$ such
that $U(x_{0})=0$. After a little algebra which we will not recapitulate here,
we wind up with the general relation

$U(x)=\int^{x}\sigma^{2}\nabla\ln\left(\rho^{m/2}\right)dz.$ (14)

If we then take $\sigma=\alpha\rho$ for $\alpha\in(0,1]$, this simplifies
further to

$U(x)=\frac{m\alpha^{2}}{4}\rho(x)^{2}.$ (15)

Thus, in order to approximate the SDE

$dX=-\nabla Udt+\sigma(x)dW$ (16)

we need a normalization which gives

$dX=-\frac{m\alpha^{2}}{4}\nabla\rho(x)^{2}dt+\alpha\rho(x)dW$ (17)

up to some scale factor.

##  IV Fokker-Planck diffusion maps

There are many different choices of kernel that will converge to the same
equation in the continuum limit. The results of Ting et al. [5] give an
infinite family of graph normalizaitons which lead to the same drift-diffusion
equation. This means that the identification of long-time macroscopic dynamics
will be the same for these different graphs, while the short-time, microscopic
dynamics will be different. For downstream applications this means that things
like boundaries between detected clusters on the graph will be different
between graphs with different kernels, while the macroscopic, coarse grained
clusters will all converge. This suggests shying away from graph-based
clustering in favor of spectral clustering.

We use a simple symmetric Gaussian kernel constructed as follows: start with a
base kernel of the form

$K(x,z)=e^{-(x-z)^{2}/\rho(x)^{2}}$ (18)

and symmetrize it. For instance, one can do the following:

$\displaystyle K_{sym}(x,y)$ $\displaystyle=\int K(x,z)K(y,z)dz$ (19)
$\displaystyle=\left(\frac{\pi\rho(x)^{2}\rho(y)^{2}}{\rho(x)^{2}+\rho(y)^{2}}\right)^{m/2}e^{-(x-y)^{2}/(\rho(x)^{2}+\rho(y)^{2})}.$
(20)

In discrete notation, this is equivalent to

$K_{sym}=KK^{T}.$ (21)

This normalization leads to the following stochastic differential equation:

$dX=(1-\frac{m}{4})\nabla\rho(x)^{2}+\rho(x)dW$ (22)

which for $m\gg 4$ agrees with (17). The downside to this symmetrization is
that it effectively turns the nearest neighbor graph into a shared nearest
neighbor graph, and so our graph may not be fully connected. To see this, note
that the initial nearest neighbor kernel can be extended to all data points by
including a Heaviside step function as in [5]:

$K_{full}(x,z)=K(x,z)\Theta(|z-x|<\rho(x)).$ (23)

When we perform the integral in eq. to symmetrize (19) we obtain

$\begin{split}&K_{sym}(x,y)=\\\ &\int
K(x,z)K(y,z)\Theta(|z-x|<\rho(x))\Theta(|z-y|<\rho(y))dz.\end{split}$ (24)

which is zero unless the neighborhoods of $x$ and $y$ are overlapping. Another
option which is guaranteed to return a connected graph is to take $K$ as:

$K_{sym}(x,y)=K(x,y)+K(y,x)$ (25)

or, in matrix notation

$K_{sym}=K+K^{T}.$ (26)

The linearity of the graph Laplacian and Taylor expansions used in [5] allows
us to add the contributions to eq. (4) from both terms separately. Doing so
yields the following SDE:

$dX=\left(1-\frac{m}{2}\right)\nabla\rho^{2}dt+\sqrt{2}\rho(x)dW.$ (27)

We can redefine $\rho\to\rho/\sqrt{2}$ to instead obtain

$dX=\frac{1}{2}\left(1-\frac{m}{2}\right)\nabla\rho^{2}dt+\rho(x)dW.$ (28)

which for $m\gg 2$ reduces to the correct result. The eigenvectors and
eigenvalues of the normalized graph Laplacians of these graphs can be used as
a basis for dimensional reduction just like the diffusion maps of [4]. They
can be used for visualizing the potential, as a basis for spectral clustering,
or for finding low dimensional dynamical systems which approximate the
original system as in [13, 17].

The local bandwidth and weight functions for diffusion maps introduced in [10]
can also be interpreted in this framework, although the initial justification
for their construction is unclear. The authors model a cell in gene expression
space as a Gaussian blob using the following:

$Y_{x}(x^{\prime})=\left(\frac{2}{\pi\sigma_{x}^{2}}\right)^{1/4}e^{-\frac{|x-x^{\prime}|^{2}}{\sigma_{x}^{2}}}.$
(29)

We believe this normalization is unjustified, since for
$x,x^{\prime}\in\mathbb{R}^{m}$ the integral over $x^{\prime}$ is

$\begin{split}\int
d^{m}x^{\prime}&\left(\frac{2}{\pi\sigma_{x}^{2}}\right)^{1/4}e^{-\frac{|x-x^{\prime}|^{2}}{\sigma_{x}^{2}}}\\\
&=\left(\frac{2}{\pi\sigma_{x}^{2}}\right)^{1/4}\left(\pi\sigma_{x}^{2}\right)^{m/2}\\\
&=2^{\frac{1}{4}}\pi^{\frac{2m-1}{4}}\sigma_{x}^{m-\frac{1}{2}}\end{split}$
(30)

which does not have any clear biological or biophysical interpretation.
Moreover, when constructing the graph kernel as as the overlap of these
Gaussian blobs, the authors of [10] use

$\begin{split}K(x,y)&=\int dzY_{x}(z)Y_{y}(z)\\\
&=\sqrt{\frac{2\sigma_{x}\sigma_{y}}{\sigma_{x}^{2}+\sigma_{y}^{2}}}e^{-\frac{|x-y|^{2}}{2(\sigma_{x}^{2}+\sigma_{y}^{2})}}\end{split}$
(31)

which is incorrect for $x,y,z\in\mathbb{R}^{m}$. The correct result should be

$\begin{split}K(x,y)&=\int dzY_{x}(z)Y_{y}(z)\\\
&=\frac{4^{1/4}\pi^{(m-1)/2}\left(\sigma_{y}^{2}\sigma_{y}^{2}\right)^{(m-2)/2}}{(\sigma_{x}^{2}+\sigma_{y}^{2})^{m/2}}e^{-\frac{|x-y|^{2}}{(\sigma_{x}^{2}+\sigma_{y}^{2})}}.\end{split}$
(32)

which reduces to their result if $m=1$. Following [4], the authors then define
$Z(x)=\sum_{y}K(x,y)$ so the final graph kernel is

$W(x,y)=Z(x)^{-1}Z(y)^{-1}K(x,y).$ (33)

The definition of $Z(x)$ employed here is the kernel density estimate of
$p(x)$ using the locally scaled kernel This follows from [4].

Biological interpretation and incorrect integrals aside, the normalization of
[10] et al. can still be interpreted within the framework of Theorem 1. A
little algebra shows that these weight and bandwidth functions yield a graph
which models the following SDE:

$dX=(m+\frac{3}{2})\nabla\sigma_{x}^{2}dt+2\sigma_{x}dW.$ (34)

Upon taking $\sigma=\rho(x)/2$ we obtain

$dX=\frac{1}{4}(m+\frac{3}{2})\nabla\rho^{2}dt+\rho(x)dW.$ (35)

Crucially the drift terms in eqs. (35) and (28) differ by a minus sign. In the
Waddington landscape view of cell development, cell types and subtypes
correspond to local minima of the potential. In analogy to real, physical
potentials, the difference in minus signs determines whether the force
experienced by a given cell will drive them towards a minimum or a maximum.
These corresponds to situations whereby $\mu(x)=-\nabla U$ versus
$\mu(x)=\nabla U$ respectively. This has the consequence that the SDE given in
(28) will drive cells towards more dense regions (those with $\nabla\rho<0$),
while the SDE in (35) will drive cells towards less dense regions (those with
$\nabla\rho>0$). The latter is good for visualization, since it will “spread
out" the data, so to speak, while the former is more appropriate for
clustering and identification of cell types, since it drives cells towards
their closest minimum.

##  V Illustrations and applications

We now turn to specific examples to illustrate our proposed diffusion map. We
first consider the dynamics of two multidimensional Ornstein-Uhlenbeck
processes–one with additive noise and one with linear, multiplicative noise–
and show how our Fokker-Planck diffusion maps can reflect the dynamics of a
single trajectory in the system. Next, we turn to a multidimensional double
well potential with orthogonal noise directions and demonstrate how an
ensemble of trajectories can induce a clear separation of minima in the
resulting diffusion map. Finally, we look at the PBMC3k set and compare to the
implementation of Haghverdi et. al. [10] to demonstrate the differences
derived above.

###  V.1 Ornstein-Uhlenbeck processes

(a) 2D projection of our 10D Ornstein-Uhlenbeck (O-U) trajectory with additive
noise from t=0 (blue) to t=1 (red)

(b) The first three diffusion coordinates of the O-U trajectory.

Figure 1: A single trajectory for an Ornstein-Uhlenbeck process in 10D with
additive noise and the corresponding diffusion map. The diffusion coordinates
correspond to products of Hermite polynomials as shown in [4]. The matrices
$\alpha$ and $\beta$ have random entries sampled from a uniform distribution
over $[0,1]$

(a) 2D projection of our 10D Ornstein-Uhlenbeck (O-U) trajectory with linear,
multiplicative noise from t=0 (blue) to t=1 (red)

(b) The first three diffusion coordinates of the O-U trajectory.

Figure 2: A single trajectory for an Ornstein-Uhlenbeck process in 10D with
linear, multiplicative noise with $\omega=.05$ and $\beta=.001$.

(a) Ensemble of trajectories of the 10D double well. The potential is
localized in 2D and orthogonal directions are subject to pure diffusion. The
colorbar denotes the value of the potential.

(b) The first three diffusion coordinates of our 10D ensemble of trajectories,
showing separation between the two minima of the potential.

Figure 3: Our double well in eq. (41) with $d=2$, $a=1$, $b=c=4$ and
$\sigma(x)=1-.2x^{2}$. Trajectories are colored by the minima they equilibrium
to. Note the presence of trajectories which start in one minima and end up in
another in both (a) and (b).

(a) The first three diffusion coordinates using the Fokker-Planck diffusion
maps.

(b) The first three diffusion coordinates using the scanpy’s diffusion maps.

Figure 4: The first three diffusion coordinates using the two different
normalizations for diffusion maps discuss. Note the FP diffusion maps tend to
clump similar cells together, while the scanpt diffusion maps spread them
apart. This is a consequence of the relative sign in the force term of the
corresponding SDE.

(a) Confusion matrix for hierarchical clustering on FP diffusion maps

(b) Confusion matrix for hierarchical clustering on scanpy diffusion maps

Figure 5: Confusion matrices for hierarchical clustering on the PBMC3k data
set. Clustering was performed on the first 6 nontrivial diffusion coordinates.
This corresponds to the number of eigenvalues before the largest spectral gap.

The Ornstein-Uhleneck (O-U) process with additive noise is the rare example of
a stochastic system where both the stochastic differential equation and the
diffusion map can be solved exactly making it a perfect sanity check for our
new normalization. Specifically, we have

$dx=-\alpha(x-\gamma)+\beta dW$ (36)

where, in general, $\alpha$ and $\beta$ are matrices. This system was treated
in [4] where it was shown that the diffusion maps can be written in terms of
(products of) Hermite polynomials. Figure 1 illustrates these diffusion maps
on a single trajectory in 10D.

Next we look at an O-U process linear multiplicative noise. Specifically, we
consider

$dx=Fxdt+GxdW.$ (37)

This can be solved if the matrices $F$ and $G$ commute, giving

$x(t)=e^{(F-\frac{1}{2}G^{2})t+GW(t)}x(0).$ (38)

In particular, we use the following matrices

$\displaystyle F$
$\displaystyle=\text{antidiag}(\omega,-\omega,\omega,\dots,-\omega)$ (39)
$\displaystyle G$ $\displaystyle=\text{diag}(\beta,\dots,\beta).$ (40)

The results are shown in figure 2.

###  V.2 Multidimensional double well

Next we consider a double well potential in multiple dimensions subject to
multiplicative noise. We use the following form of the potential

$U(x)=x_{0}^{2}(ax_{0}^{2}-b)+\frac{b^{2}}{4a}+\frac{1}{2}\sum_{i=0}^{d-2}c_{i}(x_{i}-x_{i+1})^{2}$
(41)

for some $d<D$ where $D$ is the dimensional of the system. This localizes the
potential into a $d$-dimensional subspace of the full phase space with pure
diffusion in the orthogonal directions. This particular form is meant to mimic
a system in which the data lives on some low dimensional subspace of the
entire phase space, with noise in orthogonal directions. We simulate multiple
trajectories in this space and run our diffusion maps on the resulting data
set. Our results are shown in figure 3, where we choose

$\sigma(x)=1-.2x^{2}.$ (42)

We see a clear separation of trajectories falling into their distinct minima
with interpolation between the two representing trajectories which jump the
potential barrier, suggesting the diffusion maps recapitulate the dynamics of
the entire ensemble.

###  V.3 Single Cell RNAseq–PBMC3k

Finally we look at our primary use case–dimensional reduction and subsequent
clustering of scRNAseq data. We compare our results to the current
implementation of diffusion maps in scanpy [11], with the method laid out in
[10]. The diffusion map for both of these is shown in figure 4. We see that
our implementation tends to group cells closer together, while the scanpy
implementation spreads them out more. This is due to the relative minus sign
in the force terms.

We then use these diffusion maps as a basis for hierarchical clustering. We
restrict our diffusion maps to be above the first spectral gap (corresponding
to the first 6 nontrivial diffusion coordinates), perform hierarchical
clustering and compare our clusters to the labels of the PBMC3k set.
Specifically, we used the processed data available through the scanpy API. Our
results are shown in figure 5. Note that overall the two methods are
comparable, with only a $2\%$ difference in the balanced accuracy. However, we
find that our method is preferable when cells are more difficult to
distinguish, as in the case of CD4 and CD8 T-cells, and the CD14 and FCGR3A+
monocytes. We believe this is due to the sign of our potential term which, as
discussed above, drive cells to the nearest minima, as opposed to the scanpy
implementation which drives cells away from the minima.

##  VI Discussion

We presented a new construction for a normalized graph Laplacian which
generates the backwards Fokker-Planck equation associated to a stochastic
differential equation of the form of eq. 1 with $\mu(x)=-\nabla U$ and show
how the associated diffusion maps (Fokker-Planck diffusion maps) can be used
to represent the dynamics. To do so we first assumed the system is in local
equilibrium such that the local density about each point is equivalent to the
equilibrium density. By equating the equilibrium density to the kNN density
estimate, we derived a simple relationship between the potential, the noise,
and the distance to the $k^{th}$ nearest neighbor. If we further assume the
noise is equal to a linear function of the distance to the $k^{th}$ nearest
neighbor, we can simplify the equation drastically, leading to a simple graph
kernel for our system. We then show the differences between our normalization
and the most commonly used normalization for diffusion maps introduced in [10]
and used in popular single cell packages like scanpy [11]. Finally, we
illustrated our graph on single trajectories of stochastic systems with
additive and multiplicative noise to show how the diffusion map can capture
the dynamics [17], and on an ensemble trajectories in a double well potential
localized to a 2-dimensional subspace of the a full 10 dimensional system. We
then showed improved hierarchical clustering of the the PBMC3k dataset from
illumina.

There are many potential directions and improvements to build upon this work,
including a more detailed interrogation of the equilibrium condition, the
effects of various noise estimates on the diffusion map, and inclusion of
anisotropic diffusion. Moreover, since our approach is based on the Fokker-
Planck equation, it provides a framework for studying the diffusion limit of
scRNAseq data. In particular, it offers a way to connect the “macroscopic"
approach of manifold learning to the more rigorous “microscopic" approach of
chemical master equations as in [1, 2]. These two approaches should coincide
in the limit where the number of counts per cell is large, and using Fokker-
Planck diffusion maps can help unite these two disparate approaches to the
field and put the analysis of scRNAseq on solid physical ground. Additionally,
we hope that our approach can offer ways to simulate the dynamics of single
cells using approaches such as those laid out in [13] and related work.

## References

  * Gorin and Pachter [2022] G. Gorin and L. Pachter, Distinguishing biophysical stochasticity from technical noise in single-cell rna sequencing using monod 10.1101/2022.06.11.495771 (2022).
  * Gorin  _et al._ [2023] G. Gorin, J. J. Vastola, and L. Pachter, Studying stochastic systems biology of the cell with single-cell genomics data 10.1101/2023.05.17.541250 (2023).
  * Coifman and Lafon [2006] R. R. Coifman and S. Lafon, Diffusion maps, Applied and Computational Harmonic Analysis 21, 5 (2006).
  * Coifman  _et al._ [2005] R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, F. Warner, and S. W. Zucker, Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps, Proceedings of the National Academy of Sciences 102, 7426 (2005).
  * Ting  _et al._ [2011] D. Ting, L. Huang, and M. Jordan, An analysis of the convergence of graph laplacians 10.48550/ARXIV.1101.5435 (2011).
  * Huang [2011] S. Huang, The molecular and mathematical basis of waddington's epigenetic landscape: A framework for post-darwinian biology?, BioEssays 34, 149 (2011).
  * Wang  _et al._ [2011] J. Wang, K. Zhang, L. Xu, and E. Wang, Quantifying the Waddington landscape and biological paths for development and differentiation, Proceedings of the National Academy of Sciences 108, 8257 (2011), publisher: Proceedings of the National Academy of Sciences.
  * Ferrell [2012] J. E. Ferrell, Bistability, bifurcations, and waddington's epigenetic landscape, Current Biology 22, R458 (2012).
  * Haghverdi  _et al._ [2015] L. Haghverdi, F. Buettner, and F. J. Theis, Diffusion maps for high-dimensional single-cell analysis of differentiation data, Bioinformatics 31, 2989 (2015).
  * Haghverdi  _et al._ [2016] L. Haghverdi, M. Büttner, F. A. Wolf, F. Buettner, and F. J. Theis, Diffusion pseudotime robustly reconstructs lineage branching, Nature Methods 13, 845 (2016).
  * Wolf  _et al._ [2018] F. A. Wolf, P. Angerer, and F. J. Theis, SCANPY: large-scale single-cell gene expression data analysis, Genome Biology 19, 10.1186/s13059-017-1382-0 (2018).
  * Stuart  _et al._ [2019] T. Stuart, A. Butler, P. Hoffman, C. Hafemeister, E. Papalexi, W. M. Mauck, Y. Hao, M. Stoeckius, P. Smibert, and R. Satija, Comprehensive integration of single-cell data, Cell 177, 1888 (2019).
  * Coifman  _et al._ [2008] R. R. Coifman, I. G. Kevrekidis, S. Lafon, M. Maggioni, and B. Nadler, Diffusion maps, reduction coordinates, and low dimensional representation of stochastic systems, Multiscale Modeling and Simulation 7, 842 (2008).
  * van Kampen [1981] N. G. van Kampen, Itô versus stratonovich, Journal of Statistical Physics 24, 175 (1981).
  * Coomer  _et al._ [2022] M. A. Coomer, L. Ham, and M. P. Stumpf, Noise distorts the epigenetic landscape and shapes cell-fate decisions, Cell Systems 13, 83 (2022).
  * Loftsgaarden and Quesenberry [1965] D. O. Loftsgaarden and C. P. Quesenberry, A nonparametric estimate of a multivariate density function, The Annals of Mathematical Statistics 36, 1049–1051 (1965).
  * Nadler  _et al._ [2006] B. Nadler, S. Lafon, R. R. Coifman, and I. G. Kevrekidis, Diffusion maps, spectral clustering and reaction coordinates of dynamical systems, Applied and Computational Harmonic Analysis 21, 113 (2006).