# Nonlinear quantum interferometric spectroscopy with entangled photon pairs

Shahaf Asban<EMAIL_ADDRESS>Department of Chemistry and Physics & Astronomy,
University of California, Irvine, California 92697-2025, USA  Vladimir Y.
Chernyak<EMAIL_ADDRESS>Department of Chemistry, Wayne State
University, 5101 Cass Ave, Detroit, Michigan 48202, USA  Department of
Mathematics, Wayne State University, 656 W. Kirby, Detroit, Michigan 48202,
USA  Shaul Mukamel<EMAIL_ADDRESS>Department of Chemistry and Physics &
Astronomy, University of California, Irvine, California 92697-2025, USA

###### Abstract

We develop closed expressions for a time-resolved photon counting signal
induced by an entangled photon pair in an interferometric spectroscopy setup.
Superoperator expressions in Liouville-space are derived that can account for
relaxation and dephasing induced by coupling to a bath. Interferometric setups
mix matter and light variables non-trivially, which complicates their
interpretation. We provide an intuitive modular framework for this setup that
simplifies its description. Based on separation between the detection stage
and the light-matter interaction processes. We show that the pair entanglement
time and the interferometric time-variables control the observed physics time-
scale. Only a few processes contribute in the limiting case of small
entanglement time with respect to the sample response, and specific
contributions can be singled out.

###  I Introduction

Interferometric setups introduce a promising new methodology for quantum
inference of matter information in quantum spectroscopy [1, 2, 3, 4, 5]. Since
quantum probes change their state upon interaction with external systems,
multiphoton coincidence detection schemes should reveal quantum correlations
induced by the sample [6, 7]. We consider the optical measurement setup shown
in Fig. 1 – that includes linear and nonlinear elements that transform the
optical field before and after it is coupled to a sample. We focus on time-
resolved detection of two photons. We use one interferometric setup to prepare
the initial quantum state of the light, followed by a second interferometer
that manipulates the arrival times of matter induced radiation to obtain
matter pathway resolution. Interferometric schemes such as the Mach-Zehnder
[8], Hong-Ou-Mandel [9], and Franson [10] provide a useful toolbox to scan the
change in photon statistics by coupling to a sample from which matter
information may be inferred using multiphoton detection in coincidence [11,
10, 12]. Quantum-enhancements of interferometric detection schemes are an
experimental reality in several fields [13]. The direct detection of
gravitational waves [14, 15, 16], unprecedented phase estimation precision
with high loss tolerance at lower photon flux [17, 18, 19, 20, 21], and in
wide-field imaging [22] are all contemporary examples. Here, we introduce a
new family of signals applied to the inference of objects at the microscopic
scale.

Liouville-space pathways break-down the density operator evolution into the
set of physical processes determined by time ordered excitation and de-
excitation processes (pathways) induced by the applied fields. Sorting them
out is important in order to develop understanding the underlying matter
dynamics. This description become crucially important when considering the
effects of the environment that break time reversal symmetry. Sorting these
pathways allows to infer the role of each process in a systematic manner.
Liouville pathway resolution enable to compare model-based theoretical
predictions with experiment [23].

Our goal is to sort out the Liouville-space pathways, by scanning
interferometric delays in the preparation and detection stages. We wish to
identify what type of information regarding the dynamics of the system and its
coupling to a bath can be inferred from these measurements. We consider two
interferometers, at the preparation and detection processes, as depicted in
Fig. 1 and further discussed in Sec. II. Two-port linear interferometers
induce transformations in the two-photon space [24, 25]. They generate
rotations in the basis of the electromagnetic field [1]. These transformations
offer a unique set of control knobs used in both the preparation and detection
stages. These control parameters provide novel spectral windows [2]. We show
that the interferometric control variables enable to conveniently scan the
temporal dynamics. We derive compact superoperator expressions for the time-
resolved coincidence-signal of two photons expanded in terms of Liouville-
space pathways. For simplicity, we consider two ideal detectors which is fast
in comparison to all relaxation times of the sample. This extends Glauber’s
celebrated theory of detection [26].

The general expressions are given in appendix B. For simplicity we have
derived an approximation for the limiting case of vanishing two-photon
temporal separation (the entanglement time $T_{e}$). The two photons then
arrive simultaneously relative to the observed dynamics timescale. This
corresponds to a vanishing moving time average of the response, thus sensitive
to the dynamics above $T_{e}$ which is in the femtosecond regime. At this time
scales we expect environment effects to be more pronounced which requires the
Liouville space approach. Moreover, using the novel time variables, we are
able to separate different processes in the evolution of the sample in time
domain, from which different transport mechanisms can be studied in detail
such as exciton-exciton scattering [27, 28, 29].

###  II The setup

Our setup shown in Fig. 1a contains two interferometers. An entangled pair is
prepared using a Michelson interferometer as introduced in [13, 30, 31] and
depicted in Fig. 1b. At the detection stage we employ a Hong-Ou-Mandel
interferometer. In this section we provide the theoretical framework required
for the inclusion of both interferometers.

Figure 1: The quantum interferometry setup.(a) The Michelson interferometer
${\cal U}$ generates a pair of photons with variable exchange-phase $\theta$,
used as a control for the photon-pair degree of distinguishability (see Eq.
2). The photons acquire the relative delay $s$. (b) State preparation protocol
in which two photons are generated with a possible relative delay, then
coupled to thee sample and finally detected in coincidence in a Hong-Ou-Mandel
interferometer (boxed area).

Each stage in the setup, introduces a superposition of fields which we denote
basis rotation or transformation. These platforms are presented in a modular
manner such that each stage is responsible for a well defined property (e.g.,
introducing a delay). Quantum mechanically, the stages are inseparable since
more than one absorption, emission and detection events occur in all possible
time orderings. However, we limit the discussion to a distant detection plane,
such that the propagation time of the emitted photons is much longer than the
typical time of the entire process. The detection and interaction are then
completely factorized temporally as depicted in Fig. 2

Figure 2: Block design for two-photon interferometric spectroscopy. Decoupling
of the interaction and detection process when both are well temporally
defined. (a) The detection pathways corresponding to HOM interferometer placed
after the sample. Here $\tau>0$ is the detection time difference of the two
photons and $T>0$ is the HOM delay. Two paths of direct propagation
$O_{I-II}$, along with $O_{III-IV}$ in which there is a coherence between
reflected and transmitted modes. (b) Ladder diagrams that correspond to the
light-matter interaction corresponding to the matter correlation functions in
Eq. 8a (for a three-level model system). Only processes in which two-photons
interact with the sample and two are detected and sketched.

####  II.1 Interferometric state preparation

We assume a broadband ultrafast pump pulse which is known to imprint
identifying spectral information on each of the optical modes. The modified
Michelson interferometer in Fig. 1a, generates control over this feature by
controlled systematization of the wavefunction in a type-II phase matching
setting. This results in engineered degree of distinguishability as discussed
below.

The pump beam is reflected by the dichroic mirror (DM) then passes the first
time through the BBO crystal. The generated entangled pair passes through the
$\nicefrac{{\lambda}}{{4}}$ plate that switches the polarization such that
$\phi\left(\omega_{a},\omega_{b}\right)\rightarrow\phi\left(\omega_{b},\omega_{a}\right)$,
then exits the interferometer from the left. Otherwise, the pump photon is
reflected by the second DM and passes through the BBO crystal for the second
time with a controlled phase introduced by the PZT device on which mirror M1
is positioned generating possibly an entangled pair. Finally the pump beam is
filtered out of the interferometer by the first DM.

The combined nonlinear interferometer transformations create a two photon
wavefunction of the form

$|\Psi_{\theta}\rangle=\int
d\omega_{a}d\omega_{b}\Phi_{\theta}\left(\omega_{a},\omega_{b}\right)a^{\dagger}\left(\omega_{a}\right)b^{\dagger}\left(\omega_{b}\right)|\text{vac}\rangle.$
(1)

The amplitude is given by

$\Phi_{\theta}\left(\omega_{a},\omega_{b}\right)=\frac{1}{\sqrt{2}}\left[\phi\left(\omega_{a},\omega_{b}\right)+e^{i\theta}\phi\left(\omega_{b},\omega_{a}\right)\right].$
(2)

The joint spectral amplitude (JSA) resulting from the direct-channel
$\phi\left(\omega_{a},\omega_{b}\right)$ (due to one pass in the BBO with
aligned polarizations) used in our calculations is given by
$\phi\left(\omega_{a},\omega_{b}\right)={\cal
A}_{p}\left(\omega_{a}+\omega_{b}\right)\varphi\left(\omega_{a},\omega_{b}\right)$.
The Gaussian pump envelope ${\cal
A}_{p}\left(\omega\right)=\exp\left[\nicefrac{{\left(\omega-\omega_{p}\right)^{2}}}{{\sigma_{p}^{2}}}\right]$
is centered around $\omega_{p}$ characterized with the bandwidth $\sigma_{p}$
[32]. The phase-matching factor
$\varphi\left(\omega_{a},\omega_{b}\right)=\text{sinc}\left[\left(\omega_{a}-\bar{\omega}_{a}\right)T_{a}+\left(\omega_{b}-\bar{\omega}_{b}\right)T_{b}\right]$,
breaks the frequency exchange symmetry, i.e.,
$\varphi\left(\omega_{a},\omega_{b}\right)\neq\varphi\left(\omega_{b},\omega_{a}\right)$.
Here $\bar{\omega}_{a/b}$ is the central frequency of the signal and idler
beams, and $T_{a/b}=L\left(v_{a/b}^{-1}-v_{p}^{-1}\right)$, where $L$ is the
nonlinear crystal length and $v$ is the inverse group velocity at the relevant
central frequency $\left(\bar{\omega}_{a/b},\omega_{p}\right)$. In this (type-
II) phase matching condition, the channels are flipped using the their
opposite and orthogonal polarization degree of freedom, e.g.,
$|\Psi_{\theta}\rangle=|HV\rangle+e^{i\theta}|VH\rangle$ where $V$ and $H$
correspond to horizontal and vertical polarizations respectively. For
different phase matching condition which give rise to identiclly polarized
biphotons, one would expect a different output, e.g.,
$|\Psi_{\theta}\rangle=|HH\rangle+e^{i\theta}|VV\rangle$.

Here, we have calculated the preparation step in the Schödinger picture,
modifying the initial amplitude rather than the fields operators. The
detection stage is computed using the Heisenberg picture as explained in Sec.
III. This way, the dynamics is calculated in the matter’s reference-frame, as
done in [2]. Each interferometer introduces additional spectroscopic control
parameters that can be used to study the joint light-matter quantum state.

###  III Detection pathways

