11institutetext: Universitat Pomepu Fabra, Barcelona, Spain
11email: {firstname<EMAIL_ADDRESS>
22institutetext: Università degli Studi di Cagliari, Cagliari, Italy
22email<EMAIL_ADDRESS>
33institutetext: Universitat Pomepu Fabra, Barcelona, Spain
33email<EMAIL_ADDRESS>
# Poisoning Attacks on Algorithmic Fairness

David Solans 11 0000-0001-6979-9330 Battista Biggio 22 0000-0001-7752-509X
Carlos Castillo 33 0000-0003-4544-0416 Anonymized for double-blind review

###### Abstract

Research in adversarial machine learning has shown how the performance of
machine learning models can be seriously compromised by small perturbations in
the training data. While the effects on model accuracy of adversarial attacks
have been widely studied in the field, their potential effects on other model
performance metrics remains to be evaluated. In this work, we introduce a
theoretical analysis of how different types of adversarial machine learning
attacks can be crafted to compromise a supervised classifiers from the
perspective of algorithmic discrimination. We then empirically test the
effects of two different adversarial attacks, crafted with the purpose of
increasing disparate impact and disparate mistreatment, two common measures
used to evaluate algorithmic discrimination. Our results show how
classification disparities can be introduced or exaggerated through small
changes in the input data and under limited previous knowledge from the
attacker.

###### Abstract

Research in adversarial machine learning has shown how the performance of
machine learning models can be seriously compromised by small perturbations in
the training data. While the effects on model accuracy of adversarial attacks
have been widely studied in the field, their potential effects on other model
performance metrics remains to be evaluated. In this work, we introduce a
theoretical analysis of how state of the art on adversarial machine learning
attacks can be adapted to compromise a supervised classifiers from the
perspective of algorithmic discrimination. We then empirically test the
effects of adversarial attacks, crafted with the purpose of increasing
disparate impact and disparate mistreatment, two common measures used to
evaluate algorithmic discrimination. Our results show how classification
disparities can be introduced or exaggerated through small changes in the
input data and under limited previous knowledge from the attacker.

###### Abstract

Research in adversarial machine learning has shown how the performance of
machine learning models can be seriously compromised by small perturbations in
the training data. While the effects on model accuracy of adversarial attacks
have been widely studied in the field, their potential effects on other model
performance metrics remains to be evaluated. In this work, we introduce a
theoretical analysis of how different types of adversarial machine learning
attacks can be crafted to compromise a supervised classifiers from the
perspective of algorithmic discrimination. We then empirically test the
effects of two different adversarial attacks, crafted with the purpose of
increasing disparate impact and disparate mistreatment, two common measures
used to evaluate algorithmic discrimination. Our results show how
classification disparities can be introduced or exaggerated through small
changes in the input data and under limited previous knowledge from the
attacker.

###### Keywords:

algorithmic discrimination algorithmic fairness poisoning attacks adversarial
machine learning machine learning security

###### Keywords:

algorithmic discrimination algorithmic fairness poisoning attacks adversarial
machine learning machine learning security

###### Keywords:

algorithmic discrimination algorithmic fairness adversarial machine learning
attacks

<<<<<<< HEAD

##  1 Introduction

TODO: revise intro and abs. to be more convincing. I would go for this
structure: – Battista

  1. 1.

Machine learning and problem of unfairness. Work that tries to make ML
algorithms fairer. TODO: We need a paragraph on fairness. – Battista

  2. 2.

Problem: when dealing with malicious actors that can tamper with training data
- e.g., if such data is collected in the wild - the classifier may be subject
to adversarial data manipulation. Poisoning attacks (affecting availability
and integrity) etc. These attacks have been widely demonstrated in adversarial
settings .

  3. 3.

In this work, we show that an attacker may also be able to worsen algorithmic
fairness by developing a novel poisoning attack. […]

  4. 4.

We show that (experiments)

  5. 5.

We discuss related work in Sect. LABEL:, while we conclude the paper with a
discussion on its main contributions and interesting future developments.

_Algoritmic Fairness_ is an emerging concern in computing that started within
the data mining community but has extended into other fields including machine
learning, information retrieval, and theory of algorithms
[hajian2016algorithmic]. It deals with the design of algorithms and decision
support systems that are non-discriminatory, i.e., that do not introduce an
unjustified disadvantage for members of a group, and particularly that do not
further place at a disadvantage members of an already disadvantaged social
group. In machine learning, the problem that has been most studied to date is
supervised classification, in which algorithmic fairness methods have been
mostly proposed to fulfill criteria related to parity (equality)
[zafar2017parity]. Most of the methods proposed to date assume belevolence
from the part of the data scientist or developer creating the classification
model: she is envisioned as an actor trying to eliminate or reduce potential
discrimination in her model.

The problem arises when dealing with malicious actors that can tamper with the
model development, for instance by tampering with training data.
Traditionally, _poisoning attacks_ have been studied in _Adversarial Machine
Learning_. These attacks are usually crafted with the purpose of increasing
the misclassification rate in a machine learning model, either for certain
samples or in an indiscriminate basis, and have been widely demonstrated in
adversarial settings (see, e.g., [biggio2018wild]).

In this work, we show that an attacker may be able to introduce algorithmic
discrimination by developing a novel poisoning attack. The purpose of this
attacker is to create or increase a disadvantage against a specific group of
individuals or samples.

NOTE: Text below comes from the previous introduction, it can be incorporated
following the schema proposed by Battista –ChaTo

Motivation. The main goal of this paper is to show the potential harm that an
attacker can cause in a machine learning system if the attacker can manipulate
its training data. For instance, the developer of a criminal recidivism
prediction tool [ProPublica_article] could sample training data in a
discriminatory manner to bias the tool against a certain group of people.
Similar harms can occur when training data is collected from public sources,
such as surveys done online that cannot be fully trusted. A minority of ill-
intentioned users could _poison_ this data with the objective of causing
certain defects in the machine learning system created from it.

In addition to these examples, there is the unintentional setting, where
inequities are introduced in the machine learning model as an undesired effect
of the data collection or data labeling. For instance, human annotators could
systematically make mistakes when assigning labels to images of people of a
certain skin color [buolamwini2018].

The methods we describe on this paper could be used to model the potential
harm to a machine learning system in the worst case scenario, demonstrating
the undesired effects that a limited amount of wrongly labeled samples can
cause, even if created in an unwanted manner.

Contributions. This work first characterizes several types of possible
adversarial attacks that could be crafted with the purpose of introducing
algorithmic discrimination in a model. Different types of attacks are defined
in terms of the attacker knowledge, capabilities, and strategy.

After this characterization, we perform experiments in two scenarios: a
“black-box” attack in which the attacker only has access to the training data,
but not the model, and a “white-box” scenario in which the attacker has access
to the model. The effects of these attacks are measured using impact
quantification metrics. The experiments show that by carefully perturbing a
limited amount of training examples, an skilled attacker has the possibility
of introducing different types of inequities for certain groups of
individuals. This, can be done without large effects on the overall accuracy
of the system, which makes these attacks harder to detect.

The rest of this paper is organized as follows. Section LABEL:sec:background,
provides background on adversarial machine learning and algorithmic fairness.
Section describes how an existing categorization of different types of
adversarial machine learning attacks can be extended to fit with the new type
of attacks introduced in this paper. Section demonstrates empirically the
feasibility of the new types of attacks on both synthetic and real-world data.
Section provides further references to related work. Section presents our
conclusions.

##  2 Poisoning Fairness

Notation. We denote the sample and label spaces with $\mathcal{X}$ and
$\mathcal{Y}\in[0,1]$, respectively. Training data is denoted as
$\mathcal{D}_{tr}=(\mathbf{x}_{i},y_{i})^{n}_{i=1}$, where $n$ is the training
set size. $\mathcal{L}(\mathcal{D},\mathbf{\theta})$ is used to denote the
loss incurred by a classifier
$f_{\mathbf{\theta}}:\mathcal{X}\rightarrow\mathcal{Y}$, parametrized by
$\mathbf{\theta}$, on the data $\mathcal{D}$. Typically, the loss is computed
by averaging the loss on each data point, e.g.,
$\mathcal{L}(\mathcal{D}_{tr},\mathbf{\theta})=\frac{1}{n}\sum_{i=1}^{n}\ell(\mathbf{x}_{i},y_{i},\mathbf{\theta})$.

TODO: any notation on fairness which may be useful here?

Attack Formulation. the optimal attack strategy can be formulated as a bilevel
optimization problem as follows:

$\displaystyle\max_{\mathbf{x}_{c}}$
$\displaystyle\sum_{k=1}^{p}\ell(\mathbf{x}_{k},y_{k},\mathbf{\theta}^{\star})-\lambda\sum_{j=1}^{m}\ell(\mathbf{x}_{j},y_{j},\mathbf{\theta}^{\star})$
(1) $\displaystyle{\rm s.t.}$
$\displaystyle\mathbf{\theta}^{\star}\in\operatornamewithlimits{\arg\,\min}_{\mathbf{\theta}}\mathcal{L}(\mathcal{D}_{\rm
tr}\cup(\mathbf{x}_{c},y_{c}),\mathbf{\mathbf{\theta}})$ (2)

The validation loss (outer max problem) consists of maximizing the error on
the minority group of $p$ points (indexed by $k$), while preserving the loss
on the other $m$ points (indexed by $j$). The trade-off parameter $\lambda$
aims to balance the contribution of the two terms in the validation loss.

We can explain here what we do to worsen disparate impact or mistreatment, and
then report the algorithm. TODO: can we show that the objective above is
related to such measures?

Attack Algorithm This is a projected gradient descent … (report it as an
Algorithm)

###  2.1 Attack Scenarios

White-box Attack.

Black-box Transfer Attack.

##  3 Experiments

This section describes the obtained results for two different types of
attacks. We have selected one example of a Perfect-Knowledge (PK) White-Box
Attack where the original model is fully known for the attacker. Then, we
tested a set up of Limited-Knowledge (PK) Gray-Box Attack where the attacker
has access to the train data, so a surrogate model can be trained and the
attack can be then optimized using the learning function of this model. This
scenarios were selected

======= >>>>>>> 88d64e042af2eccb9ab415d1041a6e17ebeae9c9

##  4 Introduction

The literature on Adversarial Machine Learning has described several types of
adversarial attacks crafted with the purpose of increasing the
misclassification rate in a machine learning model, either for certain samples
or in an indiscriminate basis (see, e.g., [biggio2018wild]). In this work, we
explore how analogous techniques can be used to compromise a machine learning
model, not to drive its accuracy down, but with the purpose of adding
algorithmic discrimination, or exaggerating it if it already exists. In other
words, the purpose of the attacker will be to create or increase a
disadvantage against a specific group of individuals or samples.

###  4.1 Motivation

The main goal of this paper is to show the potential harm that an attacker can
cause in a machine learning system if the attacker can manipulate its training
data. For instance, the developer of a criminal recidivism prediction tool
[ProPublica_article] could sample training data in a discriminatory manner to
bias the tool against a certain group of people. Similar harms can occur when
training data is collected from public sources, such as surveys done online
that cannot be fully trusted. A minority of ill-intentioned users could
_poison_ this data with the objective of causing certain defects in the
machine learning system created from it.

In addition to these examples, there is the unintentional setting, where
inequities are introduced in the machine learning model as an undesired effect
of the data collection or data labeling. For instance, human annotators could
systematically make mistakes when assigning labels to images of people of a
certain skin color [buolamwini2018].

The methods we describe on this paper could be used to model the potential
harm to a machine learning system in the worst case scenario, demonstrating
the undesired effects that a limited amount of wrongly labeled samples can
cause, even if created in an unwanted manner.

###  4.2 Contributions

This work first characterizes several types of possible adversarial attacks
that could be crafted with the purpose of introducing algorithmic
discrimination in a model. Different types of attacks are defined in terms of
the attacker knowledge, capabilities, and strategy.

After this characterization, we perform experiments in two scenarios: a
“black-box” attack in which the attacker only has access to the training data,
but not the model, and a “white-box” scenario in which the attacker has access
to the model. The effects of these attacks are measured using impact
quantification metrics. The experiments show that by carefully perturbing a
limited amount of training examples, an skilled attacker has the possibility
of introducing different types of inequities for certain groups of
individuals. This, can be done without large effects on the overall accuracy
of the system, which makes these attacks harder to detect.

The rest of this paper is organized as follows. Section LABEL:sec:background,
provides background on adversarial machine learning and algorithmic fairness.
Section describes how an existing categorization of different types of
adversarial machine learning attacks can be extended to fit with the new type
of attacks introduced in this paper. Section demonstrates empirically the
feasibility of the new types of attacks on both synthetic and real-world data.
Section provides further references to related work. Section presents our
conclusions and future work.


# Angle-Resolved Pair Photoemission Theory for Correlated Electrons

Thomas P. Devereaux  Stanford Institute for Materials and Energy Sciences,
SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park,
California 94025, USA  Geballe Laboratory for Advanced Materials, Stanford
University, Stanford, California 94305, USA  Department of Materials Science
and Engineering, Stanford University, Stanford, California 94305, USA  Martin
Claassen  Department of Physics and Astronomy, University of Pennsylvania,
Philadelphia, PA 19104, USA  Xu-Xin Huang  Stanford Institute for Materials
and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill
Road, Menlo Park, California 94025, USA  Michael Zaletel  Department of
Physics, University of California, Berkeley, California 94720, USA  Materials
Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California
94720, USA  Joel E. Moore  Department of Physics, University of California,
Berkeley, California 94720, USA  Materials Sciences Division, Lawrence
Berkeley National Laboratory, Berkeley, California 94720, USA  Dirk Morr
University of Illinois at Chicago, Chicago, Illinois 60607, USA  Fahad Mahmood
University of Illinois at Urbana Champaign, Champaign, Illinois 61801, USA
Peter Abbamonte  University of Illinois at Urbana Champaign, Champaign,
Illinois 61801, USA  Zhi-Xun Shen  Stanford Institute for Materials and Energy
Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo
Park, California 94025, USA  Geballe Laboratory for Advanced Materials,
Stanford University, Stanford, California 94305, USA  Departments of Physics
and Applied Physics, Stanford University, Stanford, CA 94305, USA

###### Abstract

In this paper we consider the possibility and conditions for pair
photoemission whereby two incident photons emit pairs of electrons from a
candidate material as a novel method to measure and visualize electronic
correlations. As opposed to “double photoemission” - where a single photon
precipitates the ejection of a pair electrons via a subsequent electron energy
loss scattering process - we show that pair photoemission need not be limited
to interference between initial photoelectrons and valence electrons, and
moreover, can occur without the energy penalty of two work functions. This
enables detection of pairs of electrons at high energy resolution that may be
correlated in the same quantum many-body states.

##  I Introduction

Over the past decades, angle-resolved photo-emission spectroscopy (ARPES) has
emerged as a paradigmatic experimental probe of electronic structure and
correlations, band topology or surface states, unconventional
superconductivity or the enigmatic pseudogap phase, granting insight to
characterize electronic behavior in new quantum materials. By measuring the
kinetic energy and angular dependence of photo-emitted electrons, ARPES
supplies information on the energy and momentum dependence of valence
electrons in a material, and is widely understood to reflect to a good
approximation the behavior of the single-particle spectral function [1].

Higher order photoemission processes have been utilized to further obtain
information beyond the single-particle density of states. In “double
photoemission” for example, a highly energetic photon causes the emission of
an electron which may cause a second electron to be photoemitted via the
Coulomb interaction if it can impart enough energy for the second electron to
escape to a detector [2]. For example, a photoemitted core electron may be
accompanied by Auger electron emission, whereby the energy emitted by Auger
decay of the core hole is utilized to cause another electron to be emitted
[3]. Such “shake-off” or “secondaries” spectra contained both photoemitted
core and Auger electrons. The fact that the energies of the two electrons can
be themselves continuous yet sum to conserve energy can show that the
electrons are correlated, and a comparison with single-particle photoemission
can be utilized to determine the so-called “exchange-correlation” hole energy
[4].

In analogy to photon- or electron-based coincidence spectroscopies, recently
an interesting proposal suggested extending ARPES to use energy and angle-
resolved coincidence detection to account for two-photon two-electron photo-
emission events and extract two-particle Bethe-Salpeter wave functions [6] of
valence electrons of the material. Here, in contrast to double photoemission
due to Coulomb drag, the coincidence signal derives from two-photon absorption
at lower photon energies. Measurement of the angle dependence of two-photon
coincidence events at the detector hence importantly permits resolving the
momenta and energies of the ejected electron pairs, without the need to
“disentangle” highly-complex “Coulomb drag” processes that complicates the
study of important low energy effects.

While the possibility to extract electronic correlations from coincidence
counts in ARPES immediately suggests a variety of applications such as
elucidating unconventional pairing mechanisms in high-temperature
superconductors or heavy-fermion compounds, a key question concerns
understanding exactly the nature of what this new probe actually measures in a
correlated electron system, and how to interpret its result in terms of more
intuitive quantities such as pair correlation functions or the superconducting
gap. Indeed, it is straightforward to see that the coincidence signal does not
map onto more readily interpretable superconducting pair correlation
functions, since pairs of detected electrons that comprise a coincidence
signal are not necessarily ejected from the sample at the same time. On the
other hand, a more microscopic description of the pair ARPES cross section is
necessary, to permit a formal accounting for final-state effects and the
detector geometry and differentiate from double photoemission.

In this work, we address these questions by developing a generic theoretical
description of angle-resolved pair photoemission and studying its behavior in
light of superconducting instabilities in the attractive and repulsive Hubbard
model on small clusters. We show that pair photoemission need not be limited
to interference between initial photoelectrons and valence electrons, and
moreover, can occur without the energy penalty of two work functions. This
enables detection of pairs of electrons at high energy resolution that may be
correlated in the same quantum many-body states.

Figure 1: Setup:  Two electrons at $(\bar{\mathbf{r}},\tau)$,
$(\bar{\mathbf{r}}^{\prime},\tau^{\prime})$ are photo-emitted upon absorbing
two photons with frequency $\omega_{\textrm{ph}}$. Coincidence detection takes
place at time $t$ and positions $\mathbf{x}$, $\mathbf{x}^{\prime}$.

Fig. 1 depicts a schematic of the pair photoemission process for a two-
dimensional sample or surface state. Two photons with energy
$\hbar\omega_{\textrm{ph}}$ eject two electrons from the sample, which are
subsequently observed at the detector at the same time $t$ with both angle and
energy resolution. For simplicity, but without loss of generality we ignore
bulk effects, and henceforth denote three-dimensional positions and momenta
using bold notation $\mathbf{r}$, whereas their two-dimensional components in
the sample plane are denoted by $\bar{\mathbf{r}}$. Suppose that sample and
emitted electrons are described by fields $\hat{\Phi}(\bar{\mathbf{r}})$ and
$\hat{\Psi}(\mathbf{r})$, respectively (we suppress implicit spin indices, for
conciseness), and are governed by a generic Hamiltonian $\hat{H}$

$\displaystyle\hat{H}_{0}$
$\displaystyle=\hat{H}_{\textrm{valence}}(\hat{\Phi})+\hat{H}_{\textrm{emitted}}(\hat{\Psi})+\hat{H}_{\textrm{v-e}}(\hat{\Phi},\hat{\Psi})$
(1)

such that emitted electrons $\hat{\Psi}$ behave as freely-propagating waves at
long distances from the sample while appropriately encapsulating final-state
effects (inverse LEED) as well as possible back actions
$\hat{H}_{\textrm{v-e}}(\hat{\Phi},\hat{\Psi})$ which would be important for
Coulomb-drag mediated double photoemission.

The photoemission process $\hat{H}_{\textrm{el-ph}}$ now takes a sample
electron $\hat{\Phi}(\bar{\mathbf{r}})$ to a propagating final state
$\hat{\Psi}(\mathbf{r})$

$\displaystyle\hat{H}_{\textrm{el-ph}}(t)$ $\displaystyle=s(t)\int
d^{3}\mathbf{r}~{}g(\mathbf{r})~{}\hat{\Psi}^{\dagger}(\mathbf{r})\hat{\Phi}_{S}(\bar{\mathbf{r}})~{}e^{-i\omega_{\textrm{ph}}t}+\textrm{h.c.}$
(2)

where $g(\mathbf{r})$ is the dipole matrix element and $s(t)$ describes a
Gaussian probe pulse envelope with

$\displaystyle s(t)=e^{-(t-t_{0})^{2}/2\sigma_{\textrm{pr}}^{2}}$ (3)

Subsequently, the photo-electron detector measures the mean momentum
$\mathds{k}$ of propagating electron wave packets, described by a photo-
current

$\displaystyle\langle\hat{J}_{\mathds{k}}\rangle=\frac{\mathds{k}}{e}\iint
d\mathbf{x}d\mathbf{x}^{\prime}~{}\phi_{\mathds{k}}^{\star}(\mathbf{x})\phi_{\mathds{k}}(\mathbf{x}^{\prime})~{}\langle\hat{\Psi}^{\dagger}(\mathbf{x})\hat{\Psi}(\mathbf{x}^{\prime})\rangle$
(4)

where $\phi_{\mathds{k}}(\mathbf{r})$ denotes a wave packet centered at the
detector location.

##  II Formalism

###  II.1 Single-Electron ARPES

The conventional “single-electron” ARPES signal now follows straightforwardly
[5] from a perturbative expansion in $\hat{H}_{\textrm{el-ph}}$ of the
measured photocurrent

$\displaystyle I_{\mathds{k}}$ $\displaystyle=\int_{-\infty}^{t}d\tau
d\tau^{\prime}s(\tau)s(\tau^{\prime})e^{i\omega_{\textrm{ph}}(\tau-\tau^{\prime})}$
$\displaystyle\times\int
d\mathbf{x}d\mathbf{x}^{\prime}\phi_{\mathds{k}}^{\star}(\mathbf{x})\phi_{\mathds{k}}(\mathbf{x}^{\prime})\int
d\bar{\mathbf{r}}d\bar{\mathbf{r}}^{\prime}g^{\star}(\mathbf{r})g(\mathbf{r}^{\prime})$
$\displaystyle\times\langle\hat{\Phi}^{\dagger}(\bar{\mathbf{r}},\tau)\hat{\Psi}(\mathbf{r},\tau)\hat{\Psi}^{\dagger}(\mathbf{x},t)\hat{\Psi}(\mathbf{x}^{\prime},t)\hat{\Psi}^{\dagger}(\mathbf{r}^{\prime},\tau^{\prime})\hat{\Phi}(\bar{\mathbf{r}}^{\prime},\tau^{\prime})\rangle$
(5)

where $\langle\cdot\rangle=\textrm{tr}\\{\cdot~{}e^{-\beta\hat{H}_{0}}\\}/Z$
denotes thermal expectation values with respect to $\hat{H}_{0}$. If back
action $\hat{H}_{\textrm{valence-emitted}}$ between emitted and valence
electrons can be neglected, this expression simplifies drastically, as
$\langle\hat{\Phi}^{\dagger}(\bar{\mathbf{r}},\tau)\hat{\Psi}(\mathbf{r},\tau)\hat{\Psi}^{\dagger}(\mathbf{x},t)\hat{\Psi}(\mathbf{x}^{\prime},t)\hat{\Psi}^{\dagger}(\mathbf{r}^{\prime},\tau^{\prime})\hat{\Phi}(\bar{\mathbf{r}}^{\prime},\tau^{\prime})\rangle=\langle\hat{\Phi}^{\dagger}(\bar{\mathbf{r}},\tau)\hat{\Phi}(\bar{\mathbf{r}}^{\prime},\tau^{\prime})\rangle\langle\hat{\Psi}(\mathbf{r},\tau)\hat{\Psi}^{\dagger}(\mathbf{x},t)\rangle\langle\hat{\Psi}(\mathbf{x}^{\prime},t)\hat{\Psi}^{\dagger}(\mathbf{r}^{\prime},\tau^{\prime})\rangle$.
Furthermore, assuming a single electronic valence band
$\hat{\Phi}(\bar{\mathbf{r}})=\sum_{\bar{\mathbf{k}}}u_{\bar{\mathbf{k}}}(\bar{\mathbf{r}})e^{i\bar{\mathbf{k}}\bar{\mathbf{r}}}\hat{c}_{\bar{\mathbf{k}}}^{\vphantom{{\dagger}}}$
with Bloch function $u_{\bar{\mathbf{k}}}(\bar{\mathbf{r}})$, and neglecting
the detector wave packet shape functions $\phi_{\mathds{k}}(\mathbf{x})\to
e^{i\mathds{k}\mathbf{x}}$ (thereby discarding time-of-flight information),
one arrives at ($\mathds{k}\to\mathbf{k}$):

$\displaystyle I_{\mathbf{k}}=-i\int_{-\infty}^{t}d\tau d\tau^{\prime}$
$\displaystyle
e^{i\omega_{\textrm{ph}}(\tau-\tau^{\prime})}s(\tau)s(\tau^{\prime})\left|M_{\mathbf{k}}\right|^{2}\times$
$\displaystyle\times~{}G^{<}_{\bar{\mathbf{k}}}(\tau,\tau^{\prime})\mathcal{G}^{\star}_{\mathbf{k}}(\tau,t)\mathcal{G}_{\mathbf{k}}(t,\tau^{\prime})$
(6)

where $M_{\mathbf{k}}$ is a matrix element evaluated from $g(\mathbf{r})$ and
the Bloch function of the single valence band,
$G^{<}_{\bar{\mathbf{k}}}(\tau,\tau^{\prime})$ is the lesser sample Green’s
function

$\displaystyle
G^{<}_{\bar{\mathbf{k}}}(\tau,\tau^{\prime})=i\langle\hat{c}_{\bar{\mathbf{k}}}^{{\dagger}}(\tau)\hat{c}_{\bar{\mathbf{k}}}^{\vphantom{{\dagger}}}(\tau^{\prime})\rangle$
(7)

and

$\displaystyle\mathcal{G}_{\mathbf{k}}(t,t^{\prime})=-i\langle\hat{\mathcal{T}}\hat{c}_{\mathbf{k}}^{\vphantom{{\dagger}}}(t)\hat{c}_{\mathbf{k}}^{{\dagger}}(t^{\prime})\rangle$
(8)

is the propagating electron Green’s function for an inverse LEED state.
Finally, a drastically simplified expression can be provided, if
$\mathcal{G}_{\mathbf{k}}$ is approximated by a free electron Green’s function
with dispersion $\epsilon_{\mathbf{k}}=\mathbf{k}^{2}/2m_{0}$. Then, defining
the energy $\omega$ observed at the detector as

$\displaystyle\omega\equiv\omega_{\textrm{ph}}-\frac{\mathbf{k}^{2}}{2m_{0}}-W$
(9)

where $W$ is the work function of the sample, and taking $t\to\infty$, one
finally arrives at

$\displaystyle I_{\mathbf{k}}=i\int_{-\infty}^{\infty}d\tau
d\tau^{\prime}e^{i\omega(\tau-\tau^{\prime})}s(\tau)s(\tau^{\prime})\left|M_{\mathbf{k}}\right|^{2}G^{<}_{\bar{\mathbf{k}}}(\tau,\tau^{\prime})$
(10)

This is the usual expression for single-particle ARPES in terms of
convolutions of shape functions, matrix elements and the lesser Green’s
function[5].

###  II.2 Angle-resolved Pair Photoemission

Similarly, a coincidence measurement signal can be defined as

$\displaystyle\langle\hat{J}_{\mathds{k}_{1}\mathds{k}_{2}}\rangle=\frac{\mathds{k}_{1}\mathds{k}_{2}}{e^{2}}\int
d\mathbf{x}_{1}d\mathbf{x}_{1}^{\prime}d\mathbf{x}_{2}d\mathbf{x}_{2}^{\prime}~{}\phi_{\mathds{k}_{1}}^{\star}(\mathbf{x}_{1})\phi_{\mathds{k}_{2}}^{\star}(\mathbf{x}^{\prime}_{1})\phi_{\mathds{k}_{2}}(\mathbf{x}^{\prime}_{2})$
$\displaystyle~{}~{}~{}~{}\times\sum_{\nu\nu^{\prime}}\phi_{\mathds{k}_{1}}(\mathbf{x}_{2})~{}\langle\hat{\Psi}_{\nu}^{\dagger}(\mathbf{x}_{1},t)\hat{\Psi}_{\nu^{\prime}}^{\dagger}(\mathbf{x}^{\prime}_{1},t)\hat{\Psi}_{\nu^{\prime}}(\mathbf{x}^{\prime}_{2},t)\hat{\Psi}_{\nu}(\mathbf{x}_{2},t)\rangle$
(11)

In complete analogy to single-electron ARPES, the photo-detection rate can now
be evaluated from a perturbative expansion in $\hat{H}_{\textrm{el-ph}}$. To
first order the response involves a single photoemission vertex and vanishes.
We note that this contribution is essential for Coulomb-mediated double
photoemission for high photon energies. Here, a perturbative expansion in
$\hat{H}_{\textrm{v-e}}(\hat{\Phi},\hat{\Psi})$ additionally accounts for the
Coulomb interaction mediated back action of the photo emitted electron,
imparting enough energy on a second sample electron to eject it, rendering the
coincidence signal non-zero. As discussed above, we are primarily interested
in two-photon two-electron pair ARPES processes at lower photon energy; in
this regime, the double emission contribution is negligible for energetic
reasons.

To second order in $\hat{H}_{\textrm{el-ph}}$, the two-photon two-electron
coincidence photo-detection signal formally reads

$\displaystyle D_{\mathds{k}_{1}\mathds{k}_{2}}$
$\displaystyle=\sum_{\begin{subarray}{c}\sigma_{1}\sigma^{\prime}_{1}\nu\\\
\sigma_{2}\sigma^{\prime}_{2}\nu^{\prime}\end{subarray}}~{}\int\displaylimits_{-\infty}^{t}d\tau_{1}d\tau_{2}\int\displaylimits_{-\infty}^{\tau_{1}}d\tau^{\prime}_{1}\int\displaylimits_{-\infty}^{\tau_{2}}d\tau^{\prime}_{2}$
$\displaystyle\times
e^{i\omega_{\textrm{ph}}(\tau_{1}+\tau^{\prime}_{1}-\tau_{2}-\tau^{\prime}_{2})}s(\tau_{1})s(\tau^{\prime}_{1})s(\tau_{2})s(\tau^{\prime}_{2})$
$\displaystyle\times\int
d\mathbf{r}_{1}d\mathbf{r}^{\prime}_{1}d\mathbf{r}_{2}d\mathbf{r}^{\prime}_{2}g^{\star}(\mathbf{r}_{1})g^{\star}(\mathbf{r}^{\prime}_{1})g(\mathbf{r}^{\prime}_{2})g(\mathbf{r}_{2})$
$\displaystyle\times\int
d\mathbf{x}_{1}d\mathbf{x}^{\prime}_{1}d\mathbf{x}_{2}d\mathbf{x}^{\prime}_{2}~{}\phi_{\mathds{k}_{1}}^{\star}(\mathbf{x}_{1})\phi_{\mathds{k}_{2}}^{\star}(\mathbf{x}^{\prime}_{1})\phi_{\mathds{k}_{2}}(\mathbf{x}^{\prime}_{2})\phi_{\mathds{k}_{1}}(\mathbf{x}_{2})$
$\displaystyle\times\langle\left<\hat{\Phi}_{\sigma_{1}}^{\dagger}(\mathbf{r}_{1},\tau_{1})\hat{\Psi}_{\sigma_{1}}(\mathbf{r}_{1},\tau_{1})\hat{\Phi}_{\sigma^{\prime}_{1}}^{\dagger}(\mathbf{r}^{\prime}_{1},\tau^{\prime}_{1})\hat{\Psi}_{\sigma^{\prime}_{1}}(\mathbf{r}^{\prime}_{1},\tau^{\prime}_{1})\right.$
$\displaystyle~{}~{}~{}\times\hat{\Psi}_{\nu}^{\dagger}(\mathbf{x}_{1},t)\hat{\Psi}_{\nu^{\prime}}^{\dagger}(\mathbf{x}^{\prime}_{1},t)\hat{\Psi}_{\nu^{\prime}}(\mathbf{x}^{\prime}_{2},t)\hat{\Psi}_{\nu}(\mathbf{x}_{2},t)$
$\displaystyle~{}~{}~{}\times\left.\hat{\Psi}_{\sigma^{\prime}_{2}}^{\dagger}(\mathbf{r}^{\prime}_{2},\tau^{\prime}_{2})\hat{\Phi}_{\sigma^{\prime}_{2}}(\mathbf{r}^{\prime}_{2},\tau^{\prime}_{2})\hat{\Psi}_{\sigma_{2}}^{\dagger}(\mathbf{r}_{2},\tau_{2})\hat{\Phi}_{\sigma_{2}}(\mathbf{r}_{2},\tau_{2})\right>$
(12)

Assuming negligible back action or Coulomb interactions between photo-emitted
electrons and low-energy sample electrons, this daunting multi-point
correlation function can be decomposed in analogy to conventional ARPES. The
coincidence detection rate can be written as

$\displaystyle
D_{\mathds{k}_{1}\mathds{k}_{2}}=\int\displaylimits_{-\infty}^{t}d\tau_{1}d\tau_{2}\int\displaylimits_{-\infty}^{\tau_{1}}d\tau^{\prime}_{1}\int\displaylimits_{-\infty}^{\tau_{2}}d\tau^{\prime}_{2}~{}s(\tau_{1})s(\tau^{\prime}_{1})s(\tau_{2})s(\tau^{\prime}_{2})$
$\displaystyle\times\int
d\mathbf{k}d\mathbf{k}^{\prime}d\mathbf{q}\sum_{\begin{subarray}{c}\sigma_{1}\sigma^{\prime}_{1}\nu\\\
\sigma_{2}\sigma^{\prime}_{2}\nu^{\prime}\end{subarray}}G_{\bar{\mathbf{k}}\bar{\mathbf{k}}^{\prime}\bar{\mathbf{q}}}^{\begin{subarray}{c}\sigma_{1}\sigma_{2}\\\
\sigma^{\prime}_{1}\sigma^{\prime}_{2}\end{subarray}}(\tau_{1},\tau^{\prime}_{1},\tau^{\prime}_{2},\tau_{2})e^{i\omega_{\textrm{ph}}(\tau_{1}+\tau^{\prime}_{1}-\tau_{2}-\tau^{\prime}_{2})}$
$\displaystyle\times\left[F_{\mathbf{k}\mathbf{q}}^{\sigma_{1}\sigma^{\prime}_{1}\nu\nu^{\prime}}(\tau_{1},\tau^{\prime}_{1})\right]^{*}F_{\mathbf{k}^{\prime}\mathbf{q}}^{\sigma_{2}\sigma^{\prime}_{2}\nu\nu^{\prime}}(\tau_{2},\tau^{\prime}_{2})$
(13)

where

$G_{\bar{\mathbf{k}}\bar{\mathbf{k}}^{\prime}\bar{\mathbf{q}}}^{\begin{subarray}{c}\sigma\nu\\\
\sigma^{\prime}\nu^{\prime}\end{subarray}}(\tau_{1},\tau^{\prime}_{1},\tau^{\prime}_{2},\tau_{2})=\langle\hat{c}_{\bar{\mathbf{k}}\sigma}^{{\dagger}}(\tau_{1})\hat{c}_{\bar{\mathbf{q}}-\bar{\mathbf{k}}\sigma^{\prime}}^{{\dagger}}(\tau^{\prime}_{1})\hat{c}_{\bar{\mathbf{q}}-\bar{\mathbf{k}}^{\prime}\nu^{\prime}}^{\vphantom{{\dagger}}}(\tau^{\prime}_{2})\hat{c}_{\bar{\mathbf{k}}^{\prime}\nu}^{\vphantom{{\dagger}}}(\tau_{2})\rangle$
(14)

is the two-particle Green’s function for the sample, and

$\displaystyle
F_{\mathbf{k}\mathbf{q}}^{\sigma\sigma^{\prime}\nu\nu^{\prime}}(\tau,\tau^{\prime})$
$\displaystyle=\int
d\mathbf{p}d\mathbf{p}^{\prime}~{}\phi_{\mathds{k}_{1}}(\mathbf{p})\phi_{\mathds{k}_{2}}(\mathbf{q}-\mathbf{p})~{}M_{\mathbf{q}-\mathbf{k}}M_{\mathbf{k}}$
$\displaystyle\times\left<0\right|\hat{\Psi}_{\mathbf{p}\nu^{\prime}}(t)\hat{\Psi}_{\mathbf{q}-\mathbf{p}\nu}(t)\hat{\Psi}^{\dagger}_{\mathbf{q}-\mathbf{k}\sigma^{\prime}}(\tau^{\prime})\hat{\Psi}^{\dagger}_{\mathbf{k}\sigma}(\tau)\left|0\right>$
(15)

is a four-point function for the photo-emitted electrons which includes the
Fourier-transformed detector shape functions $\phi_{\mathds{k}}(\cdot)$ and
evaluated with respect to the vacuum state. Furthermore, $M_{\mathbf{k}}$
denote photo-excitation matrix elements defined in terms of the dipole matrix
element and valence electron Bloch functions introduced above.
$F_{\mathbf{k}\mathbf{q}}(\tau,\tau^{\prime})$ encodes both propagation and
time-of-flight information, as well as interactions between the two photo-
emitted electrons. A drastic simplification follows from treating emitted
electrons as free fermions. In this case,
$F_{\mathbf{k}\mathbf{q}}(\tau,\tau^{\prime})$ factorizes to

$\displaystyle
F_{\mathbf{k}\mathbf{q}}(\tau,\tau^{\prime})=M_{\mathbf{q}-\mathbf{k}}M_{\mathbf{k}}~{}e^{-i[\epsilon_{\mathbf{k}}(t-\tau)+\epsilon_{\mathbf{q}-\mathbf{k}}(t-\tau^{\prime})]}$
$\displaystyle\times\left[\phi_{\mathds{k}_{1}}(\mathbf{k})\phi_{\mathds{k}_{2}}(\mathbf{q}-\mathbf{k})\delta_{\sigma,\nu^{\prime}}\delta_{\sigma^{\prime},\nu}-\phi_{\mathds{k}_{2}}(\mathbf{k})\phi_{\mathds{k}_{1}}(\mathbf{q}-\mathbf{k})\delta_{\sigma,\nu}\delta_{\sigma,\nu}\right]$
(16)

where $\epsilon_{\mathbf{k}}=\mathbf{k}^{2}/2m_{0}$ is the dispersion of the
photo-emitted electrons. In analogy to the theory for conventional ARPES [5],
we can now make the assumption that the detector wave packet momentum width
can be neglected
$\phi_{\mathds{k}}(\mathbf{k})\to\delta(\mathds{k}-\mathbf{k})$, discarding
again time-of-flight information and dependence on the detector position.
Denote the energies observed at the two detectors minus the photon energy as

$\displaystyle\omega_{1,2}\equiv\omega_{\textrm{ph}}-\frac{\mathds{k}_{1,2}^{2}}{2m_{0}}-W$
(17)

with $W$ the work function of the sample, and taking $t\to\infty$, the
coincidence detection rate can be written as

$\displaystyle D^{(0)}_{\mathbf{k}_{1}\mathbf{k}_{2}}$
$\displaystyle=\int\displaylimits_{-\infty}^{\infty}d\tau_{1}d\tau^{\prime}_{1}d\tau_{2}d\tau^{\prime}_{2}~{}s(\tau_{1})s(\tau^{\prime}_{1})s(\tau_{2})s(\tau^{\prime}_{2})$
$\displaystyle\times\sum_{\sigma\sigma^{\prime}}\langle\hat{\mathcal{T}}\hat{c}_{\bar{\mathbf{k}}_{1}\sigma}^{{\dagger}}(\tau_{1})\hat{c}_{\bar{\mathbf{k}}_{2}\sigma^{\prime}}^{{\dagger}}(\tau^{\prime}_{1})\hat{\mathcal{T}}\hat{c}_{\bar{\mathbf{k}}_{2}\sigma^{\prime}}^{\vphantom{{\dagger}}}(\tau^{\prime}_{2})\hat{c}_{\bar{\mathbf{k}}_{1}\sigma}^{\vphantom{{\dagger}}}(\tau_{2})\rangle$
$\displaystyle\times
e^{i\left[\omega_{1}(\tau_{1}-\tau_{2})+\omega_{2}(\tau^{\prime}_{1}-\tau^{\prime}_{2})\right]}$
(18)

where $\hat{\mathcal{T}}$ denotes time ordering, and we additionally omitted
the photo-excitation matrix elements $M_{\mathbf{k}}$ for conciseness.

###  II.3 Fermi’s Golden Rule

This expression can be recast in a more familiar Fermi’s Golden rule by
inserting complete sets of the states for the $N,N-1,$ and $N-2$ particle
sectors. Also if we neglect the time dependence of the shape functions so that
we only concentrate on frequency resolution, the time integrals can be
performed and the following expression is obtained:

$\displaystyle D^{(0)}_{\mathbf{k}_{1},\mathbf{k}_{2}}(\omega_{1},\omega_{2})$
$\displaystyle=\sum_{n}\mid
M_{0,n}(\mathbf{k}_{1},\mathbf{k}_{2},\omega_{1},\omega_{2})\mid^{2}$
$\displaystyle\times\delta(E_{n}-E_{0}+\omega_{1}+\omega_{2})$ (19)

with $E_{0}$ denoting the $N$ particle ground state energy and $E_{n}$ the
eigenenergies of the $N-2$ particle sector - viz., the expression is simply a
matrix element squared times a term that enforces energy conservation. The
matrix element reads

$\displaystyle M_{0,n}(\mathbf{k}_{1},\mathbf{k}_{2},\omega_{1},\omega_{2})=$
$\displaystyle\sum_{m,\sigma_{1},\sigma_{2}}\left\\{\frac{\left<n\right|{\hat{c}_{\mathbf{k}_{2}\sigma_{2}}^{\vphantom{{\dagger}}}}\left|m\right>\left<m\right|{\hat{c}_{\mathbf{k}_{1}\sigma_{1}}^{\vphantom{{\dagger}}}}\left|0\right>}{E_{m}-E_{0}+\omega_{1}-i\eta}\right.$
$\displaystyle-$
$\displaystyle\left.\frac{\left<n\right|\hat{c}_{\mathbf{k}_{1}\sigma_{1}}^{\vphantom{{\dagger}}}\left|m\right>\left<m\right|\hat{c}_{\mathbf{k}_{2}\sigma_{2}}^{\vphantom{{\dagger}}}\left|0\right>}{E_{m}-E_{0}+\omega_{2}-i\eta}\right\\}$
(20)

with $E_{m}$ the $N-1$ particle sector eigenvalues.

Note that this expression bears a strong resemblance to the Kramers-Heisenberg
expression for resonant inelastic x-ray scattering (RIXS) in which the
manifold of $N-1$ states $\\{\mid m\rangle\langle m\mid\\}$ play the role of
intermediate $N+1$ core hole states whereby a core electron is photoexcited
into the valence band[7]. While for RIXS the final states have the same number
of electrons $N$ as the initial state as the core hole is refilled via photo-
deexcitation, the pair photoemission final states have two less electrons
$N-2$. Despite their apparent differences, the functional form of Eq. (20)
indicates that we would expect resonant pair photoemission whenever one or
both of the frequencies $\omega_{1,2}$ correspond to the $N-1$ removal state
energies observed in photoemission rather than the core-valence transition
energies as in RIXS.

To illustrate the differences between pair photoemission and uncorrelated
single particle photoemission and how information can be obtained from both,
we start by reminding that the ”pairing energy” $\Delta({\bf k_{1},k_{2}})$
for two momentum states ${\mathbf{k}_{1,2}}$ is $\Delta({\bf
k_{1},k_{2}})=E_{2}({\bf k_{1},k_{2}})-E_{1}({\bf k_{1}})-E_{1}({\bf
k_{2}})+E_{0}$ where $E_{N}$ denotes the energies of the $N$ particle removal
states, viz., where single particle photoemission yielding $E_{1}-E_{0}$ and
pair photoemission $E_{2}-E_{0}$. By inspection of Eqs. (19) and (20), we can
see that $E_{2}-E_{0}$ is determined by the overall energy conservation by
$\omega_{1}+\omega_{2}$. In other words, this is given by the slope of the
line connecting $\omega_{1}$ and $\omega_{2}$ for pair photoemission when
plotted as a function of both frequencies. The resonance denominator of Eq.
(20) shows that the intensity on this line is modulated when $\omega_{1,2}$
coincide with the single particle energies $E_{1}-E_{0}$ observed in single
particle photoemission.

###  II.4 Retarded pairing correlator

Suppose that the measured pair emission signal is obtained as a function of
sum and difference frequencies

$\displaystyle\omega=\omega_{1}+\omega_{2}~{},~{}~{}~{}\Delta\omega=\omega_{1}-\omega_{2}$
(21)

By inspection of Eq. (18) one can see that the difference frequency
$\Delta\omega$ parameterizes the “retardation” of the pair emission process
from the sample, i.e. the time delay between emission of the first and second
electron of an observed pair. Integration over the difference frequency
$\Delta\omega$ then yields

$\displaystyle D^{(0)}_{\mathbf{k}_{1}\mathbf{k}_{2}}(\omega)$
$\displaystyle=\int\displaylimits_{-\infty}^{\infty}d\Delta\omega~{}D^{(0)}_{\mathbf{k}_{1}\mathbf{k}_{2}}\left(\frac{\omega+\Delta\omega}{2},\frac{\omega-\Delta\omega}{2}\right)$
$\displaystyle=\int\displaylimits_{-\infty}^{\infty}dt$
$\displaystyle~{}e^{-i\omega
t}\int\displaylimits_{-\infty}^{\infty}d\tau\langle\hat{c}_{\mathbf{k}_{1}}^{{\dagger}}(t)\hat{c}_{\mathbf{k}_{2}}^{{\dagger}}(t+\tau)\hat{c}_{\mathbf{k}_{1}}^{\vphantom{{\dagger}}}(\tau)\hat{c}_{\mathbf{k}_{1}}^{\vphantom{{\dagger}}}(0)\rangle$
(22)

which can be straightforwardly expressed as a spectral decomposition

$\displaystyle D^{(0)}_{\mathbf{k}_{1}\mathbf{k}_{2}}(\omega)$
$\displaystyle=\sum_{n,m}$
$\displaystyle\delta(\omega+E_{0}-E_{n})\left|\left<n\right|\hat{c}_{\mathbf{k}_{1}}^{\vphantom{{\dagger}}}\left|m\right>\left<m\right|\hat{c}_{\mathbf{k}_{2}}^{\vphantom{{\dagger}}}\left|0\right>\right|^{2}$
(23)

Thus coincidence pair ARPES can yield the dynamic superconducting pairing
susceptibility. In a BCS superconductor, the resulting response has a peak at
finite $\omega$ that corresponds to the momentum-dependent superconducting
gap. Spin-resolved pairing as well as pairing that can occur at finite momenta
corresponding to a pair density wave was recently examined in Ref. 8.
Importantly, coincidence pair ARPES can also provide measurements for the
dynamic pair susceptibility in materials at temperatures above the ordered
phase or for systems that may be highly frustrated or condense into a
different, non-superconducting pair state. While dynamic pairing correlations
have been measured via different numerical methods, such as determinant
Quantum Monte Carlo for example [9, 10], susceptibility measurements have been
lacking.

##  III Applications

###  III.1 Free electrons

If the valence electrons within the sample are free, the four-point function
factorizes to

$\displaystyle
G_{\mathbf{k}_{1},\mathbf{k}_{2}}^{(2)}(\tau_{1},\tau^{\prime}_{1},\tau^{\prime}_{2},\tau_{2})\to
G^{<}_{\mathbf{k}_{1}}(\tau_{1}-\tau_{2})~{}G^{<}_{\mathbf{k}_{2}}(\tau^{\prime}_{1}-\tau^{\prime}_{2})$
(24)

for $\mathbf{k}_{1}\neq\mathbf{k}_{2}$ and the coincidence detection rate
becomes a product of single-particle ARPES detection rates

$\displaystyle
D_{\mathbf{k}_{1}\mathbf{k}_{2}}^{(0)}(\omega_{1},\omega_{2})\to
P_{\mathbf{k}_{1}}(\omega_{1})P_{\mathbf{k}_{2}}(\omega_{2})$ (25)

which only contributes if the quantum numbers of the photodetected electrons
are not identical due to Pauli exclusion. This is a useful check to determine
the overall magnitude of the pair ARPES compared to single particle ARPES, and
can help to assess the spectral intensities of two-particle collective modes
separately from the single particle continuum.

###  III.2 BCS theory

If the system of interest is well-described by a BCS mean field ansatz, the
valence band is again composed of free Bogoliubov fermions. In this case, the
four-point function factorizes, and the coincidence pair photoemission signal
additionally includes a pairing term:

$\displaystyle
D_{\mathbf{k}_{1}\mathbf{k}_{2}}^{(0)}(\omega_{1},\omega_{2})=P_{\mathbf{k}_{1}}(\omega_{1})P_{\mathbf{k}_{2}}(\omega_{2})+\left|P_{\mathbf{k}_{1}\mathbf{k}_{2}}^{\textrm{pair}}(\omega_{1},\omega_{2})\right|^{2}$
(26)

where

$\displaystyle
P^{\textrm{pair}}_{\mathbf{k}_{1}\mathbf{k}_{2}}(\omega_{1,2})=\int d\tau
d\tau^{\prime}s(\tau)s(\tau^{\prime})e^{i(\omega_{1}\tau+\omega_{2}\tau^{\prime})}\langle\mathcal{T}\hat{c}_{\mathbf{k}_{1}}^{{\dagger}}(\tau)\hat{c}_{\mathbf{k}_{2}}^{{\dagger}}(\tau^{\prime})\rangle$
(27)

is a weighted time average of the time-ordered anomalous Green’s function. For
the case where the shape functions $s(t)=1$ the BCS singlet pair wavefunction
gives the value

$\displaystyle P^{\textrm{pair}}_{\mathbf{k}_{1}\mathbf{k}_{2}}(\omega_{1,2})$
$\displaystyle=\delta(\mathbf{k}_{1}+\mathbf{k}_{2})\delta(\sigma_{1}+\sigma_{2})\delta(\omega_{1}+\omega_{2})$
$\displaystyle\times\frac{\Delta_{\mathbf{k}_{1}}}{(\omega_{1}-i\eta)^{2}-E_{\mathbf{k}_{1}}^{2}}$
(28)

with the Bogoliubov energy given by
$E^{2}_{\mathbf{k}}=\epsilon_{\mathbf{k}}^{2}+\Delta^{2}_{\mathbf{k}}$ for
free particle dispersion $\epsilon_{\mathbf{k}}$, and $\eta$ is a small real
quantity [6].

We note that Eq. (28) yields a sharp peak at the Fermi level
($\omega_{1,2}=0$) when the delta functions are satisfied, indicating that
pair ARPES can be used to detect the underlying Cooper pair structure in terms
of center of mass spin (i.e., singlet versus triplet) and momentum (i.e.,
Fulde-Ferrel or pair density-wave state) as has been noted previously [6, 8].
Moreover, the fermion momentum dependence of the energy gap
$\Delta(\mathbf{k})$ can be scanned and directly measured.

###  III.3 Hubbard models for correlated electrons

The single band Hubbard model may provide a simple way to characterize the
behavior of pair photoemission for correlated electrons in systems without
superconducting long-range order. Specifically we will utilize eigenstates of
the particle-hole symmetric Hubbard model

$\displaystyle H=-t\sum_{\langle
i,j\rangle,\sigma}c^{\dagger}_{i,\sigma}c_{j,\sigma}+U\sum_{i}(n_{i,\uparrow}-\frac{1}{2})(n_{i,\downarrow}-\frac{1}{2})$
(29)

on an 8A (diamond) Betts cluster [12] to construct pair ARPES [Fig. 2(a)].
Here $c_{i,\sigma},c^{\dagger}_{i,\sigma}$ removes, adds a particle at site
$i$ with spin $\sigma$, $n_{i,\sigma}$ is the particle density per spin at
site $i$, $t$ denotes hybridization between nearest neighbor sites $i$ and
$j$, and $U$ is a measure of the local interaction between opposite spins.
Throughout we assume units where $\hbar=1$.

While much work has been performed via density matrix group renormalization
techniques (DMRG) for example to ascertain whether the Hubbard model in the
thermodynamic limit harbors superconductivity, our goal is more modest. By
examining the eigenstates and constructing pair ARPES on finite clusters,
which cannot have a bona fide phase transition, we may be able to highlight
how coincidence spectroscopy can be used to quantitatively measure pair field
susceptibilities in systems where $U(1)$ gauge symmetry is not broken but
fluctuating order may be inferred.

Figure 2: Pair ARPES for the attractive and repulsive Hubbard model. (a)
Schematics of the eight-site Betts cluster. (b) Pair binding energy as a
function of repulsive and attractive Hubbard interactions at half filling. (c)
Single-particle spectrum $A(\omega)$ as a function of interactions and
momentum (blue, green, red correspond to
$\mathbf{k}=0,(\pi/2,\pi/2),(\pi,\pi)$, respectively). (d) Top and bottom row
panels show opposite-spin $D_{\mathbf{k},-\mathbf{k}}(\omega_{1},\omega_{2})$
for repulsive and attractive interactions, respectively, from $|U|=0.5$ (left)
to $|U|=3$ (right). Dashed lines (center and difference frequencies) are
guides to the eye. All energies are quoted in units of $t=1$. While the main
photoemission peak with $\omega_{1}=\omega_{2}$ identically tracks attractive
and repulsive Coulomb interactions, for $U<0$ pair ARPES reveals the pair-
breaking intermediate state via the frequency difference spectrum
$\omega_{1}=\omega_{2}$ for $\omega_{1}+\omega_{2}=0$ (bottom row).

Pairing has been long investigated in exact diagonalization studies of the
Hubbard model on small clusters [13, 14, 15]. The pair binding energy $\Delta$
is defined as the energy difference between the ground state energies of $N$
and $N-2$ particle systems minus twice the energy of the $N-1$ system:

$\Delta=E_{N}+E_{N-2}-2E_{N-1}$ (30)

A negative $\Delta$ indicates an effective electron pair attraction.

The pair binding energy $\Delta$ obtained for the repulsive and attractive
Hubbard model at half-filling $N_{electrons}=8=N$ is shown in Fig. 2(b). For
repulsive $U$, $\Delta$ is negative for $U/t\lesssim 8$ and becomes positive
for larger values. The ground state of the attractive Hubbard model ($U<0$)
can be well modeled as a BCS superconducting paired state [15], and possesses
a pair binding energy that increases with $\mid U\mid$.

The ARPES spectra are identical for repulsive or attractive $U$ via particle-
hole symmetry. Fig. 2(b) depicts the spectral functions as a function of
$U/t$, with peaks corresponding to the three unique momenta on the 8A Betts
cluster $(0,0),(\pi/2,\pi/2)$ (6-fold degenerate), and $(\pi,\pi)$, shown in
blue, green and red. As noted previously [11], spectral peaks move to deeper
binding energies as $\mid U\mid$ is increased, and the development of the
lower Hubbard band can be more clearly observed. While a pairing gap is
clearly observable for attractive $U$, a superconductor cannot be
distinguished from a Mott gap in single-electron ARPES. Indeed, the spectra
for positive and negative $U$ are identical by virtue of particle hole
symmetry. This further motivates investigating pair photoemission, which
intrinsically discriminates between pair and density excitations.

Figure 3: Comparing Single and Pair Photoemission. (a) and (b) depict Pair
ARPES $D_{\mathbf{k},-\mathbf{k}}(\omega_{1},\omega_{2})$ and uncorrelated
pairs of single photoemission events $P_{\mathbf{k}}(\omega_{1})\times
P_{-\mathbf{k}}(\omega_{2})$, respectively, with line cuts for center and
difference frequencies shown in (c). Depicted responses are computed for
attractive interactions $U/t=-3$ with shape function broadening $\sigma=4/t$.

By inspection of the denominators in Eq. (20), one can expect that for pair
ARPES the largest intensity will be obtained for a given ${\bf k_{1},k_{2}}$
when the energies $\omega_{1},\omega_{2}$ are tuned to the respective energy
positions of ARPES removal spectra, giving roughly the similar pattern as that
obtained by simply multiplying the two independent ARPES spectral functions
for photoemitted electrons with opposite spin.

We focus only on momentum states lying closest to the chemical potential and
consider two-particle removal ARPES spectra for opposite spins and momenta
${\bf k_{1,2}}$ drawn from the six degenerate momentum points
$(\pm\pi/2,\pm\pi/2),(\pi,0),$ and $(0,\pi)$. Obtaining the eigenstates for
the sectors containing $N_{electrons}=8,7,6$ allows for the construction of
Fermi golden rule pair ARPES spectral functions $D_{\bf
k_{1},k_{2}}(\omega_{1},\omega_{2})$ via Eqs. (19) and (20), or equivalently
via Eq. (35) upon inclusion of the probe shape functions. We focus on spin-
resolved pair photoemission spectra; the spin-agnostic response follows via
summing equal- and opposite-spin contributions.

The resulting pair photoemission spectra are shown in Fig. 2(d) for
$\mathbf{k}_{1}=-\mathbf{k}_{2}=(\pi/2,\pi/2)$ and opposite spins
$\sigma_{1}=\uparrow,\sigma_{2}=\downarrow$, for both repulsive and attractive
interactions. While both cases show a primary peak at equal pair photoemission
energies $\omega_{1}=\omega_{2}$, corresponding to the particle-hole-symmetric
Hubbard gap, a key new feature is the emergence of a pair of additional peaks
for the attractive Hubbard model, with $\omega_{1}+\omega_{2}=0$. These
directly probe the pair breaking intermediate state and can be understood as a
two-step process: First, a photon breaks a Cooper pair to photoemit an
electron, while leaving an unpaired electron with pair breaking energy
$2\Delta$ in the sample. (2) Subsequently the second photon photoemits this
unpaired electron while removing the extra intermediate state energy from the
sample. As final state with two electrons removed from a fully-paired
superconductor has the same energy as the initial state, the total energy
$\omega_{1}+\omega_{2}$ left in the sample by the photoemission process must
equal to zero; the intermediate pair breaking state remains encoded in the
energy difference $\omega_{1}-\omega_{2}$. These observations can be confirmed
by comparing the pair photoemission response to uncorrelated pairs of single
photoemission processes, depicted in Fig. 3.

Figure 4: Spin and momentum dependence of pair ARPES. Top (a)-(d) rows depict
the pair ARPES response for for the attractive Hubbard model $U/t=-3$, for
equal/opposite spin and momentum combinations of the momenta
$(\pm\pi/2,\pm\pi/2)$ used in Fig. 2. Bottom rows (e)-(h) depict the
subdominant pair ARPES response for other momenta (scaled by a factor 1000
with respect to the top row), for equal and opposite emitted spins.

The momentum dependence on each of the fermion momenta as well as the net
total momentum ${\bf q}={\bf k_{1}}+{\bf k_{2}}$ and net spin
$\sigma=\sigma_{1}+\sigma_{2}$ can reveal further information of the pair wave
function. Fig. (4) plots pair ARPES for different combinations of photoemitted
wavevectors and spins, for $U/t=3$, as a function of center
$\omega_{1}+\omega_{2}$ and relative $\omega_{1}-\omega_{2}$ frequencies. As
expected for singlet pairing in the attractive Hubbard model, one immediately
finds that equal-spin photoemitted electrons [Fig. 4(b)] lack the pair
breaking peaks at zero center frequency of the opposite-spin response in Fig.
2(d). In contrast, observations of equal-spin pair breaking peaks in the
correlated pair ARPES response would be suggestive of triplet pairing
instabilities. A similar argument follows for photoemitted pairs of electrons
with equal momentum $\mathbf{k}_{1}=\mathbf{k}_{2}$, as shown in Fig. 4(c) and
(d); here, an observation of zero-frequency side peaks would be indicative of
finite-momentum pairing [16]. Pair photoemission for other momenta remains
strongly suppressed [Fig. 4(e)-(h)] for attractive interactions.

It is expected that these results will be affected by the finite size and
geometry of the small cluster, as well as adding symmetry breaking terms, such
as $t^{\prime}$, that can break momentum degeneracies. For example, the pair-
field correlator obtained for the same Hamiltonian on a $4\times 2$ cluster
that breaks $C4$ symmetry, increasing the number of non-degenerate momentum
points from 3 in the 8A cluster to 6, has quantitatively the same results for
attractive and repulsive $\mid U\mid=4t$. The largest low frequency
contribution is for pair momenta ${\bf q}=(\pi,0)$ and $(\pi/2,\pi)$. By
including a negative next nearest hopping $t^{\prime}=-0.25t$, the low energy
pair field correlations are largest for ${\bf q}=(\pi,0)$ and $(0,\pi)$ for
$U=4t$, while for $U=-4t$, ${\bf q}=(0,0)$ still is largest. These effects are
larger clusters and different geometries should be further addressed.

Lastly, here we have restricted consideration to zero temperature pair ARPES.
One key application of pair ARPES could be to approach ordered phases from
high temperature to measure how pair field correlations develop, either
through towards a true superconducting transition, or averted by the onset of
another competing order, such as charge and/or spin density waves. As these
phases all appear to have nearly the same ground state energies in simulations
of the Hubbard model, an experimental investigation may provide finer insight
into which terms may be missing from the Hubbard model that could formulate a
closer contact to materials such as the high temperature superconductors.

In summary, we have presented a theory for pair ARPES whereby two photons
produce two photoelectrons detected in coincidence, resolved in both energy
and momentum. The corresponding two-particle removal spectra can thus be
exploited to determine the effects of electron correlations in a direct way.
The calculated pair response for the attractive and repulsive Hubbard model at
half filling for an 8 site Betts cluster shows spectroscopically how prominent
ordering tendencies of superconductivity and the net pair momentum and spin
can be inferred directly from experiments.

Experimentally, to make data as close to superconducting pair correlation
function as possible, one should try to eject the electron pair from the
sample at the same time. In such an “instantaneous event”, two photons eject
two electrons in an “interacting volume”, for examples a Cooper pair in a
superconductor, or a pair in a Mott insulator that are sufficiently entangled.
In a Cooper pair, this means two electrons within the superconducting
coherence volume. In a Mott insulator, assuming that Hubbard model is a
reasonable starting point, this means two electrons not far from each other so
that a cascade of local interactions can entangle the electrons. Our
theoretical calculation was carried out in a small cluster such that the
entanglement is naturally strong.

Such pair photoemission is an experiment with many technical challenges.
However, several recent technological advances make it realistic. The first is
the emergence of much improved and suitable light sources, such as UV lasers,
high harmonic generation, free electron lasers, and photon focusing schemes.
Photons within a very short pulse, such as tens of femtoseconds, can be
considered as identical and instantaneous within time of flight spectrometers
having picosecond resolution. The second is the time-of-flight (TOF) based
three dimensional ARPES platform, such as the momentum microscope and its spin
filtered variant. The third is the development of two-dimensional ultrafast
multichannel detectors. With time and through an integration of these
important new technologies, enhanced by timing, energy, momentum
discrimination schemes and machine learning algorithms to improve the signal
to noise ratio, this new spectroscopy may be developed in the near future.

###### Acknowledgements.

Authors would like to thank Frank Marsiglio and Joseph Orenstein for
insightful discussions. MZ, JEM, DM, FM, PA, ZXS and TPD acknowledge support
for the work from the U.S. Department of Energy (DOE), Office of Basic Energy
Sciences, Division of Materials Sciences and Engineering, and through Contract
No. DE-AC02-05-CH11231 via the Quantum Materials program (KC2202) (JEM).

## References

  * [1] Jonathan A. Sobota, Yu He, and Zhi-Xun Shen Rev. Mod. Phys. 93, 025006 (2021).
  * [2] J. Berakdar, Phys. Rev. B 58, 9808 (1998).
  * [3] G. A. Sawatzky, Treatise on Materials Science & Technology, 30, 167 (1988).
  * [4] G.A. van Riessen, Z. Wei, R.S. Dhaka, C. Winkler, F.O. Schumann and J. Kirschner, J. Phys.: Condens. Matter 22, 092201 (2010).
  * [5] J. K. Freericks, H. R. Krishnamurthy, and Th. Pruschke, Theoretical Description of Time-Resolved Photoemission Spectroscopy: Application to Pump-Probe Experiments, Phys. Rev. Lett. 102, 136401 (2009).
  * [6] Y. Su, and C. Zhang, Coincidence angle-resolved photoemission spectroscopy: Proposal for detection of two-particle correlations, Phys. Rev. B 101, 205110 (2020).
  * [7] Luuk J. P. Ament, Michel van Veenendaal, Thomas P. Devereaux, John P. Hill, and Jeroen van den Brink, Rev. Mod. Phys. 83, 705 (2011).
  * [8] Fahad Mahmood, Thomas Devereaux, Peter Abbamonte, and Dirk K. Morr, Phys. Rev. B 105, 064515 (2022).
  * [9] T. A. Maier and D. J. Scalapino, npj Quantum Mater. 4, 30 (2019).
  * [10] Edwin W. Huang, Wen O. Wang, Jixun K. Ding, Tianyi Liu, Fangze Liu, Xu-Xin Huang, Brian Moritz, and Thomas P. Devereaux, J. Phys. Soc. Jpn. 90, 111010 (2021).
  * [11] E. Dagotto, R. Joynt, A. Moreo, S. Bacci, and E. Gagliano, Phys. Rev. B 41, 9049 (1990); E. Dagotto, A. Moreo, F. Ortolani, J. Riera, and D. J. Scalapino, Phys. Rev. Lett. 67, 1918 (1991); E. Dagotto, F. Ortolani, and D. J. Scalapino, Phys. Rev. B 46, 3183 (1992).
  * [12] D. Betts, H. Lin, and J. Flynn, Can. J. Phys. 77, 353 (1999).
  * [13] H. Q. Lin, J. E. Hirsch, and D. J. Scalapino, Phys. Rev. B 37, 7359 (1988).
  * [14] J. A. Riera and A. P. Young, Phys. Rev. B 39, 9697 (1989).
  * [15] F. Marsiglio, Phys. Rev. B 55, 575 (1997).
  * [16] See, e.g., Daniel F. Agterberg et al., Annual Review of Condensed Matter Physics 11, 231 (2020).
  * [17] J. T. Anderson and A. M. Goldman, Phys. Rev. Lett. 25, 743 (1970); N. Bergeal, J. Lesueur, M. Aprili, G. Faini, J. P. Contour, and B. Leridon, Nature Physics 4, 608 (2008).

##  Appendix A Two-Photon Pair Photoemission

The main text presents a formal expression of the pair photoemission response
which assumes negligible backaction between photoemitted electrons and sample
electrons. Starting from the second-order perturbative expression

$\displaystyle D_{\mathds{k}_{1}\mathds{k}_{2}}$
$\displaystyle=\int\displaylimits_{-\infty}^{t}d\tau_{1}d\tau_{2}\int\displaylimits_{-\infty}^{\tau_{1}}d\tau^{\prime}_{1}\int\displaylimits_{-\infty}^{\tau_{2}}d\tau^{\prime}_{2}~{}s(\tau_{1})s(\tau^{\prime}_{1})s(\tau_{2})s(\tau^{\prime}_{2})\int
d\mathbf{r}_{1}d\mathbf{r}^{\prime}_{1}d\mathbf{r}_{2}d\mathbf{r}^{\prime}_{2}\int
d\mathbf{x}_{1}d\mathbf{x}^{\prime}_{1}d\mathbf{x}_{2}d\mathbf{x}^{\prime}_{2}~{}\phi_{\mathds{k}_{1}}^{\star}(\mathbf{x}_{1})\phi_{\mathds{k}_{2}}^{\star}(\mathbf{x}^{\prime}_{1})\phi_{\mathds{k}_{2}}(\mathbf{x}^{\prime}_{2})\phi_{\mathds{k}_{1}}(\mathbf{x}_{2})$
$\displaystyle\times
e^{i\omega_{\textrm{ph}}(\tau_{1}+\tau^{\prime}_{1}-\tau_{2}-\tau^{\prime}_{2})}~{}g^{\star}(\mathbf{r}_{1})g^{\star}(\mathbf{r}^{\prime}_{1})g(\mathbf{r}^{\prime}_{2})g(\mathbf{r}_{2})\sum_{\sigma_{1}\sigma^{\prime}_{1}\nu}\sum_{\sigma_{2}\sigma^{\prime}_{2}\nu^{\prime}}\left<\hat{\Phi}_{\sigma_{1}}^{\dagger}(\mathbf{r}_{1},\tau_{1})\hat{\Psi}_{\sigma_{1}}(\mathbf{r}_{1},\tau_{1})\hat{\Phi}_{\sigma^{\prime}_{1}}^{\dagger}(\mathbf{r}^{\prime}_{1},\tau^{\prime}_{1})\hat{\Psi}_{\sigma^{\prime}_{1}}(\mathbf{r}^{\prime}_{1},\tau^{\prime}_{1})\right.$
$\displaystyle\times\left.\hat{\Psi}_{\nu}^{\dagger}(\mathbf{x}_{1},t)\hat{\Psi}_{\nu^{\prime}}^{\dagger}(\mathbf{x}^{\prime}_{1},t)\hat{\Psi}_{\nu^{\prime}}(\mathbf{x}^{\prime}_{2},t)\hat{\Psi}_{\nu}(\mathbf{x}_{2},t)\hat{\Psi}_{\sigma^{\prime}_{2}}^{\dagger}(\mathbf{r}^{\prime}_{2},\tau^{\prime}_{2})\hat{\Phi}_{\sigma^{\prime}_{2}}(\mathbf{r}^{\prime}_{2},\tau^{\prime}_{2})\hat{\Psi}_{\sigma_{2}}^{\dagger}(\mathbf{r}_{2},\tau_{2})\hat{\Phi}_{\sigma_{2}}(\mathbf{r}_{2},\tau_{2})\right>$
(31)

Neglecting backaction on the sample permits a decomposition

$\displaystyle D_{\mathds{k}_{1}\mathds{k}_{2}}$
$\displaystyle=\sum_{\sigma_{1}\sigma^{\prime}_{1}\nu}\sum_{\sigma_{2}\sigma^{\prime}_{2}\nu^{\prime}}\int\displaylimits_{-\infty}^{t}d\tau_{1}d\tau_{2}\int\displaylimits_{-\infty}^{\tau_{1}}d\tau^{\prime}_{1}\int\displaylimits_{-\infty}^{\tau_{2}}d\tau^{\prime}_{2}~{}s(\tau_{1})s(\tau^{\prime}_{1})s(\tau_{2})s(\tau^{\prime}_{2})e^{i\omega_{\textrm{ph}}(\tau_{1}+\tau^{\prime}_{1}-\tau_{2}-\tau^{\prime}_{2})}$
$\displaystyle\times\int
d\mathbf{r}_{1}d\mathbf{r}^{\prime}_{1}d\mathbf{r}_{2}d\mathbf{r}^{\prime}_{2}~{}g^{\star}(\mathbf{r}_{1})g^{\star}(\mathbf{r}^{\prime}_{1})g(\mathbf{r}^{\prime}_{2})g(\mathbf{r}_{2})\left<\hat{\Phi}_{\sigma_{1}}^{\dagger}(\mathbf{r}_{1},\tau_{1})\hat{\Phi}_{\sigma^{\prime}_{1}}^{\dagger}(\mathbf{r}^{\prime}_{1},\tau^{\prime}_{1})\hat{\Phi}_{\sigma^{\prime}_{2}}(\mathbf{r}^{\prime}_{2},\tau^{\prime}_{2})\hat{\Phi}_{\sigma_{2}}(\mathbf{r}_{2},\tau_{2})\right>$
$\displaystyle\times\int
d\mathbf{x}_{1}d\mathbf{x}^{\prime}_{1}~{}\phi_{\mathds{k}_{1}}^{\star}(\mathbf{x}_{1})\phi_{\mathds{k}_{2}}^{\star}(\mathbf{x}^{\prime}_{1})~{}\left<\hat{\Psi}_{\sigma_{1}}(\mathbf{r}_{1},\tau_{1})\hat{\Psi}_{\sigma^{\prime}_{1}}(\mathbf{r}^{\prime}_{1},\tau^{\prime}_{1})\hat{\Psi}_{\nu}^{\dagger}(\mathbf{x}_{1},t)\hat{\Psi}_{\nu^{\prime}}^{\dagger}(\mathbf{x}^{\prime}_{1},t)\right>$
$\displaystyle\times\int
d\mathbf{x}_{2}d\mathbf{x}^{\prime}_{2}~{}\phi_{\mathds{k}_{2}}(\mathbf{x}^{\prime}_{2})\phi_{\mathds{k}_{1}}(\mathbf{x}_{2})~{}\left<\hat{\Psi}_{\nu^{\prime}}(\mathbf{x}^{\prime}_{2},t)\hat{\Psi}_{\nu}(\mathbf{x}_{2},t)\hat{\Psi}_{\sigma^{\prime}_{2}}^{\dagger}(\mathbf{r}^{\prime}_{2},\tau^{\prime}_{2})\hat{\Psi}_{\sigma_{2}}^{\dagger}(\mathbf{r}_{2},\tau_{2})\right>$
(32)

Suppose now that the sample electrons near the Fermi energy are confined to a
single valence band
$\hat{\Phi}_{\sigma}(\mathbf{r},\tau)\approx\frac{1}{\sqrt{L}}\sum_{\bar{\mathbf{k}}}u_{\bar{\mathbf{k}}}(\mathbf{r})\hat{c}_{\bar{\mathbf{k}}\sigma}^{\vphantom{{\dagger}}}(\tau)e^{i\bar{\mathbf{k}}\bar{\mathbf{r}}}$
with Bloch function $u_{\bar{\mathbf{k}}}(\mathbf{r})$. Expanding the photo-
emitted electron fields in a plane wave basis, and Fourier transforming the
detector form factors $\phi_{\mathds{k}}(\mathbf{x})=\int
d\mathbf{p}e^{-i\mathbf{p}\mathbf{x}}\phi_{\mathds{k}}(\mathbf{p})$, one
obtains

$\displaystyle D_{\mathds{k}_{1}\mathds{k}_{2}}$
$\displaystyle=\sum_{\sigma_{1}\sigma^{\prime}_{1}\nu}\sum_{\sigma_{2}\sigma^{\prime}_{2}\nu^{\prime}}\int\displaylimits_{-\infty}^{t}d\tau_{1}d\tau_{2}\int\displaylimits_{-\infty}^{\tau_{1}}d\tau^{\prime}_{1}\int\displaylimits_{-\infty}^{\tau_{2}}d\tau^{\prime}_{2}~{}s(\tau_{1})s(\tau^{\prime}_{1})s(\tau_{2})s(\tau^{\prime}_{2})e^{i\omega_{\textrm{ph}}(\tau_{1}+\tau^{\prime}_{1}-\tau_{2}-\tau^{\prime}_{2})}$
$\displaystyle\times\int
d\mathbf{k}d\mathbf{k}^{\prime}d\mathbf{q}~{}\left<\hat{c}_{\bar{\mathbf{k}},\sigma_{1}}^{{\dagger}}(\tau_{1})\hat{c}_{\bar{\mathbf{q}}-\bar{\mathbf{k}},\sigma^{\prime}_{1}}^{{\dagger}}(\tau^{\prime}_{1})\hat{c}_{\bar{\mathbf{q}}-\bar{\mathbf{k}}^{\prime},\sigma^{\prime}_{2}}^{\vphantom{{\dagger}}}(\tau^{\prime}_{2})\hat{c}_{\bar{\mathbf{k}}^{\prime},\sigma_{2}}^{\vphantom{{\dagger}}}(\tau_{2})\right>$
$\displaystyle\times\int_{\textrm{BZ}}d\mathbf{r}_{1}d\mathbf{r}^{\prime}_{1}d\mathbf{r}_{2}d\mathbf{r}^{\prime}_{2}~{}g^{\star}(\mathbf{r}_{1})g^{\star}(\mathbf{r}^{\prime}_{1})g(\mathbf{r}^{\prime}_{2})g(\mathbf{r}_{2})~{}u^{*}_{\bar{\mathbf{k}}}(\mathbf{r}_{1})u^{*}_{\bar{\mathbf{q}}-\bar{\mathbf{k}}}(\mathbf{r}^{\prime}_{1})u_{\bar{\mathbf{q}}-\bar{\mathbf{k}}^{\prime}}(\mathbf{r}^{\prime}_{2})u_{\bar{\mathbf{k}}^{\prime}}(\mathbf{r}_{2})$
$\displaystyle\times\int
d\mathbf{p}~{}\phi_{\mathds{k}_{1}}^{\star}(\mathbf{q}-\mathbf{p})\phi_{\mathds{k}_{2}}^{\star}(\mathbf{p})~{}\left<\hat{\Psi}_{\sigma_{1}}(\mathbf{k},\tau_{1})\hat{\Psi}_{\sigma^{\prime}_{1}}(\mathbf{q}-\mathbf{k},\tau^{\prime}_{1})\hat{\Psi}_{\nu}^{\dagger}(\mathbf{q}-\mathbf{p},t)\hat{\Psi}_{\nu^{\prime}}^{\dagger}(\mathbf{p},t)\right>$
$\displaystyle\times\int
d\mathbf{p}^{\prime}~{}\phi_{\mathds{k}_{2}}(\mathbf{p}^{\prime})\phi_{\mathds{k}_{1}}(\mathbf{q}-\mathbf{p}^{\prime})~{}\left<\hat{\Psi}_{\nu^{\prime}}(\mathbf{p}^{\prime},t)\hat{\Psi}_{\nu}(\mathbf{q}-\mathbf{p}^{\prime},t)\hat{\Psi}_{\sigma^{\prime}_{2}}^{\dagger}(\mathbf{q}-\mathbf{k}^{\prime},\tau^{\prime}_{2})\hat{\Psi}_{\sigma_{2}}^{\dagger}(\mathbf{k}^{\prime},\tau_{2})\right>$
(33)

Rewriting the third line in terms of matrix elements recovers the expression
presented in the main text.

##  Appendix B Spectral representation of Pair ARPES

A useful representation of pair ARPES that accounts for the probe shape
function follows from a spectral decomposition of Eq. (18)

$\displaystyle
D^{(0)}_{k_{1},k_{2}}(\omega_{1},\omega_{2})=\sum_{n}\left|\int\displaylimits_{-\infty}^{\infty}dt\int\displaylimits_{0}^{\infty}dt^{\prime}\right.~{}s(t)s(t-t^{\prime})$
$\displaystyle\times~{}\sum_{m}\left[\vphantom{\frac{1}{2}}\right.\left<n\right|\hat{c}_{\mathbf{k}_{1}}^{\vphantom{{\dagger}}}\left|m\right>\left<m\right|{\hat{c}_{\mathbf{k}_{2}}^{\vphantom{{\dagger}}}}\left|0\right>e^{i\omega_{1}t^{\prime}}-\left.\left<n\right|\hat{c}_{\mathbf{k}_{2}}^{\vphantom{{\dagger}}}\left|m\right>\left<m\right|\hat{c}_{\mathbf{k}_{1}}^{\vphantom{{\dagger}}}\left|0\right>e^{i\omega_{2}t^{\prime}}\right]\left.e^{i[(\epsilon_{n}-\epsilon_{0}-\omega_{1}-\omega_{2})t-(\epsilon_{n}-\epsilon_{m})t^{\prime}]}\right|^{2}$
(34)

with $s(t)$ the Gaussian shape functions. The time integrals can now be
evaluated and one arrives at

$D^{(0)}_{k_{1},k_{2}}(\omega_{1},\omega_{2})=\sum_{n}\left|\sum_{m}\left[\left<n\right|\hat{c}_{\mathbf{k}_{1}}^{\vphantom{{\dagger}}}\left|m\right>\left<m\right|{\hat{c}_{\mathbf{k}_{2}}^{\vphantom{{\dagger}}}}\left|0\right>X_{nm}(\omega_{1},\omega_{2})\right.\right.-~{}\left<n\right|\left.\left.\hat{c}_{\mathbf{k}_{2}}^{\vphantom{{\dagger}}}\left|m\right>\left<m\right|\hat{c}_{\mathbf{k}_{1}}^{\vphantom{{\dagger}}}\left|0\right>X_{nm}(\omega_{2},\omega_{1})\right]\vphantom{\frac{1}{2}}\right|^{2}$
(35)

where

$X_{nm}(\omega_{1},\omega_{2})=\sigma^{2}e^{-\frac{\sigma^{2}}{4}\left(\omega_{1}+\omega_{2}+\epsilon_{0}-\epsilon_{n}\right)^{2}}\left\\{\frac{1}{2}e^{-\frac{\sigma^{2}}{4}\left(\omega_{2}-\omega_{1}+\epsilon_{0}+\epsilon_{n}-2\epsilon_{m}\right)^{2}}+\frac{i}{\sqrt{\pi}}\mathcal{D}\left[\frac{\sigma}{2}\left(\omega_{2}-\omega_{1}+\epsilon_{0}+\epsilon_{n}-2\epsilon_{m}\right)\right]\right\\}$
(36)

Here, $\sigma$ is the Gaussian broadening of the shape function and
$\mathcal{D}(x)$ denotes the Dawson function, which is related to the error
function $\mathcal{D}(x)=\sqrt{\pi/4}e^{-x^{2}}\textrm{erfi}(x)$.