The detection process involves a HOM interferometer as shown in the boxed area
of Fig. 1b. We define the signal in time domain. We consider a sample
described by the Hamiltonian $H_{\mu}$, that is coupled to field degrees of
freedom $H_{\phi}$ by the dipolar interaction
$H_{\text{int}}=\boldsymbol{\mu}\left(\boldsymbol{r},t\right)\cdot\boldsymbol{E}\left(\boldsymbol{r},t\right)$.
Here
$\boldsymbol{\mu}\left(\boldsymbol{r},t\right)=\boldsymbol{V}\left(\boldsymbol{r},t\right)+\boldsymbol{V}^{\dagger}\left(\boldsymbol{r},t\right)$
is the dipole operator, $\boldsymbol{V}=\sum_{i<j}\mu_{ij}|i\rangle\langle j|$
is a lowering transition operator acting in the molecular Hilbert space with
the corresponding matrix element $\mu_{ij}$,
$\boldsymbol{E}\left(\boldsymbol{r},t\right)=\sum_{\sigma}\left[E_{\sigma}\left(\boldsymbol{r},t\right)+\text{H.c.
}\right]$ is the electric field operator given by
$E_{\sigma}\left(\boldsymbol{r},t\right)=\sum_{\boldsymbol{k}}\sqrt{\frac{2\pi
k}{\Omega_{Q}}}\hat{\epsilon}_{\sigma}\left(\boldsymbol{k}\right)a_{\boldsymbol{k},\sigma}\left(t\right)e^{i\boldsymbol{k}\cdot\boldsymbol{r}}$,
where $\hat{\epsilon}_{\sigma}\left(\boldsymbol{k}\right)$ is the
$\sigma$-polarization vector, $\Omega_{Q}$ is the quantization volume ($c=1$),
and
$a_{\boldsymbol{k},\sigma}\left(a_{\boldsymbol{k},\sigma}^{\dagger}\right)$
are (bosonic) photon annihilation (creation) operators obeying
$\left[a_{\boldsymbol{k},\sigma},a_{\boldsymbol{k}^{\prime},\sigma^{\prime}}^{\dagger}\right]=\delta_{\sigma,\sigma^{\prime}}\delta_{\boldsymbol{k},\boldsymbol{k}^{\prime}}$.
Hereafter we assume that the applied field is near resonance with a molecular
transition, such that the rotating wave approximation may be applied by
setting $H_{\text{int}}=EV^{\dagger}+E^{\dagger}V$.

The two-photon coincidence signal is defined by

$\displaystyle\hat{{\cal O}}\left(t,\tau\right)$
$\displaystyle=E_{a^{\prime},R}^{\dagger}\left(\boldsymbol{r}_{a},t\right)E_{b^{\prime},R}^{\dagger}\left(\boldsymbol{r}_{b},t+\tau\right)$
(3) $\displaystyle\times
E_{b^{\prime},L}\left(\boldsymbol{r}_{b},t+\tau\right)E_{a^{\prime},L}\left(\boldsymbol{r}_{a},t\right).$

Here, $E_{R}$ and $E_{L}$ are electric field superoperators that act from the
right $E_{R}\rho\equiv\rho E$, and left $E_{L}\rho\equiv E\rho$ of the density
operator. We further define their linear combinations $A_{\pm}O=A_{L}O\pm
A_{R}O$ (corresponding to Hilbert-space commutator and anti commutator). These
will be useful for the compact description in the interaction picture below.
The primes reflect the fact that the detection and interaction planes are
described in different basis sets due to the HOM transformation, as explained
below.

####  III.1 The HOM detection interferometer

In our setup in Fig. 1, two beams are combined on a beam splitter (BS) which
is mounted on a movable stage. This enables to scan variable propagation times
of the reflected pathways of the beam with respect to the transmitted ones,
introducing the HOM time-delay $T$. This process is described by a linear
transformation since the photons at the detection plane are represented in a
different basis. The transformation (Jordan-Schwinger map) can be represented
as an $SU\left(2\right)$ rotation in the frequency-domain [33, 34, 24, 35,
25], resulting in the input-output relation for the field operators. Writing
the field in vector notation given by
$\boldsymbol{E}\left(\boldsymbol{r},\omega\right)=\left(E_{a}\left(\boldsymbol{r},\omega\right),E_{b}\left(\boldsymbol{r},\omega\right)\right)^{T}$,
the HOM rotation the detected field is given by
$\boldsymbol{E}|_{\text{detection}}\left(\boldsymbol{r},\omega\right)=\hat{{\cal
R}}_{T}\boldsymbol{E}|_{\text{interaction}}\left(\boldsymbol{r},\omega\right)$
where $\hat{{\cal R}}_{T}$ is given by

$\hat{{\cal R}}_{T}=\left(\begin{array}[]{cc}t&ire^{i\omega T}\\\
ire^{-i\omega T}&t\end{array}\right).$ (4)

Here $t$ and $r$ are the transmission and reflection coefficients obeying
$\left|t\right|^{2}+\left|r\right|^{2}=1$. For the 50:50 BS considered here,
$t=r=\nicefrac{{1}}{{\sqrt{2}}}$. In the following, we express all field
operators in the matter interaction-domain basis
$\boldsymbol{E}|_{\text{interaction}}\equiv\boldsymbol{E}\left(\boldsymbol{r},\omega\right)$,
which requires the inverse rotation of the observable in Eq. 3 [1].

####  III.2 The observable

Glauber’s $G^{\left(2\right)}$$\left(\tau\right)$ coincidence signal is
formally given by the expectation value of this observable, evolved using the
total density matrix of the field and matter in the interaction picture

$\displaystyle{\cal C}\left(t,\tau\right)$ $\displaystyle=\left\langle{\cal
T}\hat{{\cal
O}}^{\prime}\left(t,\tau\right)e^{-\frac{i}{\hbar}\underset{-\infty}{\overset{t^{*}}{\int}}dsH{}_{\text{int},-}\left(s\right)}\right\rangle,$
(5)

here ${\cal T}$ is the time ordering superoperator that maintains the
bookkeeping of the interaction events, e.g., ${\cal
T}A\left(t_{1}\right)B\left(t_{2}\right)=\theta\left(t_{1}-t_{2}\right)A\left(t_{1}\right)B\left(t_{2}\right)+\theta\left(t_{2}-t_{1}\right)B\left(t_{2}\right)A\left(t_{1}\right)$,
the Heaviside step-function is defined by $\theta\left(t\right)=1,$ $\forall
t\geq 0$ and $\theta\left(t\right)=0$, $\forall t<0$. Note that the
interaction Hamiltonian superoperator is represented by the field modes prior
to the transformation. HOMI introduces the time delay $T$ as an additional
control parameter to the observable superoperator
$O\left(t,\tau\right)\rightarrow O\left(t,\tau,T\right)$. For the HOM
detection setup, we should transform the observable in Eq. 3 according to the
HOM transformation (Eq. 4), resulting in 16 detection pathways. Only the four
– in which one photon of each mode are detected – contribute to our signal
(see also [9]), reducing Eq. 3 to

$\displaystyle\hat{{\cal O}}^{\prime}\left(t,\tau,T\right)$
$\displaystyle=O_{I}+O_{II}+O_{III}+O_{IV}$ (6)

where the detection pathways are given by

$\displaystyle O_{I}$
$\displaystyle=E_{a,R}^{\dagger}\left(\boldsymbol{r}_{a},t\right)E_{b,R}^{\dagger}\left(\boldsymbol{r}_{b},t+\tau\right)E_{b,L}\left(\boldsymbol{r}_{b},t+\tau\right)E_{a,L}\left(\boldsymbol{r}_{a},t\right)$
(7a) $\displaystyle O_{II}$
$\displaystyle=E_{b,R}^{\dagger}\left(\boldsymbol{r}_{a},t\right)E_{a,R}^{\dagger}\left(\boldsymbol{r}_{b},t+\tau\right)E_{a,L}\left(\boldsymbol{r}_{b},t+\tau\right)E_{b,L}\left(\boldsymbol{r}_{a},t\right)$
(7b) $\displaystyle O_{III}$
$\displaystyle=-E_{a,R}^{\dagger}\left(\boldsymbol{r}_{a},t\right)E_{b,R}^{\dagger}\left(\boldsymbol{r}_{b},t+\tau\right)$
(7c) $\displaystyle\times
E_{a,L}\left(\boldsymbol{r}_{b},t+T+\tau\right)E_{b,L}\left(\boldsymbol{r}_{a},t-T\right)$
$\displaystyle O_{IV}$
$\displaystyle=-E_{b,R}^{\dagger}\left(\boldsymbol{r}_{a},t-T\right)E_{a,R}^{\dagger}\left(\boldsymbol{r}_{b},t+T+\tau\right)$
(7d) $\displaystyle\times
E_{b,L}\left(\boldsymbol{r}_{b},t+\tau\right)E_{a,L}\left(\boldsymbol{r}_{a},t\right).$

Note that since Eq. 6 is given in the basis of the interaction domain
(different from Eq. 3), none of the quantities are primed in the definition of
$\hat{{\cal O}}$ as well as the field operators $E_{m,X}$
($m\in\left\\{a,b\right\\}$; $X\in\left\\{L,R\right\\}$). We have explicitly
included the HOM delay $T$ variable to the coincidence observable, which is
expressed in $O_{III-IV}$ and a $\left(-\right)$ sign [Eq. 6]. All four
combinations depicted in Fig. 2a (top-right) contribute to the interferometric
coincidence signal. When the BS is removed, the ordinary coincidence detection
setup can also be recovered by only keeping the $O_{I}$ contributions.

###  IV The interaction pathways

We expand the signal in Eq. 5 pertubatively to $4^{\text{th}}$ order in
$H_{\text{int}}$ such that each photon interacts twice with the sample.
Generally, 4 interactions generate 16 left-right Liouville pathways. In
addition, each arrow may point inward/outward resulting in a total of 256
possible pathways. As depicted in Figs. 2a and b, this number is significantly
reduced mainly due to the coincidence detection and the initial Fock state.
Note that Fig. 2b contains half of the contributions and their complex
conjugates should be added. It is possible to further reduce the number of
contributions, by the following considerations. Near resonance, the rotating
wave approximation (RWA) can be invoked, resulting in the simplified
interaction Hamiltonian $H_{\text{int}}=EV^{\dagger}+E^{\dagger}V$. We further
consider the three level model systems depicted in Fig. 3, initially in the
ground state $\rho_{0}=|g\rangle\langle g|$, so that the first interaction can
only be excitation (no de-excitation). This eliminates contributions in which
an emission event occurs after a single photon is detected. The expectation
value $\left\langle\hat{{\cal O}}\right\rangle$ is and expectation value and
thus real (note that the diagrams in Fig. 3 are symmetric with respect to
exchange of L-R and taking the complex conjugate). By convention, we only
include pathways in which the last interaction is taken from the left side
with an outgoing arrow (generated a detected photon). The contributions in
which the last interaction is from the right are related to these by
conjugation and interchanging L-R. The full signal is finally given by
$2\mathfrak{Re}\left\langle\hat{{\cal O}}p_{i}\right\rangle$ where $p_{i}$
denotes all pathways terminated at the left. The top interaction from the left
must point outwards, otherwise this mode is not occupied, hence not detected.
Photon number conservation implies equal number of inward/outward arrows. Only
diagrams in which two photons interact with the sample and two photons are
detected contribute to the signal. Since two-photon population is detected,
diagrams in which there is a single arrow in one of the sides are eliminated.

Figure 3: Interaction pathways for three-level exciton model. The five
interaction pathways with the respective exciton dynamics for the selected
three-level exciton model system. Each energy manifold
$|\left\\{\boldsymbol{g}\right\\}\rangle,|\left\\{\boldsymbol{e}\right\\}\rangle,|\left\\{\boldsymbol{f}\right\\}\rangle$
is composed of several states.

The five surviving pathways are depicted in Fig. 3, all contain two field
modes from each side of the density operator. Their complex conjugates should
be added as well. These processes are labeled $P_{i}$ where
$i\in\left[1,5\right]$ in Fig. 3, correspond to the superoperator correlation
functions denoted ${\cal F}_{i}$ given by

$\displaystyle{\cal F}_{1}\left(\tau_{1},\tau_{2},\tau_{3}\right)$
$\displaystyle=\left\langle V_{L}{\cal
G}\left(\tau_{1}\right)V_{R}^{\dagger}{\cal
G}\left(\tau_{2}\right)V_{L}^{\dagger}{\cal
G}\left(\tau_{3}\right)V_{R}\right\rangle$ (8a) $\displaystyle{\cal
F}_{2}\left(\tau_{1},\tau_{2},\tau_{3}\right)$ $\displaystyle=\left\langle
V_{L}{\cal G}\left(\tau_{1}\right)V_{L}^{\dagger}{\cal
G}\left(\tau_{2}\right)V_{R}^{\dagger}{\cal
G}\left(\tau_{3}\right)V_{R}\right\rangle$ (8b) $\displaystyle{\cal
F}_{3}\left(\tau_{1},\tau_{2},\tau_{3}\right)$ $\displaystyle=\left\langle
V_{L}{\cal G}\left(\tau_{1}\right)V_{R}^{\dagger}{\cal
G}\left(\tau_{2}\right)V_{R}{\cal
G}\left(\tau_{3}\right)V_{L}^{\dagger}\right\rangle$ (8c) $\displaystyle{\cal
F}_{4}\left(\tau_{1},\tau_{2},\tau_{3}\right)$ $\displaystyle=\left\langle
V_{L}{\cal G}\left(\tau_{1}\right)V_{L}^{\dagger}{\cal
G}\left(\tau_{2}\right)V_{L}{\cal
G}\left(\tau_{3}\right)V_{L}^{\dagger}\right\rangle$ (8d) $\displaystyle{\cal
F}_{5}\left(\tau_{1},\tau_{2},\tau_{3}\right)$ $\displaystyle=\left\langle
V_{L}{\cal G}\left(\tau_{1}\right)V_{L}{\cal
G}\left(\tau_{2}\right)V_{L}^{\dagger}{\cal
G}\left(\tau_{3}\right)V_{L}^{\dagger}\right\rangle,$ (8e)

(plus their complex conjugates). Here
$\left\langle\cdots\right\rangle\equiv\text{tr}\left\\{\cdots\rho_{\mu}\left(-\infty\right)\right\\}$
where $\rho_{\mu}\left(-\infty\right)$ is the initial state of the matter, and
Liouville-space Green’s function is given by ${\cal
G}\left(t\right)=-\frac{i}{\hbar}\theta\left(t\right)e^{-i{\cal L}_{\mu}t-\eta
t}$.

The matter correlation functions ${\cal F}_{j}(\tau_{1},\tau_{2},\tau_{3})$,
with $j=1,\ldots,5$, presented in Eq. 8a and illustrated in Fig. 3, provide a
useful microscopic insight into the capabilities of the entangled-photon
spectroscopy to retrieve detailed information on ultrafast photoinduced
dynamics of various chemical systems, inter alia molecular aggregates, whose
dynamics is determined by the electronic interactions induced Frenkel exciton
scattering and exciton-phonon interactions. Since the actual measured
correlated signals are represented by convolutions of ${\cal F}_{j}$ with the
doorway and window functions, the latter containing relevant information on
the entangled photon sources, as well as interferometric supplement, and
having nothing to do with the dynamics of the system under study, one should
in first place understand what kind of a matter dynamical information is
contained in the aforementioned four-point correlators.

Since all matter dynamical information, available via coherent four-wave
mixing spectroscopy is fully contained in the third-order nonlinear response
function

$\displaystyle{\cal F}_{\rm coh}=\left\langle\hat{V}_{+}{\cal
G}(\tau_{1})\hat{V}_{-}{\cal G}(\tau_{2})\hat{V}_{-}{\cal
G}(\tau_{3})\hat{V}_{-}\right\rangle,$ (9)

with $\hat{V}=V+V^{\dagger}$, the differences in the spectroscopic
information, provided by entangled-photon versus four-wave mixing
spectroscopies, originates just from the different Liouville-space structure
of the four-point matter correlators ${\cal F}$. Although the latter is
apparently very different for ${\cal F}$, presented in Eq. 8a versus Eq. 9,
the issue allows for a clear, simple, and conceptually explicit analysis, as
follows.

Indeed, as shown, using Liouville space Green functions techniques [36, 37]
and further by means of the Nonlinear Exciton Equations (NEE) [38], in the
case of weak to moderate exciton-phonon coupling, when the dynamics is
dominated by excitonic effects, and effects of polaron formation/self-trapping
are not substantial, there are three phenomena that contribute to optical
response: (i) exciton-exciton scattering, described in terms of the exciton-
exciton scattering matrix $\bar{\Gamma}(\omega)$, (ii) exciton-photon coupling
mediated transport, described by the exciton transport correlation function
$G^{(2)}(\omega)$, and combined effects, expressed in terms of a convolutions
of the two above; of course, the expressions for the response contain the one-
exciton Green function $G(\omega)$ that contains of the exciton-phonon
coupling induced exciton dephasing.

The Green function approach can be extended to analyze the Liouville space
correlators ${\cal F}_{j}$ in Eq. (8); this analysis will be addressed in
detail in a separate publication, here we just present its main outcome, to
relate it to the Feynman diagrams (Fig. 3). Despite a very different
Liouville-space structure of the correlators in the entangled photon versus
four-wave mixing case, due to the specific features of the Frenkel exciton
model with moderate exciton-photon coupling, the ingredients that enter the
final expressions, namely $\bar{\Gamma}(\omega)$, $G^{(2)}(\omega)$, and
$G(\omega)$ stay the same, just the final expressions get modified.
Translating the results, presented in Eq. 8a from the Sum-over-States (SOS) to
the exciton-scattering language, we can combine the Feynman diagrams (4) and
(5) in Fig. 3 to obtain type (i) effects, i.e., pure exciton-excitons
scattering; it is well known that combining these two diagrams we take care of
the so-called cancellation of the $N^{2}$ terms problem, which in the exciton-
scattering approach happens automatically. Combining the diagrams (1), (2),
and (3) we obtain a type (ii) contribution that reflects exciton transport
effects, since in diagrams (1) and (3) the system is in the
population/exciton-exciton coherence state during the $\tau_{2}$ time segment.
In diagram 2 the system is in the $|g\rangle\langle g^{\prime}|$, which means
in the ground electronic state with a different phonon structures, it plays a
proper cancellation role for exciton transport, in the way how the diagram (4)
operates for exciton scattering. Note that both diagrams (2) and (4) also
describe slight modifiication/renormalization of the excitons scattering
matrix due to exciton-phonon coupling.

Very importantly type (iii) effects (combined exciton-exciton scattering and
exciton transport) in the four-wave mixing, on the SOS language originate from
the diagrams, when the system is in the population/exciton-exciton coherence
and two-exciton ground state coherence ($|f\rangle\langle g|$ or
$|g\rangle\langle f|$) during the time periods $\tau_{2}$ and $\tau_{3}$,
respectively. Since such diagrams never appear in the entangled-photon
spectroscopy case, type (iii) effects do not contribute to the signals in the
latter case.

Summarizing, unlike the coherent four-wave mixing spectroscopy, the entangled-
photon spectroscopy studies only effects of exciton-exciton scattering and
exciton transport, with no combined contributions, thus provided a better
separation of dynamical phenomena that contribute to spectroscopic data.

####  IV.1 The signal

In Fig. 1a we consider the HOMI detection setup. The preparation alters the
primary source, resulting in the JSA given in Eq. 2. The HOMI detection
transforms the $G^{\left(2\right)}\left(\tau\right)$ signal in Eq. 3 into Eq.
6.

The total signal involves summing over all the product combinations of
detection pathways in Fig. 2 a with all the interaction pathways depicted in
Fig. 2b (see appendix B for the expressions of all combinations for a general
preparation process). The signal contains many terms corresponding to all
combinations of Liouville pathways $i\in\left[1,5\right]$, with all detection
pathways defined in Eq. 6 $O_{\nu}$, where $\nu=I,...,IV$. One way to think
about it is that each process $P_{i}$ is obtained by a coherent superposition
of all HOM detection pathways ${\cal S}_{\nu,i}$, and the signal is given by
their superposition. For an illustrative example of the derivation of an
contribution, see appendix B. The coincidence signal in Eq. 10 is finally
given by

$\displaystyle{\cal C}\left(\tau,T,s\right)$
$\displaystyle=2\mathfrak{Re}\left\\{\sum_{\nu,i}{\cal S}_{\nu,i}\right\\},$
(10)

where ${\cal
S}_{\nu,i}\equiv\text{tr}\left\\{O_{\nu}P_{i}\rho\left(-\infty\right)\right\\}$
and the detection and interaction pathways are labeled in Fig. 2. The
detection pathways are given in Eq. 6 (see appendices for detailed
expressions). Note that the density matrix is given by a product of the matter
and field respectively
$\rho\left(-\infty\right)\equiv\rho_{\mu}\otimes\rho_{\varphi}$. The field is
traced with respect to
$\rho_{\varphi}\equiv|\Psi_{\theta}\rangle\langle\Psi_{\theta}|$ using
$\Phi_{\theta}\left(t_{1},t_{2}\right)=\frac{1}{\sqrt{2}}\left[\phi\left(t_{1},t_{2}\right)+e^{-i\theta}\phi\left(t_{2},t_{1}\right)\right]$
with respect to Eq. 1 due to the Michelson interferometer.

### The short entanglement-time limit

We no invoke an approximation which greatly simplifies this signal. Consider a
symmetric joint spectral amplitude, obtained by either using $\theta=0$ in Eq.
2, or via a narrowband pump. In either case, the entanglement time $T_{e}$
represents the time window in which both photons arrive [30]. We also consider
the characteristic timescale for the matter dynamics to be bound from above by
$\tau_{\text{R}}$. We focus on the regime $\tau_{R}\gg T_{e}$ such that both
photons arrive simultaneously. The relative delay $s$ introduced in Fig. 1a
between the pair, now sets the time interval in which all interactions occur.
In this limit, the amplitude is approximated by a narrow distribution

$\Phi\left(t_{1},t_{2}\right)\rightarrow\delta\left(t_{1}-t_{2}-s\right).$
(11)

Consequently, processes with vanishing time intervals between interactions do
not contribute to this order $\left(\tau_{1},\tau_{2},\tau_{3}\right)\neq 0$.
Note that $\Phi\left(t_{1},t_{2}\right)$ is symmetric to $t_{1}\leftrightarrow
t_{2}$ exchange under this approximation which is consistent with $\theta=0$.
Here the Michelson interferometer is used for rectification of the exchange
phase $\theta$ when an ultrafast pump is used. Alternatively, a narrowband
pump can be used in which the exchange phase correction is no longer
essential. Since $\left(T,\tau,s\right)$ are measured on a finite grid, we
define the discrete time delta distribution $\delta_{t_{1},t_{2}}$ which
attains the value $1$ when $t_{1}=t_{2}$ within our setup. Plugging Eq. 11, in
the signal $\left(\text{Eq. }\text{\ref{eq:
signal_sup}},\text{appendix}\right)$ for $\tau\geq-T$ and $s>0$ we obtain

$\displaystyle{\cal C}\left(\tau>-T,T,s\right)=$ $\displaystyle
2\text{$\mathfrak{Re}$$\left\\{-\theta\left(\tau+T\right){\cal
F}_{1}\left(T,\tau+T-s,2s-T\right)|_{E}\right.$}$
$\displaystyle-\theta\left(\tau+T\right){\cal
F}_{2}\left(2T+\tau-s,s-T-\tau,\tau+s\right)|_{E}$
$\displaystyle-\theta\left(\tau+T\right)\left[{\cal
F}_{3}\left(T,\tau+s,T-2s\right)+{\cal
F}_{3}\left(T+\tau,s-2T-\tau,T+\tau\right)\right]|_{E}$
$\displaystyle+2\delta_{\tau,s}\int d\tau_{3}{\cal
F}_{5}\left(\left|\tau\right|,\tau_{3},\left|\tau\right|\right)|_{D}\left.-2\delta_{2T+\tau,s}\int
d\tau_{3}{\cal
F}_{5}\left(\left|\tau\right|,\tau_{3},2T+\tau\right)|_{E}\right\\}$ (12)

where the coincidence contribution ${\cal C}_{4}$ does not appear in this
limit. We have introduced the notation $|_{E/D}$ corresponding to the direct
$\left(O_{I,II}\right)$ and exchange $\left(O_{III,IV}\right)$ paths of the
HOM interferometer. From this we see that only process $i=5$ in Fig. 3
contributes to both the direct paths $O_{I/II}$, the rest are limited to the
exchanges $O_{III/IV}$. In this limit we already appreciate the degree of
control offered by the interferometric setups, offering novel temporal
inference tool-box. For example, for $s=0$, only exchange path processes may
contribute since $s$ set up the scale in which all the interactions with the
sample occur. HOM exchange paths are not restricted by this due to the
ambiguity in the arrival times. We thus single out ${\cal F}_{3}$
contribution, ${\cal
C}\left(\tau>-T,T,s=0\right)\propto-\text{$\theta\left(\tau+T\right)$}\mathfrak{Re}\left\\{{\cal
F}_{3}\left(T,\tau,T\right)\right\\}$. Also, for $\tau=0$ we obtain

$\displaystyle{\cal C}\left(\tau=0,s<T<2s,s\right)$
$\displaystyle\propto\theta\left(T\right)\text{$\mathfrak{Re}$$\left\\{{\cal
F}_{1}\left(T,T-s,2s-T\right)\right\\}$ },$ (13a) $\displaystyle{\cal
C}\left(\tau=0,\frac{s}{2}<T<s,s\right)$
$\displaystyle\propto\theta\left(T\right)\text{$\mathfrak{Re}$$\left\\{{\cal
F}_{2}\left(2T-s,s-T,s\right)\right\\}$ },$ (13b) $\displaystyle{\cal
C}\left(\tau=0,2s<T,s\right)$
$\displaystyle\propto\theta\left(T\right)\text{$\mathfrak{Re}$$\left\\{{\cal
F}_{3}\left(T,s,T-2s\right)\right\\}$ }.$ (13c)

Eqs. 13a-c demonstrate that by following the multidimensional data, it is
possible to isolate certain contributions in the time domain.

###  V Discussion

The setup above presents several types of control variable over the signal.
These can be categorized in three groups: (a) classical pump, (b) preparation
and (c) detection parameters. The pump related parameters includes the central
frequency of the pump and its spectral width
$\left(\omega_{p},\sigma_{p}\right)$. The preparation setup is rich with
parameters including the central frequencies of the daughter photons and their
respective time $\left(\bar{\omega}_{a},\bar{\omega}_{b},T_{a},T_{b}\right)$,
using the phase matching conditions and dispersion properties of the nonlinear
crystal at these frequencies (see Sec. II). These parameters were not scanned
here and offer a rich playground for future studies. The detection parameters
include the number of detected photons (here two) and the HOM Delay $T$.

Interferometric spectroscopy with quantum light has several merits. Due to the
specified number of interacting and detected photons, certain pathways which
contribute to classical signals, are eliminated. This feature has strictly
quantum origin since we are using Fock states. Due to the application and
detection of fixed number of photons, the signal records only processes that
lie within the two-photon subspace. This greatly reduce the number of
Liouville-space pathways. One can define an entropic measure from the
Liouville pathways probability functional which depends on the preparation and
detection details. Ultimately, it is possible to identify the activated
pathways by shaping this probability with the available control parameters.
This represents the quantum information gain obtained by the protocol. One way
to see that is by defining a pathway related entropy ${\cal
S}_{0}=-\sum_{i=1}^{L}P_{i}\log P_{i}$ where $L$ is the number of Liouville
pathways with a given probe and $P_{i}$ is the probability of the
$i^{\text{th}}$ pathway. Then compare it to the entropy of the interferometric
setup ${\cal S}=-\sum_{i=1}^{L}Q_{i}\log Q_{i}$ where $Q_{i}$ is the overall
probability of the $i^{\text{th}}$ pathway in the manipulated scheme.
Ultimately, one can quantify the quantum inference due to the detection
process alone using an identical probe and different detection by calculating
the Kullback-Liebler divergence ${\cal
D}\left[P||Q\right]=\sum_{i=1}^{L}p_{i}\log\frac{P_{i}}{Q_{i}}$ [39].

Each of the delays in this setup affects the signal differently, which allows
to control the pathways in the time domain. For example, by taking
$\Phi\left(t_{1},t_{2}\right)\rightarrow\delta\left(t_{1}-t_{2}-s\right),$ the
entire process duration is set by $s$ for the direct detection pathways. This
stems from the fact that one photon is absorbed and emitted while its
entangled partner goes straight to the detector in some pathways. One observes
a superposition of processes in which both photon interact with the sample
with a conjugate process in which both had not. Therefore, the two-photon
coherence time sets the characteristic interaction duration for the direct
pathways. The exchange detection pathways have more freedom due to the
superposition in time domain introduced by the HOMI. When the material system
under study possesses several characteristic time-scales $\tau_{R}$, they can
be studied separately by adjusting $s$.

Coincidence-detection further involves unique scaling relations between the
applied intensity $I_{p}$, the light-sample coupling and the detected signal.
This allows to avoid damaging the sample by using weak quantum fields [2].

###### Acknowledgements.

The support of the National Science Foundation (NSF) Grant CHE-1953045 is
gratefully acknowledged. The research was supported by the U.S. Department of
Energy (DOE), Office of Science, Basic Energy Sciences, under Award Number DE-
SC0022134. S.M. and V.C. were supported by the DOE award.

### References

  * [1] Shahaf Asban, Konstantin E. Dorfman, and Shaul Mukamel.  Interferometric spectroscopy with quantum light: Revealing out-of-time-ordering correlators.  The Journal of Chemical Physics, 154(21):210901, 2021.
  * [2] Shahaf Asban and Shaul Mukamel.  Distinguishability and which pathway information in multidimensional interferometric spectroscopy with a single entangled photon-pair.  Science Advances, 7(39):eabj4566, 2021.
  * [3] Konstantin E. Dorfman, Shahaf Asban, Bing Gu, and Shaul Mukamel.  Hong-ou-mandel interferometry and spectroscopy using entangled photons.  Communications Physics, 4(1):49, Mar 2021.
  * [4] Audrey Eshun, Bing Gu, Oleg Varnavski, Shahaf Asban, Konstantin E. Dorfman, Shaul Mukamel, and Theodore Goodson.  Investigations of molecular optical properties using quantum light and hong–ou–mandel interferometry.  Journal of the American Chemical Society, 143(24):9070–9081, Jun 2021.
  * [5] Szilard Szoke, Hanzhe Liu, Bryce P. Hickam, Manni He, and Scott K. Cushing.  Entangled light–matter interactions and spectroscopy.  Journal of Materials Chemistry C, 8(31):10732–10741, 2020.
  * [6] Shaul Mukamel, Matthias Freyberger, Wolfgang Schleich, Marco Bellini, Alessandro Zavatta, Gerd Leuchs, Christine Silberhorn, Robert W Boyd, Luis Lorenzo Sánchez-Soto, André Stefanov, et al.  Roadmap on quantum light spectroscopy.  Journal of Physics B: Atomic, Molecular and Optical Physics, 53(7):072002, 2020.
  * [7] Shahaf Asban, Konstantin E. Dorfman, and Shaul Mukamel.  Quantum phase-sensitive diffraction and imaging using entangled photons.  Proceedings of the National Academy of Sciences, 116(24):11673–11678, 2019.
  * [8] JG Rarity, PR Tapster, E Jakeman, T Larchuk, RA Campos, MC Teich, and BEA Saleh.  Two-photon interference in a mach-zehnder interferometer.  Physical review letters, 65(11):1348, 1990.
  * [9] C. K. Hong, Z. Y. Ou, and L. Mandel.  Measurement of subpicosecond time intervals between two photons by interference.  Phys. Rev. Lett., 59:2044–2046, Nov 1987.
  * [10] MG Raymer, Andrew H Marcus, Julia R Widom, and Dashiell LP Vitullo.  Entangled photon-pair two-dimensional fluorescence spectroscopy (epp-2dfs).  The Journal of Physical Chemistry B, 117(49):15559–15575, 2013\.
  * [11] AA Kalachev, DA Kalashnikov, AA Kalinkin, TG Mitrofanova, AV Shkalikov, and VV Samartsev.  Biphoton spectroscopy in a strongly nondegenerate regime of spdc.  Laser Physics Letters, 5(8):600–602, 2008.
  * [12] Jonathan Lavoie, Tiemo Landes, Amr Tamimi, Brian J. Smith, Andrew H. Marcus, and Michael G. Raymer.  Phase-modulated interferometry, spectroscopy, and refractometry using entangled photon pairs.  Advanced Quantum Technologies, 3(11):1900114, 2020.
  * [13] W. P. Grice and I. A. Walmsley.  Spectral information and distinguishability in type-ii down-conversion with a broadband pump.  Phys. Rev. A, 56:1627–1634, Aug 1997.
  * [14] Carlton M. Caves.  Quantum-mechanical noise in an interferometer.  Phys. Rev. D, 23:1693–1708, Apr 1981.
  * [15] Carlton M. Caves and Bonny L. Schumaker.  New formalism for two-photon quantum optics. i. quadrature phases and squeezed states.  Phys. Rev. A, 31:3068–3092, May 1985.
  * [16] M. Tse, Haocun Yu, N. Kijbunchoo, A. Fernandez-Galiana, P. Dupej, L. Barsotti, C. D. Blair, D. D. Brown, S. E. Dwyer, A. Effler, M. Evans, P. Fritschel, V. V. Frolov, A. C. Green, G. L. Mansell, F. Matichard, N. Mavalvala, D. E. McClelland, L. McCuller, T. McRae, J. Miller, A. Mullavey, E. Oelker, I. Y. Phinney, D. Sigg, B. J. J. Slagmolen, T. Vo, R. L. Ward, C. Whittle, R. Abbott, C. Adams, R. X. Adhikari, A. Ananyeva, S. Appert, K. Arai, J. S. Areeda, Y. Asali, S. M. Aston, C. Austin, A. M. Baer, M. Ball, S. W. Ballmer, S. Banagiri, D. Barker, J. Bartlett, B. K. Berger, J. Betzwieser, D. Bhattacharjee, G. Billingsley, S. Biscans, R. M. Blair, N. Bode, P. Booker, R. Bork, A. Bramley, A. F. Brooks, A. Buikema, C. Cahillane, K. C. Cannon, X. Chen, A. A. Ciobanu, F. Clara, S. J. Cooper, K. R. Corley, S. T. Countryman, P. B. Covas, D. C. Coyne, L. E. H. Datrier, D. Davis, C. Di Fronzo, J. C. Driggers, T. Etzel, T. M. Evans, J. Feicht, P. Fulda, M. Fyffe, J. A. Giaime, K. D. Giardina, P. Godwin, E. Goetz, S. Gras, C. Gray, R. Gray, Anchal Gupta, E. K. Gustafson, R. Gustafson, J. Hanks, J. Hanson, T. Hardwick, R. K. Hasskew, M. C. Heintze, A. F. Helmling-Cornell, N. A. Holland, J. D. Jones, S. Kandhasamy, S. Karki, M. Kasprzack, K. Kawabe, P. J. King, J. S. Kissel, Rahul Kumar, M. Landry, B. B. Lane, B. Lantz, M. Laxen, Y. K. Lecoeuche, J. Leviton, J. Liu, M. Lormand, A. P. Lundgren, R. Macas, M. MacInnis, D. M. Macleod, S. Márka, Z. Márka, D. V. Martynov, K. Mason, T. J. Massinger, R. McCarthy, S. McCormick, J. McIver, G. Mendell, K. Merfeld, E. L. Merilh, F. Meylahn, T. Mistry, R. Mittleman, G. Moreno, C. M. Mow-Lowry, S. Mozzon, T. J. N. Nelson, P. Nguyen, L. K. Nuttall, J. Oberling, R. J. Oram, B. O’Reilly, C. Osthelder, D. J. Ottaway, H. Overmier, J. R. Palamos, W. Parker, E. Payne, A. Pele, C. J. Perez, M. Pirello, H. Radkins, K. E. Ramirez, J. W. Richardson, K. Riles, N. A. Robertson, J. G. Rollins, C. L. Romel, J. H. Romie, M. P. Ross, K. Ryan, T. Sadecki, E. J. Sanchez, L. E. Sanchez, T. R. Saravanan, R. L. Savage, D. Schaetzl, R. Schnabel, R. M. S. Schofield, E. Schwartz, D. Sellers, T. J. Shaffer, J. R. Smith, S. Soni, B. Sorazu, A. P. Spencer, K. A. Strain, L. Sun, M. J. Szczepańczyk, M. Thomas, P. Thomas, K. A. Thorne, K. Toland, C. I. Torrie, G. Traylor, A. L. Urban, G. Vajente, G. Valdes, D. C. Vander-Hyde, P. J. Veitch, K. Venkateswara, G. Venugopalan, A. D. Viets, C. Vorvick, M. Wade, J. Warner, B. Weaver, R. Weiss, B. Willke, C. C. Wipf, L. Xiao, H. Yamamoto, M. J. Yap, Hang Yu, L. Zhang, M. E. Zucker, and J. Zweizig.  Quantum-enhanced advanced ligo detectors in the era of gravitational-wave astronomy.  Phys. Rev. Lett., 123:231107, Dec 2019.
  * [17] F. Hudelist, Jia Kong, Cunjin Liu, Jietai Jing, Z. Y. Ou, and Weiping Zhang.  Quantum metrology with parametric amplifier-based photon correlation interferometers.  Nature Communications, 5:3049, Jan 2014.
  * [18] Dong Li, Chun-Hua Yuan, Z Y Ou, and Weiping Zhang.  The phase sensitivity of an SU(1,1) interferometer with coherent and squeezed-vacuum light.  New Journal of Physics, 16(7):073020, jul 2014.
  * [19] Brian E. Anderson, Prasoon Gupta, Bonnie L. Schmittberger, Travis Horrom, Carla Hermann-Avigliano, Kevin M. Jones, and Paul D. Lett.  Phase sensing beyond the standard quantum limit with a variation on the su(1,1) interferometer.  Optica, 4(7):752–756, Jul 2017.
  * [20] Mathieu Manceau, Gerd Leuchs, Farid Khalili, and Maria Chekhova.  Detection loss tolerant supersensitive phase measurement with an su(1,1) interferometer.  Phys. Rev. Lett., 119:223604, Nov 2017.
  * [21] Yaakov Shaked, Yoad Michael, Rafi Z. Vered, Leon Bello, Michael Rosenbluh, and Avi Pe’er.  Lifting the bandwidth limit of optical homodyne measurement with broadband parametric amplification.  Nature Communications, 9(1):609, Feb 2018.
  * [22] G. Frascella, E. E. Mikhailov, N. Takanashi, R. V. Zakharov, O. V. Tikhonova, and M. V. Chekhova.  Wide-field su(1,1) interferometer.  Optica, 6(9):1233–1236, Sep 2019.
  * [23] Shaul Mukamel.  Principles of Nonlinear Optical Spectroscopy.  Oxford University Press, 1995.
  * [24] R D Mota, M A Xicoténcatl, and V D Granados.  Jordan–schwinger map, 3d harmonic oscillator constants of motion, and classical and quantum parameters characterizing electromagnetic wave polarization.  Journal of Physics A: Mathematical and General, 37(7):2835–2842, feb 2004.
  * [25] R. D. Mota, D. Ojeda-Guillén, M. Salazar-Ramírez, and V. D. Granados.  Su(1,1) approach to stokes parameters and the theory of light polarization.  J. Opt. Soc. Am. B, 33(8):1696–1701, Aug 2016.
  * [26] Roy J. Glauber.  The quantum theory of optical coherence.  Phys. Rev., 130:2529–2539, Jun 1963.
  * [27] Darius Abramavicius, Benoit Palmieri, Dmitri V. Voronine, František Šanda, and Shaul Mukamel.  Coherent multidimensional optical spectroscopy of excitons in molecular aggregates; quasiparticle versus supermolecule perspectives.  Chemical Reviews, 109(6):2350–2408, Jun 2009.
  * [28] Gregory D. Scholes, Graham R. Fleming, Lin X. Chen, Alán Aspuru-Guzik, Andreas Buchleitner, David F. Coker, Gregory S. Engel, Rienk van Grondelle, Akihito Ishizaki, David M. Jonas, Jeff S. Lundeen, James K. McCusker, Shaul Mukamel, Jennifer P. Ogilvie, Alexandra Olaya-Castro, Mark A. Ratner, Frank C. Spano, K. Birgitta Whaley, and Xiaoyang Zhu.  Using coherence to enhance function in chemical and biophysical systems.  Nature, 543(7647):647–656, Mar 2017.
  * [29] Hohjai Lee, Yuan-Chung Cheng, and Graham R. Fleming.  Coherence dynamics in photosynthesis: Protein protection of excitonic coherence.  Science, 316(5830):1462–1465, 2007.
  * [30] D. Branning, W. P. Grice, R. Erdmann, and I. A. Walmsley.  Engineering the indistinguishability and entanglement of two photons.  Phys. Rev. Lett., 83:955–958, Aug 1999.
  * [31] David Branning, Warren Grice, Reinhard Erdmann, and I. A. Walmsley.  Interferometric technique for engineering indistinguishability and entanglement of photon pairs.  Phys. Rev. A, 62:013814, Jun 2000.
  * [32] C. K. Law, I. A. Walmsley, and J. H. Eberly.  Continuous frequency entanglement: Effective finite hilbert space and entropy control.  Phys. Rev. Lett., 84:5304–5307, Jun 2000.
  * [33] Bernard Yurke, Samuel L. McCall, and John R. Klauder.  Su(2) and su(1,1) interferometers.  Phys. Rev. A, 33:4033–4054, Jun 1986.
  * [34] F. Rohrlich J. M. Jauch.  The Theory of Photons and Electrons.  Springer, Berlin, Heidelberg, 1976.
  * [35] R D Mota, M A XicotÃ©ncatl, and V D Granados.  Two-dimensional isotropic harmonic oscillator approach to classical and quantum stokes parameters.  Canadian Journal of Physics, 82(10):767–773, 2004.
  * [36] Vladimir Chernyak, Ningjun Wang, and Shaul Mukamel.  Four-wave mixing and luminescence of confined excitons in molecular aggregates and nanostructures. many-body green function approach.  Physics Reports, 263(4):213–309, 1995.
  * [37] T. Meier, V. Chernyak, and S. Mukamel.  Multiple exciton coherence sizes in photosynthetic antenna complexes viewed by pump-probe spectroscopy.  The Journal of Physical Chemistry B, 101(37):7332–7342, Sep 1997\.
  * [38] Vladimir Chernyak, Wei Min Zhang, and Shaul Mukamel.  Multidimensional femtosecond spectroscopies of molecular aggregates and semiconductor nanostructures: The nonlinear exciton equations.  The Journal of Chemical Physics, 109(21):9587–9601, 1998.
  * [39] S. Kullback and R. A. Leibler.  On information and sufficiency.  The Annals of Mathematical Statistics, 22(1):79–86, 1951.
  * [40] Shaul Mukamel.  Partially-time-ordered schwinger-keldysh loop expansion of coherent nonlinear optical susceptibilities.  Phys. Rev. A, 77:023801, Feb 2008.

##

###  Appendix A Hilbert-space approach to photon counting

In this section we propose an alternative derivation in which the entire
calculation is computed in Hilbert space. This representation has the
advantage of offering more compact expressions reflected in less diagrams.
However, when external degrees of freedom are included to account for
inaccessible processes as a result of possible coupling to the environment,
this is no longer possible.

Our observable is expressed via interaction of an field mode with the
detector, changing its polarization. This can be described using perturbative
expansion of the interaction Hamiltonian with the detector degrees of freedom
such that each interaction event contributes a single interaction to the
wavefunction that describes the light, sample and detector. An $N$ photon
measurement operator corresponds to

${\cal M}^{\prime}={\cal T}\int\cdots\int dt_{1}\cdot\dots
dt_{N}d\boldsymbol{r}_{1}\cdot\dots
d\boldsymbol{r}_{N}H_{\text{int}}\left(t_{1}\right)\cdot\dots
H_{\text{int}}\left(t_{N}\right)$ (14)

where each active pixel correspond to a single interaction Hamiltonian
$H_{\text{int}}\left(t\right)=\boldsymbol{E}\cdot\boldsymbol{V}$. When the
detector’s dipole response is taken to be small and fast, we approximate it as
delta distribution in space-time and obtain the known Glauber detection
scheme. We consider an ordered measurement scheme without loss of generality
for the two-photon coincidence scheme. We also assume that the detection plane
is far from the sample such that the time ordering operator does not mix the
two photon detection with the light-matter coupling. The wavefunction of the
sample and the detectors at positions
$\left\\{\boldsymbol{r}_{i},t_{i}\right\\}$ is separable and therefore takes
the form

$|\Theta\left(t\right);\left\\{\boldsymbol{r}_{i},t_{i}\right\\}\rangle_{\mu\varphi
M}=|\Psi\left(t\right)\rangle_{\mu\varphi}|g,g\rangle_{M},$

here he subscripts $M,\varphi,\mu$ describes the detectors, field and sample
respectively. After the interaction with the detectors

$\displaystyle|\Theta\left(t\right);\left\\{\boldsymbol{r}_{i},t_{i}\right\\}\rangle_{\mu\varphi
M}$ $\displaystyle\underset{\text{2
detectors}}{=}E\left(\boldsymbol{r}_{a},t_{a}\right)E\left(\boldsymbol{r}_{b},t_{b}\right)|\Psi\left(t\right)\rangle_{\mu\varphi}|e_{1},e_{2}\rangle_{M}$
$\displaystyle={\cal
M}\left(\boldsymbol{r}_{a},\boldsymbol{r}_{b};t_{a},t_{b}\right)|\Psi\left(t\right)\rangle_{\mu\varphi}|e_{1},e_{2}\rangle_{M},$
(15)

Developing the light-sample wavefunction pertubatively we obtain

$\displaystyle|\Psi\left(t\right)\rangle_{\mu\varphi}$
$\displaystyle=\sum_{k=0}|\Psi^{\left(k\right)}\left(t\right)\rangle_{\mu\varphi}$
(16) $\displaystyle|\Psi^{\left(k\right)}\left(t\right)\rangle_{\mu\varphi}$
$\displaystyle=\left(-\frac{i}{\hbar}\right)^{k}{\cal T}_{+}\int\cdots\int
dt_{1}\cdot\dots dt_{k}d\boldsymbol{r}_{1}\cdot\dots
d\boldsymbol{r}_{k}H_{\text{int}}\left(t_{1}\right)\cdot\dots
H_{\text{int}}\left(t_{k}\right)|\Psi\left(-\infty\right)\rangle_{\mu\varphi},$
(17)

where ${\cal T}_{+}$ denotes the time ordering operator forward in time ${\cal
T}_{+}A\left(t_{1}\right)B\left(t_{2}\right)=\theta\left(t_{1}-t_{2}\right)A\left(t_{1}\right)B\left(t_{2}\right)+\theta\left(t_{2}-t_{1}\right)B\left(t_{2}\right)A\left(t_{1}\right)$.
Note that we introduced the $\left(+\right)$ subscript to the time ordering
since the Hermitian conjugate evolves formally backwards in time using ${\cal
T}_{-}$ to the left of the observable. We calculate the probability of this
detection setup by taking the modulus square of this amplitude, resulting in

$\displaystyle{\cal P}\left(\left\\{\boldsymbol{r}_{i},t_{i}\right\\}\right)$
$\displaystyle=\langle\Theta\left(t\right);\left\\{\boldsymbol{r}_{i},t_{i}\right\\}|\Theta\left(t\right);\left\\{\boldsymbol{r}_{i},t_{i}\right\\}\rangle_{\mu\varphi
M},$ (18) $\displaystyle=\langle\Psi\left(t\right){\cal
M}^{\dagger}\left(\left\\{\boldsymbol{r}_{i},t_{i}\right\\}\right){\cal
T}_{-}|{\cal T}_{+}{\cal
M}\left(\left\\{\boldsymbol{r}_{i},t_{i}\right\\}\right)\Psi\left(t\right)\rangle_{\mu\varphi}$
(19)

explicitly

${\cal
P}\left(\boldsymbol{r}_{a},\boldsymbol{r}_{b};t_{a},t_{b}\right)=\sum_{k,l=0}^{\infty}\langle\Psi^{\left(l\right)}\left(t\right)E^{\dagger}\left(\boldsymbol{r}_{b},t_{b}\right)E^{\dagger}\left(\boldsymbol{r}_{a},t_{a}\right){\cal
T}_{-}|{\cal
T}_{+}E\left(\boldsymbol{r}_{a},t_{a}\right)E\left(\boldsymbol{r}_{b},t_{b}\right)\Psi^{\left(k\right)}\left(t\right)\rangle_{\mu\varphi}.$
(20)

This equation is exact in the light-matter interaction and pertubative in the
interaction with the detector. Note that until this stage, all field-source
contractions are permitted such that emission of a photon can occur after the
detection of another. From this point, we assume that the detectors are placed
far from the interaction area such that one can assume the time ordering
applies for the light-matter interaction solely and the detection events are
ordered by definition $\left(t_{a}>t_{b}\right)$

${\cal
P}\left(\boldsymbol{r}_{a},\boldsymbol{r}_{b};t_{a},t_{b}\right)=\sum_{k,l=0}\langle\Psi^{\left(l\right)}\left(t\right)|E^{\dagger}\left(\boldsymbol{r}_{b},t_{b}\right)E^{\dagger}\left(\boldsymbol{r}_{a},t_{a}\right)E\left(\boldsymbol{r}_{a},t_{a}\right)E\left(\boldsymbol{r}_{b},t_{b}\right)|\Psi^{\left(k\right)}\left(t\right)\rangle_{\mu\varphi}.$
(21)

Figure 4: Diagrammatic representation of typical term in Eq. 21. The left
(right) branch of the diagram represents the forward (backward) evolution of
the joint light-matter wavefunction. The number of interactions in each branch
correspond to the respective order of perturbation of either the bra and the
ket in the interaction picture. The interaction and detection interval in this
calculation do not overlap although in principle the may.

Each term in Eq.21 can be represented using a fully time ordered loop diagram
as depicted in Fig. 4. this diagram represents the forward (in time) evolution
of the ket (left) pertubatively in the interaction picture to the
$k^{\text{th}}$ order, and the bra backwards in time to the $l^{\text{th}}$
order along the time contour $\tau_{C}$. Between interaction events, the
sample and the electromagnetic field are evolved using their free Hamiltonian.
The lowest order in which two photons interact with the sample and detected is
the $4^{\text{th}}$. Contributions with odd number of photons from the left or
right at the detection are naturally eliminated. This corresponds to

${\cal
P}\left(\boldsymbol{r}_{a},\boldsymbol{r}_{b};t_{a},t_{b}\right)\approx\sum_{k,l=0}^{4}\langle\Psi^{\left(l\right)}\left(t\right)|E_{-}^{\dagger}\left(\boldsymbol{r}_{b},t_{b}\right)E_{-}^{\dagger}\left(\boldsymbol{r}_{a},t_{a}\right)E_{+}\left(\boldsymbol{r}_{a},t_{a}\right)E_{+}\left(\boldsymbol{r}_{b},t_{b}\right)|\Psi^{\left(k\right)}\left(t\right)\rangle_{\mu\varphi},$
(22)

where the $\pm$ subscripts highlight the operation direction of the field with
respect to the time contour $\tau_{C}$, for positive and negative time
direction. Equation 22 gives rise to two kinds of contributions: (1) four
interactions in one side and non at the other and (2) two interactions from
each side. The Hilbert space description – while equivalent to the alternative
Liouville space – results in partial time ordering. The time ordering is
maintained along the contour $\tau_{C}$, thus, along the left and right
branches of the diagram in Fig. 4individually. Alternatively, if one is
interested in evolving the density matrix (in Liouville-space), the relative
time-ordering of left and right branches is also important, resulting in
absolute time ordering. This difference become important in the interpretation
of pertubative treatments of light-matter coupling and essential when one
considers coupling to reservoirs. As a result, in the wavefunction approach
(Hilbert space) the relative coherence during the evolution is not expressed,
only the the final phase accumulated along the entire evolution of the bra and
ket separately such that over-all coherence is accounted for as shown in Fig.
5a. Alternatively, in the Liouville-space approach the change in coherence is
instantaneously monitored in the calculation process as demonstrated in
Fig.5b. For time dependent perturbation theory that includes terms which break
time reversal (bath), using the Liouville space approach is inevitable and
thus invoked in this paper.

Figure 5: Partial Vs. full time ordering. a) a typical $4^{\text{th}}$ order
term in the wavefunction pertubative treatment. b.) Similar terms in Liouville
space expansion. Breaking the loop contour (forward and backward) into a
single time evolution for both the bra and the ket.

###  Appendix B The signal

The signal in Eq. 10 is composed of all processes $P_{i}$ evaluated with all
observables $O_{j}$. One way to think about it is that each process $P_{i}$ is
obtained by coherent superposition of all HOM detection pathways ${\cal
C}_{i}\left(t,\tau,T\right)$, and the signal is given by their superposition.
The coincidence signal in Eq. 10 can be written accordingly

${\cal C}\left(t,\tau,T\right)=2\mathfrak{Re}\left\\{\sum_{i=1}^{5}{\cal
S}_{\nu,i}\right\\}$ (23)

and solved for each detection-interaction pathway combination below
separately.

####  B.1 Example of one process-observable combination

We now illustrate how to combine the preparation and observation boxes for a
single term from the total signal. We chose ${\cal S}_{I,1}$ as shown in Fig.
6. This contribution introduces four combinations of field modes corresponding
to coupling with $a$ and $b$ modes, $\left(aa,bb,ab,ba\right)$ . We consider
the realization in which $a$ is coupled from the left and $b$ from the right.

Figure 6: The resulting combination of $\nu=I$ direct detection pathway and
$i=1$ interaction pathway.

$\displaystyle{\cal S}_{I,1}|_{a_{L}b_{R}}$ $\displaystyle=\int
dt_{4}dt_{3}dt_{2}dt_{1}\theta\left(t_{1}t_{2}\right)\theta\left(t_{2}t_{3}\right)\theta\left(t_{3}t_{4}\right)$
$\displaystyle\times\text{tr}\left\\{\left[E_{b,R}^{\dagger}\left(t+\tau\right)E_{b,L}\left(t+\tau\right)E_{a,R}^{\dagger}\left(t\right)E_{a,L}\left(t\right)\right]E_{a,L}^{\dagger}\left(t_{1}\right)E_{b,R}\left(t_{2}\right)E_{a,L}\left(t_{3}\right)E_{b,R}^{\dagger}\left(t_{4}\right)\rho_{\varphi}\left(-\infty\right)\right\\}$
$\displaystyle\times\text{tr}\left\\{V_{L}\left(t_{1}\right)V_{R}^{\dagger}\left(t_{2}\right)V_{L}^{\dagger}\left(t_{3}\right)V_{R}\left(t_{4}\right)\rho_{\mu}\left(-\infty\right)\right\\}.$

The operators in the square brackets correspond to the detection process and
thus last. Initially, only two fields modes are populated and thus we assume
that after the detection process, the field returns to its ground-state
(vacuum) and obtain

$\displaystyle{\cal S}_{I,1}|_{b_{L}a_{R}}$ $\displaystyle=\int
dt_{4}dt_{3}dt_{2}dt_{1}\theta\left(t_{1}t_{2}\right)\theta\left(t_{2}t_{3}\right)\theta\left(t_{3}t_{4}\right)$
$\displaystyle\times\langle\text{vac}|E_{b,L}\left(t+\tau\right)E_{a,L}\left(t\right)E_{b,L}^{\dagger}\left(t_{1}\right)E_{b,L}\left(t_{3}\right)|\Psi_{\theta}\rangle_{\varphi}\langle\Psi_{\theta}|E_{a,R}^{\dagger}\left(t_{4}\right)E_{a,R}\left(t_{2}\right)E_{a,R}^{\dagger}\left(t\right)E_{b,R}^{\dagger}\left(t+\tau\right)|\text{vac}\rangle_{\varphi}$
$\displaystyle\times\left\langle
V_{L}\left(t_{1}\right)V_{R}^{\dagger}\left(t_{2}\right)V_{L}^{\dagger}\left(t_{3}\right)V_{R}\left(t_{4}\right)\right\rangle_{\mu}.$

where we have plugged in the explicit expression for the source term
$\rho_{\varphi}\equiv|\Psi_{\theta}\rangle\langle\Psi_{\theta}|$. using
$\Phi_{\theta}\left(t_{1},t_{2}\right)=\frac{1}{\sqrt{2}}\left[\phi\left(t_{1},t_{2}\right)+e^{-i\theta}\phi\left(t_{2},t_{1}\right)\right]$.
Following Ref. [40], we change the integration time variables to time
differences between interaction events and obtain

$\displaystyle{\cal S}_{I,1}|_{b_{L}a_{R}}$ $\displaystyle=\int
d\tau_{3}d\tau_{4}$
$\displaystyle\times\Phi_{\theta}\left(t,t-\tau_{3}\right)\Phi_{\theta}^{*}\left(t,t-\tau_{3}-\tau_{4}\right)$
$\displaystyle\times\left\langle V_{L}{\cal
G}\left(\tau_{0}\right)V_{R}^{\dagger}{\cal
G}\left(\tau_{2}\right)V_{L}^{\dagger}{\cal
G}\left(\tau_{3}\right)V_{R}\right\rangle_{\mu},$

where the Liouville-space Green’s function is given by ${\cal
G}\left(t\right)=-\frac{i}{\hbar}\theta\left(t\right)e^{-i{\cal L}_{\mu}t-\eta
t}$. This terms was obtained from Eq. 10, using one additional approximation:
the free-photon propagator from the sample to the detector is taken to be
$\left\langle
E_{a,L}\left(t\right)E_{a,L}^{\dagger}\left(t_{1}\right)\right\rangle\approx\delta\left(\Delta
t-\nicefrac{{L}}{{c}}\right)$ where $\nicefrac{{L}}{{c}}$ is the distance
between the detector and the sample, $c$ is the speed of light and $\Delta t$
is taken to be the time difference between the emission and detection.

###  Appendix C The final signal – all combinations

Here we combine all possible contributions that correspond to all contributing
configurations of the detection and interaction pathways. All pathways are
summed as shown in Eq 10.

####  C.1 ${\cal S}_{\nu,1}$

$\displaystyle{\cal S}_{I,1}=$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t-\tau_{3:4},t+\tau_{0}\right)\Phi\left(t,t-\tau_{3}\right){\cal
F}_{1}\left(\tau_{0},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{II,1}:$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+\tau_{0},t-\tau_{3:4}\right)\Phi\left(t-\tau_{3},t\right){\cal
F}_{1}\left(\tau_{0},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{III,1}=-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t,t+\tau_{0}-\tau_{3:4}\right)\Phi\left(t+\tau_{0}-\tau_{3},t-T\right){\cal
F}_{1}\left(T,\tau_{3},\tau_{4}\right)$ $\displaystyle-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t-\tau_{3:4},t+\tau_{0}\right)\Phi\left(t-\tau_{3},t-T\right){\cal
F}_{1}\left(T+\tau_{0},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{IV,1}=-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+T+\tau,t-T-\tau_{3:4}\right)\Phi\left(t-T-\tau_{3},t+\tau_{0}\right){\cal
F}_{1}\left(T,\tau_{3},\tau_{4}\right)$ $\displaystyle-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+T+\tau,t-T-\tau_{3:4}\right)\Phi\left(t,t-T-\tau_{3}\right){\cal
F}_{1}\left(\tau_{0}+T,\tau_{3},\tau_{4}\right)$

####  C.2 ${\cal S}_{\nu,2}$

$\displaystyle{\cal S}_{I,2}=$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t-\tau_{4},t+\tau_{0}\right)\Phi\left(t,t+\tau_{3}\right){\cal
F}_{2}\left(\tau_{0}-\tau_{3},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{II,2}=$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+\tau_{0},t-\tau_{4}\right)\Phi\left(t+\tau_{3},t\right){\cal
F}_{2}\left(\tau_{0}-\tau_{3},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{III,2}$ $\displaystyle=-\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t,t+\tau_{0}-\tau_{4}\right)\Phi\left(t+\tau_{0}+\tau_{3},t-T\right){\cal
F}_{2}\left(T-\tau_{3},\tau_{3},\tau_{4}\right)$ $\displaystyle-\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t-\tau_{4},t+\tau_{0}\right)\Phi\left(t+\tau_{3},t-T\right){\cal
F}_{2}\left(\tau_{0}+T-\tau_{3},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{IV,2}=-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+T+\tau,t-T-\tau_{4}\right)\Phi\left(t-T+\tau_{3},t+\tau_{0}\right){\cal
F}_{2}\left(T-\tau_{3},\tau_{3},\tau_{4}\right)$ $\displaystyle-$
$\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+T+\tau,t-T-\tau_{4}\right)\Phi\left(t,t-T+\tau_{3}\right){\cal
F}_{2}\left(\tau_{0}+T-\tau_{3},\tau_{3},\tau_{4}\right)$

####  C.3 ${\cal S}_{\nu,3}$

$\displaystyle{\cal S}_{I,3}=$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t-\tau_{3},t+\tau_{0}\right)\Phi\left(t,t-\tau_{3:4}\right){\cal
F}_{3}\left(\tau_{0},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{II,3}=$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+\tau_{0},t-\tau_{3}\right)\Phi\left(t-\tau_{3:4},t\right){\cal
F}_{3}\left(\tau_{0},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{III,3}=-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t,t+\tau_{0}-\tau_{3}\right)\Phi\left(t+\tau_{0}-\tau_{3:4},t-T\right){\cal
F}_{3}\left(T,\tau_{3},\tau_{4}\right)$ $\displaystyle-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t-\tau_{3},t+\tau_{0}\right)\Phi\left(t-\tau_{3:4},t-T\right){\cal
F}_{3}\left(\tau_{0}+T,\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{IV,3}=-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+T+\tau,t-T-\tau_{3}\right)\Phi\left(t-T-\tau_{3:4},t+\tau_{0}\right){\cal
F}_{3}\left(T,\tau_{3},\tau_{4}\right)$ $\displaystyle-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+T+\tau,t-T-\tau_{3}\right)\Phi\left(t,t-T-\tau_{3:4}\right){\cal
F}_{3}\left(T+\tau_{0},\tau_{3},\tau_{4}\right)$

####  C.4 ${\cal S}_{\nu,4}$

$\displaystyle{\cal S}_{I,4}=$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t,t+\tau_{0}\right)\Phi\left(t-\tau_{4},t+\tau_{3}\right){\cal
F}_{4}\left(\tau_{0}-\tau_{3},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{II,4}$ $\displaystyle=\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+\tau_{0},t\right)\Phi\left(t+\tau_{0}+\tau_{3},t+\tau_{0}-\tau_{4}\right){\cal
F}_{4}\left(\tau_{0}-\tau_{3},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{III,4}=-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t,t+\tau_{0}\right)\Phi\left(t-T+\tau_{3},t-T-\tau_{4}\right){\cal
F}_{4}\left(2T+\tau_{0}-\tau_{3},\tau_{3},\tau_{4}\right)$

$\displaystyle{\cal S}_{IV,4}=-$ $\displaystyle\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+T+\tau,t-T\right)\Phi\left(t-\tau_{4},t+\tau_{3}\right){\cal
F}_{4}\left(\tau_{0}-\tau_{3},\tau_{3},\tau_{4}\right)$

####  C.5 ${\cal S}_{\nu,5}$

Here, each contribution is naturally symmetrized.

${\cal S}_{I,5}=\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t,t+\tau_{0}\right)\left[\Phi_{ab}\left(t-\tau_{3},t-\tau_{3:4}\right)+\Phi_{ba}\right]{\cal
F}_{5}\left(\tau_{0},\tau_{3},\tau_{4}\right)$

${\cal S}_{II,5}=\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+\tau_{0},t\right)\left[\Phi_{ab}\left(t-\tau_{3},t-\tau_{3:4}\right)+\Phi_{ba}\right]{\cal
F}_{5}\left(\tau_{0},\tau_{3},\tau_{4}\right)$

${\cal S}_{III,5}=-\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t,t+\tau_{0}\right)\left[\Phi_{ab}\left(t-T-\tau_{3:4},t-T-\tau_{3}\right)+\Phi_{ba}\right]{\cal
F}_{5}\left(\tau_{0}+2T,\tau_{3},\tau_{4}\right)$

${\cal S}_{IV,5}=-\int
d\tau_{3}d\tau_{4}\Phi^{*}\left(t+T+\tau_{0},t-T\right)\left[\Phi_{ab}\left(t-\tau_{3},t-\tau_{3:4}\right)+\Phi_{ba}\right]{\cal
F}_{5}\left(\tau_{0},\tau_{3},\tau_{4}\right)$


Results results All the experiments were performed on a single NVIDIA RTX
A6000 with 10752 CUDA cores and 48 GB RAM. In all the experiments different
values of seed are used to see how the results of the models change with the
split.

ogbl-ddi In this section, we test the performance of SNN on the ogbl-ddi
dataset. During the comparison, the learning rate and the weight decay are
fixed and equal to $0.003$ and $0.0005$, respectively. The model is trained
with Adam optimizer, and 300 epochs are performed, completing the training
phase in about 18 hours. Training our model takes quite a long time due to the
inherent complex structure of SNNs, and the operations performed make
computing the backpropagation algorithm less time efficient than other
architectures.

It is important to notice that the variance of our results, compared to our
competitors, decreases consistently. Our network is less affected by
variations in the parameters than the other architectures making our method
stable and, therefore, fitter to be used in different problems.

We run statistical tests on the results obtained by our models on this task to
assess the significance of the differences. We used the student’s unpaired
t-test between our predictions and the predictions of our competitors with
independent and identically distributed samples. By selecting $p=0.05$ we
demonstrated that our results are statistically significant if we compare them
against GraphSage, GCN, and Node2Vec.

Model | Validation Hits@20 | Test Hits@20
---|---|---
Node2Vec | 0.329$\pm$0.012 | 0.233$\pm$0.021
SEAL | 0.285$\pm$0.027 | 0.306$\pm$0.039
GCN† | 0.555$\pm$0.021 | 0.370$\pm$0.051
GraphSAGE † | 0.626$\pm$0.037 | 0.539$\pm$0.047
SNN (our) † | 0.632$\pm$0.036 | 0.555$\pm$0.044
|  |
Table : Results of the evaluation using the ogbl-ddi dataset. SNN has a 32%
improvement on the test hits with respect to the baseline and a 1.6%
improvement with respect to the best model. † means statistically
significantly (Student’s unpaired t-test, $\mathbf{p<0.05}$)

##  MovieLens 1M

In this section, we test the performance of SNN on MovieLens 1M. We used a
90-10 train-test split. 90% of the samples are used for training and 10% for
testing. This is the most common example of a split for this dataset. During
the comparison, the learning rate and the weight decay are fixed and equal to
$0.005$ and $0.0001$. The model is trained with a computational time of 210
minutes. The selected size for the embedding vector is 64; consequently, the
input size of SNN is 64. We use as a loss function the BPR loss, as explained
in chapter LABEL:method. For the sake of completeness, we also experimented
with different loss functions. In recommender systems, the most common example
of a loss function is the Root Mean Squared Error (RMSE) loss. We also tested
our approach using this loss function, but our experimental results show that
the F1-Score was lower. The main advantage of BPR over RMSE is on the
computational time. Table  shows that BPR is up to 50% faster than RMSE. We
also tested Binary Cross-Entropy (BCE) loss. In this case, the efficiency of
the method was comparable, but the performance in terms of F1-Score was much
worse. Thus, we also did not consider this loss when training our models.

Loss | #Layers | F1-Score@10 | Time
---|---|---|---
RMSE | 2 | 0.112 | 283
RMSE | 5 | 0.148 | 371
BCE | 2 | 0.087 | 198
BCE | 5 | 0.101 | 262
BPR | 2 | 0.163 | 163
BPR | 5 | 0.192 | 210
|  |  |
Table : Performance of the same model using different loss functions and
numbers of layers. The computational time is in minutes, $\downarrow$ is
better.

The choice of F1-Score@10 and F1-Score@100 as evaluation metrics is not
random. In fact, by choosing $N=10$ and $N=100$ we are working with tho
different sizes for the sets of relevant items for the users. The most common
recommender systems based on graphs have a high F1-Score@N for small values of
N. In our case also with high values of N the F1-Score is really high. In
Figure  it is possible to compare our result with other state-of-the-art
systems and we demonstrate that our solution is able to achieve good
performance also with high values of $N$. In this Image was not possible to
represent the values of F1-Score@N for values of $N>50$ because the other
approaches obtained very bad results. Finally, in Table  we can see that SNN
outperforms all the baselines in terms of F1-Score and Recall.

[width=]figure/grafico_tesi_2.png

Figure : F1-Score values on MovieLens 1M compared with state-of-the-art
solutions. The value of F1-Score for CAKR [cross_att], RippleNet [Wang_2018]
and Ripp-MKR [ripp_mkr] goes down extremely fast. The value of the F1-Score of
our model goes down less fast.

In Figure  we have tested the value of F1 when increasing the number of
returned recommendations, $K$. We can observe that F1-Score starts decreasing
for $K>60$. This is different from our competitors, where this value ranges
between 20 and 40. This experiment shows that this architecture can achieve
good results also when we produce a high number of recommended, and
potentially relevant for the users, items. For example, consider MovieLens.
The MovieLens dataset contains movie ratings, therefore, the higher the number
of movies that we are able to recommend the more satisfied the user.

[width=]figure/grafico_tesi_4.png

Figure : The value of F1-Score of our solution starts decreasing with a value of K between 60 and 70. Model | Precision@10 | Recall@10 | F1-Score@10 | Precision@100 | Recall@100 | F1-Score@100
---|---|---|---|---|---|---
MVAE† | 0.090 | 0.091 | 0.091 | 0.054 | 0.414 | 0.096
SVAE† | 0.144 | 0.125 | 0.134 | 0.069 | 0.494 | 0.121
CoFM† | 0.321 | 0.130 | 0.178 | 0.354 | 0.135 | 0.189
SNN (our) | 0.284 | 0.145 | 0.192 | 0.122 | 0.520 | 0.198
|  |  |  |  |  |
Table : Results of the evaluation using the MovieLens 1M dataset. SNN has a
5.4% improvement on F1@10 with respect to the baseline and a 0.7% improvement
with respect to the best model.
† is used when the results are taken from the original paper.

##  MovieLens 100K

Model | Precision@10 | Recall@10 | F1-Score@10 | Precision@20 | Recall@20 | F1-Score@20
---|---|---|---|---|---|---
LinUCB† | 0.321 | 0.135 | 0.192 | 0.286 | 0.203 | 0.238
SVD† | 0.343 | 0.163 | 0.222 | 0.282 | 0.245 | 0.272
DeepMF† | 0.344 | 0.160 | 0.219 | 0.314 | 0.244 | 0.262
MDP† | 0.357 | 0.168 | 0.228 | 0.313 | 0.251 | 0.278
SNN (our) | 0.210 | 0.287 | 0.243 | 0.265 | 0.309 | 0.285
|  |  |  |  |  |
Table : Results of the evaluation using the MovieLens 100k dataset. SNN has a
5.1% improvement on F1@10 with respect to the baseline and a 1.5% improvement
with respect to the best model.
† is used when the results are taken from the original paper.

In this section, we test the performance of SNN on MovieLens 100k. We used an
80-20 train-test split. 80% of the samples are used for training and 20% for
testing. The value of the hyperparameters is different from the previous
experiment: the learning rate is fixed at $0.001$ while the weight decay is
$0.0002$. The model is trained with Adam optimizer and 50 epochs are performed
with a computational time of 30 minutes. It is possible to see how our
solution outperforms all the baselines in terms of F1-Score. The results can
be seen in Table . One of the advantages of SNN is the possibility of
preventing over-smoothing. Deeply stacking the number of layers does not
decrease the performance. In a lot of GNNs the features become progressively
smoother with increased depth [https://doi.org/10.48550/arxiv.1905.1094]. In
our case, the best results are obtained using 5 layers, as shown in Figure .

[width=]figure/number_of_layers_variation.png

Figure : Values of F1-Score on MovieLens 100k depending on the number of
layers. Our network is robust with respect to over-smoothing problems that are
really common in GNNs

##  Book-Crossing

In this section, we test the performance of SNN on Book-Crossing. During the
comparison, the learning rate and the weight decay are fixed and equal to
$0.001$ and $0.0005$. The model is trained for 70 epochs with a computational
time of 251 minutes using as loss function the BPR loss. The computational
time is a bit higher with respect to the previous experiments because the
dataset is bigger and also contains a lot of side information, so the
processing time is higher. In this experiment, we obtained better results by
setting the selected size for the embedding vector to 32. As we can see in
Table  we outperform the state-of-the-art results both on Recall@10 and
Recall@100.

Model | Recall@10 | Recall@100
---|---|---
SVD† | 0.046 | 0.109
LibFM† | 0.062 | 0.124
RippleNet† | 0.074 | 0.127
KGNN-LS† | 0.082 | 0.149
SNN (our) | 0.093 | 0.177
|  |
Table : Results of the evaluation using the Book-Crossing dataset. SNN has a
2.8% improvement on Recall@100 with respect to the baseline. † is used when
the results are taken from the original paper.