# Spectrally-resolved four-photon interference of time-frequency entangled
photons

Sofiane Merkouche∗1, Valérian Thiel1, and Brian J. Smith1 Corresponding
author<EMAIL_ADDRESS>1 Department of Physics and Oregon Center for Optical,
Molecular, and Quantum Science, University of Oregon, Eugene, Oregon 97403,
USA

###### Abstract

Pairs of photons entangled in their time-frequency degree of freedom are of
great interest in quantum optics research and applications, due to their
relative ease of generation and their high capacity for encoding information.
Here we analyze, both theoretically and experimentally, the behavior of phase-
insensitive spectrally-resolved interferences arising from two pairs of time-
frequency entangled photons. At its core, this is a multimode entanglement
swapping experiment, whereby a spectrally resolved joint measurement on the
idler photons from both pairs results in projecting the signal photons onto a
Bell state whose form depends on the measurement outcome. Our analysis is a
thorough exploration of what can be achieved using time-frequency entanglement
and spectrally-resolved Bell-state measurements.

###### pacs:

##  I Introduction

With the advent of the quantum information age, it is well-established by now
that the encoding of quantum information into the degrees of freedom of light
is the key component of quantum communication networks Kimble (2008). Whereas
polarization and spatial-mode encoding benefit from ease of implementation,
they are prone to scrambling from environmental noise and optical-fiber
transmission, which undermines their suitability for long-distance
communication. Meanwhile, the time-frequency (TF) degree of freedom is more
robust in this regard (distinct colors stay distinct), and TF-encoding has
long been the standard for classical communication for this reason Brecht  _et
al._ (2015).

In the quantum regime, entanglement plays a key role in many protocols for
computation Jozsa and Linden (2003) and communication Ekert (1991); de
Riedmatten  _et al._ (2005). Furthermore, entanglement of photons has recently
been of great interest to the metrology and spectroscopy communities due to
its promise of enhancements in sensitivity beyond what is attainable in the
classical domain Raymer  _et al._ (2013). It comes as no surprise, then, that
the generation of pairs of photons in well-defined TF-entangled states is a
widely-researched area of quantum optics, and great strides have been made
over this terrain over the past two decades Ansari  _et al._ (2018). In
addition to state generation, harnessing the full capabilities of quantum
entanglement also requires the ability to projectively measure onto entangled
states. Indeed, such entangled measurements, of which the Bell-state
measurement is the prototype Vértesi and Navascués (2011), are nearly as
ubiquitous in quantum protocols as entangled states, most notably in quantum
teleportation Bennett  _et al._ (1993) and entanglement swapping Żukowski  _et
al._ (1993).

In this work, we describe theoretically, and demonstrate experimentally, a
novel and versatile TF entanglement-swapping scheme. Our scheme relies on the
multimode nature of TF-entanglement in pairs of photons generated in SPDC,
commonly denoted signal and idler. The central component of the setup is a
multimode, frequency-resolved Bell-state measurement (BSM), performed on idler
photons from two independent pairs. The BSM heralds the surviving signal
photons onto a pulsed Bell state whose central frequencies depend on the
result of the BSM. In this way we are able to herald multiple orthogonal Bell
pairs and verify entanglement in each pair, all derived from the same source
state and measurement scheme. In section II we outline the theory underlying
our work. In section III we describe the experimental setup. Then, in section
IV we describe our results, which show remarkable agreement with a simple and
intuitive Gaussian model using pure quantum states, before concluding in
section V. Finally, in the Appendix we cover the more technical details of our
work which would otherwise encumber the account of our main results.

##  II Theory

###  II.1 Four photon state

The entanglement swapping setup for this work is depicted schematically in
Fig. 1. The setup consists of two independent spontaneous parametric down
conversion (SPDC) sources, each of which generates pairs of photons into paths
labeled by the bosonic operators $\hat{a}_{n}$ for the signal and
$\hat{b}_{n}$ for the idler, where $n\in\\{1,2\\}$ labels the two sources.
Here we consider pulsed collinear type II SPDC sources, where the signal and
idler modes are distinguished through polarization, and all the light is
collected into single-mode optical fibers, so that only the time-frequency of
freedom is relevant.

In the low gain regime, the output state of $n^{\textrm{th}}$ SPDC source is
given by

$\displaystyle\ket{\psi_{n}}=\sum_{k=0}^{\infty}\frac{\sqrt{\eta_{n}}^{k}}{k!}\left(\int\textrm{d}\omega_{\text{S}}\textrm{d}\omega_{\text{I}}f_{n}(\omega_{\text{S}},\omega_{\text{I}})\hat{a}^{\dagger}_{n}(\omega_{\text{S}})\hat{b}^{\dagger}_{n}(\omega_{\text{I}})\right)^{k}\ket{\mathrm{vac}},$
(1)

The function $f_{n}(\omega_{\text{S}},\omega_{\text{I}})$ is the joint
spectral amplitude (JSA), given by

$\displaystyle
f_{n}(\omega_{\text{S}},\omega_{\text{I}})=u_{n}(\omega_{\text{S}}+\omega_{\text{I}})\textrm{sinc}\left[\frac{\Delta
k_{n}(\omega_{\text{S}},\omega_{\text{I}})L}{2}\right],$ (2)

where $u_{n}$ represents the spectral mode function of the pump, $\Delta
k_{n}$ is the wave-vector mismatch between the pump, signal and idler modes,
and $L$ is the length of the interaction medium. Finally, $\eta_{n}$ is the
gain of the parametric process, which depends on the length $L$, the non-
linear strength of the material and the number of photons in the pump beam.

The SPDC state due to the two independent sources can be written as a tensor
product $\ket{\psi_{1}}\otimes\ket{\psi_{2}}$. We expand this and keep only
terms of order $\eta$, which are responsible for the four-photon contribution,
obtaining the following normalized state:

$\begin{gathered}\ket{\psi_{\mathrm{\eta}}}=\frac{\left(\sqrt{\eta_{1}\eta_{2}}\ket{\psi_{12}}+\frac{\eta_{1}}{2}\ket{\psi_{11}}+\frac{\eta_{2}}{2}\ket{\psi_{22}}\right)}{\sqrt{\eta_{1}\eta_{2}+\eta_{1}^{2}/4+\eta_{2}^{2}/4}},\end{gathered}$
(3)

where

$\displaystyle\ket{\psi_{nm}}=$ $\displaystyle\int\textrm{d}^{4}\omega\
f_{n}(\omega_{\text{S}},\omega_{\text{I}})f_{m}(\omega_{\text{S}}^{\prime},\omega_{\text{I}}^{\prime})\times$
$\displaystyle\hat{a}^{\dagger}_{n}(\omega_{\text{S}})\hat{a}^{\dagger}_{m}(\omega_{\text{S}}^{\prime})\hat{b}^{\dagger}_{n}(\omega_{\text{I}})\hat{b}^{\dagger}_{m}(\omega_{\text{I}}^{\prime})\ket{\text{vac}}.$
(4)

is the four-photon state arising from either a photon pair from each source,
or two pairs from one source and none from the other.

In order to facilitate the discussion of our results, we will make a few
simplifying assumptions, and we will relax these assumptions whenever
appropriate in the main text and discuss them further in the Appendix. First
we will assume that both sources are identical, such that
$f_{1}(\omega_{\text{S}},\omega_{\text{I}})=f_{2}(\omega_{\text{S}},\omega_{\text{I}})=f(\omega_{\text{S}},\omega_{\text{I}})$
and $\eta_{1}=\eta_{2}=\eta$. This is because in the experiment, the two
sources are derived from double-pumping the same crystal as in Ref. Pan  _et
al._ (1998), and are matched to a great degree as discussed in the
experimental section. Further, any source mismatch does not reduce the quality
of the entanglement in the swapped state, just the visibility of quantum
interference in the method we use to verify entanglement.

Second, we will assume that the relative phases between the three terms in (3)
are random, that is, these terms are not mutually coherent do not contribute
to any quantum interference. In reality there is such a coherence, which is
due to the phase of the pump, and we observe this in both two- and four-photon
interference, as we show in Appendix B. However, this phase drifts over the
course of the interference measurements we will describe, which are several
hours long, and thus the $\ket{\psi_{11}}$ and $\ket{\psi_{22}}$ contribute to
a constant background in these measurements.

Third, in the main text we will focus our attention solely on the
$\ket{\psi_{12}}$ term of (3), which corresponds to each source producing a
pair of entangled photons. In most entanglement swapping and quantum
teleportation experiments relying on SPDC sources, the other two terms (the
double-pair terms), though present, do not contribute when the four photons
are detected in coincidence, so that heralded states are post-selected
Wagenknecht  _et al._ (2010). Similarly, the double-pair terms only contribute
in our setup in the aforementioned interference measurement, where they
constitute a constant background which we measure and subtract.

Finally, as most of the experiments that we performed rely on performing a
spectrally-resolved Bell state measurement (BSM) on the idler photon, we will
assume for simplicity that our spectral resolution is unlimited. Although it
introduces some slight abuses of notation, this limit has the critical benefit
of providing a simple and intuitive model with which to understand the physics
in terms of pure quantum states (the pure-state approximation). A more
complete model taking into account the finite spectral resolution of the BSM
and the resulting states is then easily constructed from this pure-state
approximation in Appendix A. As we shall see, the pure-state approximation is
sufficient to account for the majority of the physics and the results of our
experiment.

Figure 1: Conceptual scheme of the experiment. Two sources, 1 and 2, emit
photon pairs, labelled $\hat{a}$ and $\hat{b}$ for the signal and idler
photons, respectively. The photon pairs from the first source undergo a delay
$\tau$ relative to their counterpart on source 1.

###  II.2 Heralded state and JSI

A Bell-state measurement is performed on the idler photons by resolving their
frequencies at the output of a 50:50 beamsplitter, with a (small) path
difference at the input giving a relative time delay $\tau_{\text{I}}$. The
beamsplitter transformation is then defined by the relations

$\begin{gathered}\hat{c}^{\dagger}(\Omega)=\frac{e^{i\Omega\tau_{\text{I}}}\hat{b}^{\dagger}_{1}(\Omega)+\hat{b}^{\dagger}_{2}(\Omega)}{\sqrt{2}},\\\
\hat{d}^{\dagger}(\Omega^{\prime})=\frac{e^{i\Omega^{\prime}\tau_{\text{I}}}\hat{b}^{\dagger}_{1}(\Omega^{\prime})-\hat{b}^{\dagger}_{2}(\Omega^{\prime})}{\sqrt{2}},\end{gathered}$
(5)

and coincidences are detected between $\hat{c}$ and $\hat{d}$. Note that for
clarity, we will label as $\Omega$ the optical frequencies of the idlers and
as $\omega$ those of the signals. Because we use the pure state approximation,
we can consider that this spectral measurement is achieved with perfect
resolution, such that the POVM element for this detection is given by

$\begin{gathered}\hat{\Pi}_{jk}^{\textrm{BSM}}=\ket{\Omega_{j},\Omega_{k}}\bra{\Omega_{j},\Omega_{k}},\\\
\textrm{with
}\ket{\Omega_{j},\Omega_{k}}=\hat{c}^{\dagger}(\Omega_{j})\hat{d}^{\dagger}(\Omega_{k})\ket{\text{vac}},\end{gathered}$
(6)

which is a projector onto the monochromatic idler frequencies $\Omega_{j}$ and
$\Omega_{k}$. We then proceed to compute the heralded signal state,
$\ket{\Psi_{jk}}$, defined by:

$\displaystyle\ket{\Psi_{jk}}\otimes\ket{\Omega_{j},\Omega_{k}}=\frac{\hat{\Pi}_{jk}^{\textrm{BSM}}\ket{\psi_{12}}}{\sqrt{p_{jk}}},$
(7)

where the norm $p_{jk}$ is given by

$\displaystyle
p_{jk}=\bra{\psi_{12}}\hat{\Pi}_{jk}^{\textrm{BSM}}\ket{\psi_{12}}.$ (8)

To compute the form of $\ket{\Psi_{jk}}$, we will introduce, for convenience,
the density matrices for the reduced states of the signal and idler photons,
given respectively by,

$\begin{gathered}\rho_{\mathrm{S}}(\omega,\omega^{\prime})=\int\textrm{d}\Omega
f(\omega,\Omega)f^{*}(\omega^{\prime},\Omega),\\\
\rho_{\mathrm{I}}(\Omega,\Omega^{\prime})=\int\textrm{d}\omega
f(\omega,\Omega)f^{*}(\omega,\Omega^{\prime}).\end{gathered}$ (9)

Since the sources are identical, the reduced density matrices obey the
relation $\rho(\omega,\omega^{\prime})=\rho^{\ast}(\omega^{\prime},\omega)$.
Using these, it is straightforward to show that

$\displaystyle
p_{jk}=\frac{1}{2}\Big{[}\rho_{\mathrm{I}}(\Omega_{j},\Omega_{j})\rho_{\mathrm{I}}(\Omega_{k},\Omega_{k})-|\rho_{\mathrm{I}}(\Omega_{j},\Omega_{k})|^{2}\cos{\theta_{jk}}\Big{]},$
(10)

with

$\displaystyle\theta_{jk}=(\Omega_{j}-\Omega_{k})\tau_{\text{I}}.$ (11)

This is the probability distribution for a coincidence between the idler
photons at $(\Omega_{j},\Omega_{k})$, or equivalently, the JSI of the idlers
at the output of the beamsplitter.

Most importantly, the heralded signal Bell-state has the simple form

$\displaystyle\ket{\Psi_{jk}}=\frac{\ket{\phi_{j}}_{1}\ket{\phi_{k}}_{2}-e^{i\theta_{jk}}\ket{\phi_{k}}_{1}\ket{\phi_{j}}_{2}}{\sqrt{2\mathcal{C}_{jk}}},$
(12)

where $\ket{\phi_{j(k)}}_{1(2)}$ is a pulse-mode normalized single photon
state Blow  _et al._ (1990) given by

$\ket{\phi_{j(k)}}_{1(2)}=\int\textrm{d}\omega\,\phi_{j(k)}(\omega)\hat{a}^{\dagger}_{1(2)}(\omega)\ket{\text{vac}},$
(13)

with

$\phi_{j(k)}(\omega)=\frac{f(\omega,\Omega_{j(k)})}{\sqrt{\rho_{\mathrm{I}}(\Omega_{j(k)},\Omega_{j(k)})}}.$
(14)

Essentially, $\ket{\phi_{j(k)}}$ is the state that the signal photon is
projected onto, when its corresponding idler photon is detected at frequency
$\Omega_{j(k)}$. Finally the normalization constant $\mathcal{C}_{jk}$ is
given by

$\displaystyle\mathcal{C}_{jk}=1-\left|\braket{\phi_{j}}{\phi_{k}}\right|^{2}\cos{\theta_{jk}}.$
(15)

The functions $\phi_{j(k)}$ are defined from the JSA and several identities
are shown in Appendix F. Note that, while normalized, they are not necessarily
orthogonal, hence the dependency on the modal overlap in the normalization
from Eq.(15). In Sec.IV, we will use a Gaussian approximation that gives a
simple expression for those functions.

Since nearly all of our measurements are conditioned upon the BSM on the idler
photons, we will use the state $\ket{\Psi_{jk}}$ to calculate our quantities
of interest. In addition, we also consider the case of a non-resolving BSM,
given by

$\displaystyle\hat{\Pi}^{\text{BSM}}$
$\displaystyle=\sum_{j,k}\hat{\Pi}^{\text{BSM}}_{jk}$ (16)
$\displaystyle=\int\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}\ket{\Omega_{j},\Omega_{k}}\bra{\Omega_{j},\Omega_{k}},$

which, upon taking
$\text{Tr}\Big{(}\hat{\Pi}^{\text{BSM}}\ket{\psi_{12}}\bra{\psi_{12}}\Big{)}$,
heralds the mixed state

$\hat{\rho}=\sum_{jk}p_{jk}\ket{\Psi_{jk}}\bra{\Psi_{jk}},$ (17)

and the quantities that arise from this state. Note that this state is not
normalized, since $\sum_{jk}p_{jk}\neq 1$, but this fact is inconsequential to
our measurements, as only the relative probabilities $p_{jk}$ are physically
relevant. At this point we note that we are using the sum $\sum_{jk}$ loosely
in place of the integral $\int\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}$
because of the finite spectral resolution of our heralding BSM. This means
that the pure state description of $\ket{\Psi_{jk}}$ is an approximation, but
one whose validity we justify throughout the paper as well as in Appendix A
and also gives similar results as the model which takes into account the
finite resolution of the BSM and the resultant impurity of the heralded state.

###  II.3 State characterization and entanglement verification

We characterize the heralded state $\ket{\Psi_{jk}}$ first by measuring its
joint spectral intensity (JSI). This measurement is defined by the POVM

$\begin{gathered}\hat{\Pi}^{\text{JSI}}=\ket{\omega_{1},\omega_{2}}\bra{\omega_{1},\omega_{2}},\\\
\textrm{with
}\ket{\omega_{1},\omega_{2}}=\hat{a}^{\dagger}_{1}(\omega_{1})\hat{a}^{\dagger}_{2}(\omega_{2})\ket{\text{vac}},\end{gathered}$
(18)

and we calculate the resultant JSI as

$\displaystyle F_{jk}$
$\displaystyle(\omega_{1},\omega_{2})=\bra{\Psi_{jk}}\hat{\Pi}^{\text{JSI}}\ket{\Psi_{jk}}=\left|\Braket{\omega_{1},\omega_{2}}{\Psi_{jk}}\right|^{2}.$
(19)

Using (12), we obtain explicitly

$\displaystyle
F_{jk}(\omega_{1},\omega_{2})=\frac{1}{2\mathcal{C}_{jk}}\Big{|}\phi_{j}(\omega_{1})\phi_{k}(\omega_{2})-e^{i\theta_{jk}}\phi_{j}(\omega_{2})\phi_{k}(\omega_{1})\Big{|}^{2}.$
(20)

In the absence of spectral resolution in the BSM, the heralded state is
$\hat{\rho}$ from Eq.(17), and the measured JSI is given by

$\displaystyle F(\omega_{1},\omega_{2})$
$\displaystyle=\text{Tr}\left(\hat{\rho}\ \hat{\Pi}^{\text{JSI}}\right)$
$\displaystyle=\sum_{jk}p_{jk}\left|\Braket{\omega_{1},\omega_{2}}{\Psi_{jk}}\right|^{2}$
$\displaystyle=\sum_{jk}p_{jk}F_{jk}(\omega_{1},\omega_{2}).$ (21)

In terms of previously defined quantities, this function takes the general
form

$\displaystyle F(\omega_{1},\omega_{2})=$
$\displaystyle\frac{1}{2}\Bigg{[}\rho_{\mathrm{S}}(\omega_{1},\omega_{1})\rho_{\mathrm{S}}(\omega_{2},\omega_{2})$
$\displaystyle-\left|\int\textrm{d}\Omega
f(\omega_{1},\Omega)f^{*}(\omega_{2},\Omega)e^{i\Omega\tau_{\text{I}}}\right|^{2}\Bigg{]}.$
(22)

To verify that the heralded state $\ket{\Psi_{jk}}$ is indeed entangled,
beyond classical correlation, two-photon interference is used in a manner
analogous to Refs. Ramelow  _et al._ (2009); Graffitti  _et al._ (2020). Here
the signal photons are combined at another 50:50 beamsplitter while scanning a
relative delay between the input modes ($\hat{a}_{1}$ and $\hat{a}_{2}$)
denoted by $\tau_{\text{S}}$ and monitoring coincidences at the output modes
($\hat{x}$ and $\hat{y}$) as depicted in Fig. 1. These modes transform as:

$\displaystyle\hat{x}^{\dagger}(\omega)=\frac{e^{i\omega\tau_{\text{S}}}\hat{a}^{\dagger}_{1}(\omega)+\hat{a}^{\dagger}_{2}(\omega)}{\sqrt{2}},$
$\displaystyle\hat{y}^{\dagger}(\omega^{\prime})=\frac{e^{i\omega^{\prime}\tau_{\text{S}}}\hat{a}^{\dagger}_{1}(\omega^{\prime})-\hat{a}^{\dagger}_{2}(\omega^{\prime})}{\sqrt{2}}.$
(23)

The POVM associated with such a coincidence detection is defined as

$\displaystyle\begin{gathered}\hat{\Pi}_{\textrm{verif}}=\int\textrm{d}\omega\textrm{d}\omega^{\prime}\ket{\omega,\omega^{\prime}}\bra{\omega,\omega^{\prime}},\\\
\textrm{with
}\ket{\omega,\omega^{\prime}}=\hat{x}^{\dagger}(\omega)\hat{y}^{\dagger}(\omega^{\prime})\ket{\text{vac}},\end{gathered}$
(26)

and the probability of detecting a coincidence heralded for the input state
$\ket{\Psi_{jk}}$ is given by

$\displaystyle P_{jk}(\tau_{\text{S}})$
$\displaystyle=\Braket{\Psi_{jk}}{\hat{\Pi}_{\textrm{verif}}}{\Psi_{jk}}$
$\displaystyle=\int\textrm{d}\omega\textrm{d}\omega^{\prime}\,\Big{|}\braket{\omega,\omega^{\prime}}{\Psi_{jk}}\Big{|}^{2}.$
(27)

Evaluating this using Eq. (12), we obtain

$\displaystyle P_{jk}(\tau_{\text{S}})=\frac{1}{2\mathcal{C}_{jk}}\Big{(}1$
$\displaystyle+\mathcal{F}_{j}(\tau_{\text{S}})\mathcal{F}_{k}(\tau_{\text{S}})\cos{\left[(\omega_{j}-\omega_{k})\tau_{\text{S}}-\theta_{jk}\right]}$
$\displaystyle-\mathcal{O}(\left|\braket{\phi_{j}}{\phi_{k}}\right|^{2})\Big{)}.$
(28)

where $\mathcal{F}_{j(k)}(\tau_{\text{S}})$ is the Fourier transform of
$|\phi_{j(k)}(\omega)|^{2}$, $\omega_{j(k)}$ is the first moment of
$\phi_{j(k)}$, and finally,

$\displaystyle\mathcal{O}(\left|\braket{\phi_{j}}{\phi_{k}}\right|^{2})$
$\displaystyle=\left|\int\textrm{d}\omega\phi_{j}^{*}(\omega)\phi_{k}(\omega)e^{i\omega\tau_{\text{S}}}\right|^{2}$
$\displaystyle+\left|\braket{\phi_{j}}{\phi_{k}}\right|^{2}\cos{\theta_{jk}}$
(29)

are terms that depend on the overlap of $\phi_{j}$ and $\phi_{k}$ and are
negligible except for when $\Omega_{j}\approx\Omega_{k}$, which is the regime
where $p_{jk}\approx 0$. The main feature of $P_{jk}(\tau_{\text{S}})$ is the
fringes from the oscillating term at the difference frequency
$(\omega_{j}-\omega_{k})$, which is a signature of frequency-bin entanglement
Ramelow  _et al._ (2009); Fedrizzi  _et al._ (2009), as these fringes arise
due to the coherence between the two terms in the state $\ket{\Psi_{jk}}$.

In the absence of frequency resolution of the BSM, we again use $\hat{\rho}$
to obtain

$\displaystyle P(\tau_{\text{S}})$ $\displaystyle=\text{Tr}(\hat{\rho}\
\hat{\Pi}_{\textrm{verif}})$
$\displaystyle=\sum_{jk}p_{jk}\int\textrm{d}\omega\textrm{d}\omega^{\prime}\left|\Braket{\omega,\omega^{\prime}}{\Psi_{jk}}\right|^{2}$
$\displaystyle=\sum_{jk}p_{jk}P_{jk}(\tau_{\text{S}}).$ (30)

In terms of the previously defined quantities, we have

$\displaystyle P(\tau_{\text{S}},\tau_{\text{I}})$
$\displaystyle=\frac{1}{4}\Bigg{(}1+\left|\int\textrm{d}\omega\textrm{d}\Omega\,\left|f(\omega,\Omega)\right|^{2}e^{i(\omega\tau_{\text{S}}+\Omega\tau_{\text{I}})}\right|^{2}$
$\displaystyle-\int\textrm{d}^{2}\Omega\,\left|\rho_{\mathrm{I}}(\Omega,\Omega^{\prime})\right|^{2}e^{i(\Omega-\Omega^{\prime})\tau_{\text{I}}}$
$\displaystyle-\int\textrm{d}^{2}\omega\,\left|\rho_{\mathrm{S}}(\omega,\omega^{\prime})\right|^{2}e^{i(\omega-\omega^{\prime})\tau_{\text{S}}}\Bigg{)},$
(31)

where we have written the dependence of $\theta_{jk}$ on $\tau_{\text{I}}$
explicitly.

Eqs.(21) and (30) show that summing over the individual “pixellized”
quantities requires a specific weight or gain. The results are naturally
equivalent to replacing the spectrometers with bucket detectors. This process
may also be interpreted in term of probabilities, where the probability of
getting a four fold coincidence is the product of the probability of a
coincidence between the heralding idler photons by the conditional probability
of a coincidence between the heralded signal photons. This concept of
weighting the average is similar to other experiments that utilize multipixel
detection in the spatial Treps  _et al._ (2005) or spectral Roslund  _et al._
(2014); Thiel  _et al._ (2017) domain.

##  III Experiment

###  III.1 Description

Figure 2: a) Experimental setup; PD: photodiode; SHG: second harmonic
generation; DM: dichroic mirror; PDC: parametric down conversion; PZT:
piezoelectric actuator; PBS: polarizing beamsplitter; TOFS: time of flight
spectrometer; SMF: single mode fiber; CFBG: chirped fiber Bragg grating; FBS:
fiber beamsplitter. b) Bell state measurement: the idler photons interferences
are spectrally-resolved and used as a herald. c-d) different configuration of
the measurement on the signal photons depending on the experiment.

The BSM is achieved by interfering the idler photons from each source on a
fiber beamsplitter which is used for heralding. Spectral resolution is then
gained by utilizing frequency-to-time conversion at the output of the
beamsplitter, thus heralding the frequency bins $\Omega_{j}$ and $\Omega_{k}$,
see Fig.2b). Finally, the heralded JSI and verification is done respectively
by routing the signal photons through the setup described in Fig.2c-d). An
adjustable delay is introduced independently on the signal and idler paths of
source 1 using delay lines to match the time of arrival of all the photons.

###  III.2 Detection

The single photon are detected utilizing superconducting nanowire single
photon detector (SNSPD) from IDQuantique (ID281) which can detect the arrival
time of photons with a resolution of 20 ps. This temporal resolution is
translated into spectral resolution using time-of-flight spectrometers (TOFS),
thanks to frequency-to-time conversion Torres-Company  _et al._ (2011); Goda
and Jalali (2013). For coarse spectral resolution, we used two spools of 500
meters-long HP780 fiber. These imprint a dispersion of approximately 50 ps/nm
and the losses per spool at 830 nm are less than 1 dB. For fine resolution, we
used two chirped fiber Bragg gratings (CFBG from Teraxion) with a dispersion
of 1000 ps/nm Davis  _et al._ (2017). This extra resolution comes with heavy
losses of over 10 dB from technical origins and from a finite spectral window
of 10 nm. The photocurrent coming out of the detectors are registered with a
time-to-digital converter (TDC, ID900 from IDQuantique). The time reference is
provided by the clock generated by an internal photodiode in the laser source,
thus ensuring that each time tag is taken with respect to a stable signal for
each pulse. With this setup, it is possible to register coincidences between
any combination of the four photons with spectral resolution. This allows for
the “pixelization” of any event into spectral bins in post-processing.

Figure 3: Acquisiton of the heralded JSI $p_{jk}$ of the idlers photon showing
the calibration and indexing convention of the frequency bins.

The data for the main experiments is acquired setting the time tagger’s base
unit at 100 ps, resulting in a spectral resolution of 0.1 nm when used with
the CFBG and 2 nm with the fiber spools. For experiments that only require
spectral resolution on the herald (such as the verification $P_{jk}$), time
tags corresponding to heralding frequencies $\Omega_{j},\Omega_{k}$ are
acquired and subsequently histogrammed at the resolution of the spectrometer,
corresponding to the probability of getting spectral coincidences in every
possible combination of frequency bins. When all four spectrometers are needed
(for instance to measure the JSI $F_{jk}$ heralded by a BSM at frequencies
$\Omega_{j}$ and $\Omega_{k}$), the measurement then consists of four sets of
time tags that can be binned into a four-dimensional histogram, such as shown
in Fig.6. For more readability, we label those heralding bins by integers $j$
and $k$, such that index $0$ corresponds to $\Omega_{0}$, the center
wavelength of the idlers’ joint spectrum. The axes in Fig.3 showcases the bin
index (called pixel) to wavelength mapping.

###  III.3 Source distinguishability

This experiment relies critically on the indinstinguishability between the two
sources as hinted in the previous section. This has to be achieved on every
degree of freedom. Since the polarization and the spatial degree of freedom
are constrained by polarization-maintaining fibers, there remains to match the
sources in frequency and in time.

The spectral indistinguishability can be estimated by measuring the joint
spectral intensity of both sources. This is steadily performed by connecting
our spectrometers to the ports corresponding to $\hat{a}^{\dagger}_{n}$ and
$\hat{b}^{\dagger}_{n}$ in Fig.2a), where $n=1,2$ is the source number, and
acquiring the two-fold spectral coincidences. The JSIs, shown in Fig.4, were
obtained using the fiber spool as spectrometers but operating the time tagger
at its maximum resolution of 13 ps, thus limiting spectral resolution by the
timing jitter of the superconducting nanowires. We can see that both sources
are very similar thanks to the dual pass configuration of the pump. This
experimental JSI is used throughout this paper to define the parameters of the
JSA $f(\omega,\Omega)$ that are used in the simulation. Note that there can be
an additional cause of distinguishability due to spectral phase mismatch
between both sources which cannot be determined with an intensity measurement.
Therefore, the dispersion was mostly matched in every path of the
interferometer by ensuring that every fiber element had the same length. Any
remaining dispersion would from free space propagation has a negligible
effect.

Figure 4: Experimental joint spectral intensity of both sources. Insets:
marginal spectra.

The time indinstinguishability is achieved by scanning both delays
$\tau_{\text{S}}$ and $\tau_{\text{I}}$ while monitoring the coincidences
between ports $\hat{c}$ and $\hat{d}$ for the idlers and $\hat{x}$ and
$\hat{y}$ for the signals, which statistics are uncorrelated. This results in
a Hong Ou Mandel dip between the thermal states obtained by tracing over the
other photon. The dips have a poor visibility (less than 10%) which is still
sufficient to coarsely match the time of arrival of the photons. Note that the
same HOM interference was also realized by heralding, increasing the
visibility. In Appendix C, we show how such measurement with and without
spectral resolution can give a lower bound on the purity of the heralded
state.

Finally, another convenient method Kim  _et al._ (2015) to match both sources
can be applied. With the configuration from Fig.1, we can monitor the
coincidences between either port of each beamsplitter, for instance $\hat{c}$
and $\hat{x}$, which share photon number correlations, as opposed to the
previous uncorrelated case. When both delays are matched, i.e.
$\tau_{\text{S}}=\tau_{\text{I}}=0$, then these coincidence oscillate at the
optical frequency. This is akin to the classical case except that the
interfering term involves a two-partite system. In Appendix B, we show that
the visibility of these interferences provides a direct measurement of the
overlap between both sources, taking into account any phase effect. This
method proved essential to accurately match delays before every experiment,
while it also provided a bound to quantify indistinguishability, with a
maximum measured contrast of $80\%$.

##  IV Simulations and results

In this section, we further model the experiment with a Gaussian approximation
of the JSA to derive analytical expressions from the quantities defined in
Sec. II. This allows for a better understanding of the dependency of the
interferences on the experimental parameters, notably on the delay between the
idler photons in the BSM. We then compare our experimental results to the
theory using this approximation. We first consider the case of a perfect BSM
with $\tau_{\text{I}}=0$ and explore the effect of minor distinguishability in
the BSM with $\tau_{\text{I}}\neq 0$.

###  IV.1 Gaussian model

It is convenient to write the JSA $f(\omega_{\text{S}},\omega_{\text{I}})$ as
a Gaussian distribution by approximating the Sinc function with a Gaussian of
the same width :

$\displaystyle f(\omega_{\text{S}},\omega_{\text{I}})=C\,\exp\Bigg{[}$
$\displaystyle-\left(\frac{\omega_{\text{S}}-\omega_{0}}{2\sigma_{\text{S}}}\right)^{2}-\left(\frac{\omega_{\text{I}}-\omega_{0}}{2\sigma_{\text{I}}}\right)^{2}$
$\displaystyle-\alpha(\omega_{\text{S}}-\omega_{0})(\omega_{\text{I}}-\omega_{0})\Bigg{]},$
(32)

where $\sigma_{s}$ ($\sigma_{\text{I}}$) is the cross-sectional width of the
JSI in the $\omega_{t}ext{S}$\- ($\omega_{\text{I}}$-) direction evaluated at
$\omega_{\text{I}}$ ($\omega_{\text{S}}$), $\omega_{0}$ is the center
frequency, $\alpha$ quantifies the amount of spectral entanglement, and
$C=\left(\int\textrm{d}^{2}\omega\left|f(\omega,\omega^{\prime})\right|^{2}\right)^{-1/2}$
is a normalization constant.

From this definition, the density matrices $\rho_{\mathrm{S}}$ and
$\rho_{\mathrm{I}}$ are given by

$\begin{gathered}\rho_{\mathrm{S}}(\omega,\omega^{\prime})=C^{2}\sqrt{2\pi}\sigma_{\text{I}}\exp\Bigg{[}-\frac{(\omega-\omega_{0})^{2}+(\omega^{\prime}-\omega_{0})^{2}}{4\sigma_{\text{S}}^{2}}\\\
+\frac{1}{2}\alpha^{2}\sigma_{\text{I}}^{2}(\omega+\omega^{\prime}-2\omega_{0})^{2}\Bigg{]},\end{gathered}$
(33)
$\begin{gathered}\rho_{\mathrm{I}}(\Omega,\Omega^{\prime})=C^{2}\sqrt{2\pi}\sigma_{\text{S}}\exp\Bigg{[}-\frac{(\Omega-\omega_{0})^{2}+(\Omega^{\prime}-\omega_{0})^{2}}{4\sigma_{\text{I}}^{2}}\\\
+\frac{1}{2}\alpha^{2}\sigma_{\text{S}}^{2}(\Omega+\Omega^{\prime}-2\omega_{0})^{2}\Bigg{]}.\end{gathered}$
(34)

Figure 5: Gaussian approximation of the spectral heralding: a detection at
frequency $\Omega_{j}$ projects the signal photon onto the mode $\phi_{j}$,
centered at $\omega_{j}$, shown on the left (dashed). The pure state model
assumes that the idler photon is detected with perfect resolution and does not
have any spectral support (top, dashed). A more accurate development shown in
Appendix A considers integration over a finite spectral window for the herald
(top, plain) and the resulting heralded mode (left, plain). In this case,
their first moment slightly differ from $\omega_{j}$ and $\Omega_{j}$ obtained
with the approximation.

Meanwhile, the $\phi_{j(k)}$ functions are given by

$\phi_{j(k)}(\omega)=\frac{1}{\sqrt{\sigma_{\text{S}}\sqrt{2\pi}}}\exp[-(\omega-\omega_{j(k)})^{2}/4\sigma_{\text{S}}^{2}],$
(35)

which are Gaussians with a width equal to that of the signal’s cross-sectional
width and a central frequency given by

$\displaystyle\omega_{j(k)}=\omega_{0}-2\alpha\ \sigma_{\text{S}}^{2}\
\Omega_{j(k)}.$ (36)

Finally, the normalization constant $\mathcal{C}_{jk}$ is given by

$\mathcal{C}_{jk}=1-\exp\Big{[}-(\omega_{j}-\omega_{k})^{2}/4\sigma_{\text{S}}^{2}\Big{]}\cos(\Omega_{j}-\Omega_{k})\tau_{\text{I}}.$
(37)

Note that, for most of our data,
$(\omega_{j}-\omega_{k})^{2}\gg\sigma_{\text{S}}^{2}$, so that
$\mathcal{C}_{jk}\simeq 1$. In Fig.5, we used the experimental JSI to simulate
its marginals with a finite resolution. The heralded mode is centered at
frequency $\omega_{j}$ and we plotted in dashed lines
$\left|\phi_{j}\right|^{2}$ showing that the Gaussian approximation is quite
accurate. We remind that we consider for the simulations that the spectral
resolution along the idlers frequency is infinite.

###  IV.2 Case of $\tau_{\text{I}}=0$ ($\theta_{jk}=0)$

We first study the case where the idler paths are exactly matched, setting
$\tau_{\text{I}}=\theta_{jk}=0$, which yields the results that we report on in
Merkouche  _et al._ (2021). Under this constraint, and using the Gaussian
model in (32), we proceed to obtain analytic forms for our measured
quantities. First, the probability $p_{jk}$ of performing a BSM at frequencies
$\Omega_{j}$ and $\Omega_{k}$ and heralding the state $\ket{\Psi_{jk}}$ is
obtained from (10):

$\displaystyle p_{jk}=\frac{1}{2}\Big{[}$
$\displaystyle\rho_{\mathrm{I}}(\Omega_{j},\Omega_{j})\rho_{\mathrm{I}}(\Omega_{k},\Omega_{k})-\left|\rho_{\mathrm{I}}(\Omega_{j},\Omega_{k})\right|^{2}\Big{]}.$
(38)

This in fact corresponds to the joint probability of detecting two photons at
the output of a beamsplitter, when they are in a separable state at the input,
with each photon described by a density matrix $\rho_{\mathrm{I}}$. Notably,
the indistinguishability of the photons manifests as a dip along the
degenerate $\Omega_{j}=\Omega_{k}$ frequencies as $p_{j=k}=0$. This implies
that the probability of measuring any quantities is zero in this case. An
experimental measurement of $p_{jk}$ is shown in Fig.3, showing the bimodal
structure of this distribution. It also shows our labelling convention for
$\Omega_{j(k)}$ where we set $j(k)=0$ to correspond to $\Omega_{0}$ which is
the center frequency of the down-converted light.

Figure 6: Left: simulation of the heralded JSI $F_{jk}$ for pure states, from
Eq.(39), using experimental parameters derived from the sources’ JSI
measurement. Right: experimental result from Merkouche  _et al._ (2021)
obtained by acquiring $30000$ spectral coincidences over 10 hours and binning
the JSI into frequency bins, labelled according to Fig.3.

Next we compute the JSI $F_{jk}$ associated with the state $\ket{\Psi_{jk}}$
from (20), which yields:

$\displaystyle
F_{jk}(\omega_{1},\omega_{2})=\frac{\left|\phi_{j}(\omega_{1})\phi_{k}(\omega_{2})-\phi_{j}(\omega_{2})\phi_{k}(\omega_{1})\right|^{2}}{2\mathcal{C}_{jk}},$
(39)

which can be approximated for $\left|\Omega_{j}-\Omega_{k}\right|\gg 0$, i.e.
for distant heralding frequencies, as

$\displaystyle
F_{jk}(\omega_{1},\omega_{2})\simeq\frac{\left|\phi_{j}(\omega_{1})\phi_{k}(\omega_{2})\right|^{2}+\left|\phi_{j}(\omega_{2})\phi_{k}(\omega_{1})\right|^{2}}{2},$
(40)

which corresponds to two Gaussian spots centered at $(\omega_{j},\omega_{k})$
and $(\omega_{k},\omega_{j})$, hence mirror symmetric about the
$\omega_{1}=\omega_{2}$ axis. For degenerate heralding events
$\Omega_{j}=\Omega_{k}$, the heralded JSI is zero as the probability of
performing a BSM in that case is null.

In the intermediate case where $\Omega_{j}\simeq\Omega_{k}$, the Gaussian
spots are close to overlapped but the cross terms obtained by expanding (39)
maintains the symmetry of the heralded Bell state. These contribute to the JSI
as a dip along the degenerate $\omega_{1}=\omega_{2}$ frequencies, akin to
Eq.(38). The resolved heralded JSI is simulated in Fig.6a) for the pure state
approximation, showing the aforementioned behaviour. The retrieve the same
quantity experimentally, we set up the experiment for the characterization
procedure (Fig.2c) and measure spectral coincidences at the output of the
signal FBS heralded by a BSM at frequencies $\Omega_{j}$,$\Omega_{k}$
(Fig.2b). The result is shown in Fig.6b) where the frequency bins are labelled
according to Fig.3. It is in good agreement with the pure state simulation,
although the degenerate diagonal $j=k$ is not perfectly zero. This is due to
the spectral resolution of the TOFS, where the degenerate frequency bins
contain events corresponding to $\Omega_{j}\simeq\Omega_{k}$. Note also that
both color maps are normalized, such that the amount of energy per bin has to
be multiplied by the probability of realizing this measurement, given by
$p_{jk}$ in Fig.3. Therefore, experimentally, the amount of counts per bin
exactly follows Fig.3. In fact, the $p_{jk}$ distribution can be obtained by
computing the total count of events detected per bin in Fig.6b).

Tracing over the heralding frequencies, we obtain the mixed-state JSI
$F(\omega_{1},\omega_{2})$ from Eq.(22), which simplifies to

$\displaystyle
F(\omega_{1},\omega_{2})=\frac{1}{2}\left[\rho_{\mathrm{S}}(\omega_{1},\omega_{1})\rho_{\mathrm{S}}(\omega_{2},\omega_{2})-\left|\rho_{\mathrm{S}}(\omega_{1},\omega_{2})\right|^{2}\right],$
(41)

where we recall that we have multiplied each element $F_{jk}$ from Eq.(39) by
the proper weight $p_{jk}$ according to Eq.(21). First we note that this is
analogous to the idler JSI from Eq.(38). Indeed, this is the joint spectral
distribution that would be obtained were the beamsplitter placed in the signal
paths rather than the idler paths, and the fact that such a distribution is
measured without the presence of a beamsplitter is evidence of the non-local
nature of this fourfold measurement. In Fig.7a), we simulate the heralded JSI
with our approximated model, showing again the binodal structure
characteristic of a Bell state.

With summing our experimental data from Fig.6b over all heralded bins, we
obtain the histogram shown in Fig.7b), which closely matches the simulation.
We note that the JSI is indeed zero for degenerate frequencies, which is a
consequence of photon-bunching. The measurements of the heralded JSI show the
validity of the pure state approximation for the resolution of the heralding
TOFS, which is sufficiently narrow.

We now proceed to evaluate the entanglement verification signal
$P_{jk}(\tau_{\text{S}})$ from Eq.(28), for this case of $\theta_{jk}=0$. We
find that $\mathcal{F}_{j}(\tau_{\text{S}})=\mathcal{F}_{k}(\tau_{\text{S}})$
since the Gaussians $\phi_{j(k)}$ only differ in their first moment
$\omega_{j(k)}$ (see Eqs.(35) and (36)), such that we have:

$\displaystyle P_{jk}(\tau)=$
$\displaystyle\frac{1}{2\left(1-|\braket{\phi_{j}}{\phi_{k}}|^{2}\right)}\times$
$\displaystyle\Big{\\{}1+e^{-\sigma_{\text{S}}^{2}\tau_{\text{S}}^{2}}\cos{\big{[}(\omega_{j}-\omega_{k})\tau_{\text{S}}\big{]}}$
$\displaystyle-|\braket{\phi_{j}}{\phi_{k}}|^{2}(1+e^{-\sigma_{\text{S}}^{2}\tau_{\text{S}}^{2}})\Big{\\}}.$
(42)

Similar to the heralded JSI, we analyzing this function depending on the
spectral distance between $\Omega_{j}$ and $\Omega_{k}$. We find that for
$\left|\Omega_{j}-\Omega_{k}\right|\gg 0$,
$\braket{\phi_{j}}{\phi_{k}}\rightarrow 0$ such that the verification signal
reduces to

$\displaystyle
P_{jk}(\tau_{\text{S}})=\frac{1+e^{-\sigma_{\text{S}}^{2}\tau_{\text{S}}^{2}}\cos{\big{[}(\omega_{j}-\omega_{k})\tau_{\text{S}}\big{]}}}{2}.$
(43)

It corresponds to a background probability of $1/2$ and to a Gaussian envelope
with a that corresponds to the transform limited temporal width of the
signal’s photon. This envelope is modulated by fringes at the difference of
the heralded frequencies $\omega_{j}-\omega_{k}$. We note that the fringes
vanish when setting $\alpha=0$ in (36), corroborating that the observation of
fringes is a witness of entanglement swapping in the considered frequency
bins.

Since in the pure state model, the probability $p_{jk}$ of obtaining a
coincidence at degenerate frequencies $\Omega_{j}=\Omega_{k}$ is null, the
verification signal (42) is not defined so we set this $P_{j=k}=0$ for the
simulation. However, we find that it has the following limit:

$\displaystyle P_{j\rightarrow
k}(\tau_{\text{S}})\rightarrow\frac{1}{2}-\frac{1}{2}\left(2\sigma_{\text{S}}^{2}\tau_{\text{S}}^{2}-1\right)e^{-\sigma_{\text{S}}^{2}\tau_{\text{S}}^{2}},$
(44)

which can be easily demonstrated by noticing that
$\braket{\phi_{j}}{\phi_{k}}$ is a Gaussian of the variable
$\omega_{j}-\omega_{k}$ under the Gaussian approximation (see for instance
Eq.(37)). The expression is similar to that reported in Graffitti  _et al._
(2020) that utilizes an engineered non-linearity to obtain spectral Bell
states. Interestingly, the full signal $P(\tau_{\text{S}})$ that we describe
next in Eq.(45) has a similar expression with more complicated dependency on
the experimental parameters and both describe the Gaussian peak sitting in a
HOM dip. The former is a proof of entanglement whilst the latter results from
HOM interferences for close-to-degenerate heralded frequencies. On Fig.8a), we
plotted the simulated $P_{jk}$ using our approximated model and parameters
obtained experimentally. This plot shows the previously described behaviour,
showing oscillations at the difference frequency which merge into a single
peak in the near degenerate case. The background colormap for these plots
represent the probability of measuring these events, i.e. $p_{jk}$ from
Eq.(38).

Figure 7: Left: simulation of the full heralded JSI $F$ defined by the sum of
Eq.(21) over all $j,k$. Right: experimental result obtained by summing the
acquisitions from Fig.6b over all bins. Figure 8: Left: simulation of the
entanglement verification signal $P_{jk}$ for pure states, see Eq.(42). Right:
experimental result from Merkouche  _et al._ (2021), obtained by binning the
verification signal acquired over 15 hours. These results are fitted to the
approximated model from (43). The frequency bins are labelled according to
Fig. 3. The colormap in the background represents the probability $p_{jk}$ of
an heralding event, theoretical (top) and experimental (bottom).

In Fig.8b), we show the experimental counterpart, with the colored background
representing the experimental $p_{jk}$ obtained by computing the total number
of counts in each bin. Again we see a nice agreement with the simulation
except for the extreme bins which have a higher oscillating frequency in the
simulated case. This is again the pure state approximation that results in
$\phi_{j(k)}$ modes that are exactly centered at $\omega_{j(k)}$ whereas in
the experimental case, some averaging occurs due to the finite resolution of
the TOFS. This is illustrated by Fig.5 where the heralding and heralded center
frequencies $\omega_{j}$ and $\Omega_{j}$ can be different depending on the
resolution of the spectral detection. Nevertheless, this shows that our
approximated model, and most notably the limit derived in Eq.(43), is
sufficient to describe the entanglement swapping verification protocol, and we
therefore use it to fit our results, shown as solid red curve in the
experimental plot.

Finally we evaluate the verification signal $P(\tau_{\text{S}})$ with
$\theta_{jk}=0$ for the non-spectrally-resolved case, given by Eq.(30). We
find that it contains four terms:

$\displaystyle P(\tau_{\text{S}})=\frac{1}{4}\Bigg{[}$ $\displaystyle
1+\left|\int\textrm{d}\omega\textrm{d}\Omega\,\left|f(\omega,\Omega)\right|^{2}e^{i\omega\tau_{\text{S}}}\right|^{2}$
$\displaystyle-\int\textrm{d}^{2}\Omega\,\left|\rho_{\mathrm{I}}(\Omega,\Omega^{\prime})\right|^{2}$
$\displaystyle-\int\textrm{d}^{2}\omega\,\left|\rho_{\mathrm{S}}(\omega,\omega^{\prime})\right|^{2}e^{i(\omega-\omega^{\prime})\tau_{\text{S}}}\Bigg{]}.$
(45)

The first one is simply a background probability, while the second one
corresponds to the overlap integral between the two sources with a relative
delay $\tau_{\text{S}}$ between the signal photons. Evaluating this term
reveals a Gaussian along $\tau_{\text{S}}$ whose width depends on the joint
temporal distribution of the sources. This is quite similar to the cross-
correlation between two classical pulses, except that in the present case, the
phase of the fringes is constant, implying that the $P_{jk}$ sum coherently to
a single peak at $\tau_{\text{S}}=0$. We will see in the next section how that
phase can be offset by introducing an additional time delay in the BSM.

The last two terms correspond respectively to the overlap integrals between
the idler and the signal density matrices of each source. The former evaluate
to a constant (which is unity when the sources are perfectly matched), while
the latter describes an unheralded HOM dip between the signals photons. Hence,
the full verification signal can be summarized as a Gaussian peak centered in
a HOM dip. In Fig. 9a), we plotted a simulation of the full signal
$P(\tau_{\text{S}})$ showing this behaviour. In Fig. 9b), we show the
experimental result obtained without binning that data from Fig.8, which shows
good agreement with the approximated model.

We stress that the presence of oscillating fringes in $P_{jk}$ or a peak in
$P$, where the coincidence probability goes above the baseline of $1/2$, is a
witness of an entangled state (see, for example, Fedrizzi  _et al._ (2009)).
Our setup is therefore capable of performing entanglement swapping between a
large amount of frequency Bell states. As we show in Merkouche  _et al._
(2021), not all of the heralded states shown in Fig.6 are mutually orthogonal,
but it is possible to select multiple subsets which form a set of mutually
orthogonal Bell states. This is described in Appendix E, where we run an
algorithm to select 6 sets of 5 orthogonal modes, all of which satisfy the
verification procedure for entanglement swapping.

Figure 9: Probability $P(\tau)$ of a coincidence heralded by a BSM without
spectral resolution. Left: simulation from Eq.(45); Right: experimental result
presented in Merkouche  _et al._ (2021) obtained by summing the individual
$P_{jk}$ that are depicted in Fig.8. The solid red curve is a sum of the
individual fits on the experimental Gaussian model from Eq.(43).

###  IV.3 Case of $\tau_{\text{I}}\neq 0$ ($\theta_{jk}\neq 0$)

We now study the case when there is a temporal delay $\tau_{\text{I}}$ between
the idler photons. The resulting HOM interferences between them therefore have
a decreased visibility. Note that when this delay becomes too large, the idler
photons become distinguishable at the output of the beamsplitter, effectively
removing the entanglement in the heralded Bell state. The formalism that we
choose for our model necessitates to start from an entangled state. Therefore,
in the scope of the paper, we will only consider a delay $\tau_{\text{I}}$
that is smaller than the inverse bandwidth $\sigma_{\text{I}}^{-1}\approx 300$
fs of the idler photons such that the state $\ket{\Psi_{jk}}$ retains some
entanglement. One limitation of the pure state model arises because the idlers
BSM (6) is achieved with perfect resolution, leading to infinite temporal
support of the quantities defined in Sec.II over $\tau_{\text{I}}$ whereas
they should necessarily be bounded by an envelope inversely proportional to
the spectral resolution. In Fig.5, we see how the width of a realistic
filtering of the idler photon differs from an infinitely narrow spectral
detection. The derivation for the mixed state case is shown in Appendix A.
Nevertheless, the pure state model predict the correct physical behaviour and
is therefore sufficient to describe the effects of slight distinguishability
in the BSM.

The heralded state in the pure state approximation is given by Eq.(12) and
shows a dephasing between the two states that depends on $\tau_{\text{I}}$.
The heralded JSI from Eq.(20) evaluates to:

$\displaystyle
F_{jk}(\omega_{1},\omega_{2})=\frac{1}{2\mathcal{C}_{jk}}\Big{[}$
$\displaystyle\left|\phi_{j}(\omega_{1})\phi_{k}(\omega_{2})\right|^{2}+\left|\phi_{k}(\omega_{1})\phi_{j}(\omega_{2})\right|^{2}$
$\displaystyle-2\gamma_{jk}(\omega_{1},\omega_{2})\cos{(\Omega_{j}-\Omega_{k})\tau_{\text{I}}}\Big{]},$
(46)

where
$\gamma(\omega_{1},\omega_{2})=\phi_{j}(\omega_{1})\phi_{k}(\omega_{2})\phi_{j}(\omega_{2})\phi_{k}(\omega_{1})$
for real modes $\phi_{j(k)}$, for simplicity. The last term is responsible for
the HOM dip along the degenerate frequencies $\omega_{1}=\omega_{2}$, marking
indistinguishability. As previously, $\gamma_{jk}\rightarrow 0$ when heralding
distant bins $\left|\Omega_{j}-\Omega_{k}\right|\gg 0$ so the delay between
the idlers photons has no influence. However, for degenerate bins
$\Omega_{j}\rightarrow\Omega_{k}$ and for
$\tau_{\text{I}}<\sigma_{\text{I}}^{-1}$,
$\gamma_{jk}\times\cos\theta_{jk}\rightarrow 0$ which causes the two Gaussian
spots to merge.

Putting these limits together, we find that a small delay between the idlers
results in a heralded JSI similar to Fig.6a) where the JSI in the bins close
to the $j\rightarrow k$ diagonal are more or less merged depending on
$\tau_{\text{I}}$. For spectral bins that are spaced further, the JSI is
unchanged since the distinguishability is marked by a relative phase between
the Gaussian spots.

The full JSI is obtained by summing $F_{jk}$ according to Eq.(21):

$\displaystyle
F(\omega_{1},\omega_{2})=\frac{\rho_{\mathrm{S}}(\omega_{1},\omega_{1})\rho_{\mathrm{S}}(\omega_{2},\omega_{2})-\Gamma(\omega_{1},\omega_{2};\tau_{\text{I}})}{2},$
(47)

where $\Gamma$ is a function that depends on the overlap of the signal density
matrices as a function of the idlers’ delay (see Appendix F). Similar to the
previous case, the overlap between $\phi_{j}$ and $\phi_{k}$ depends on the
indistinguishably in time of the idlers.

We next consider the entanglement verification signal from Eq.(28), which
depends on both $\tau_{\text{S}}$ and $\tau_{\text{I}}$ through $\theta_{jk}$.
Evaluating it with the Gaussian model, we then obtain:

$\displaystyle P_{jk}$
$\displaystyle(\tau_{\text{S}},\tau_{\text{I}})=\frac{1}{2(1-|\braket{\phi_{j}}{\phi_{k}}|^{2}\cos{\theta_{jk}})}\times$
$\displaystyle\Bigg{(}1+e^{-\sigma_{\text{S}}^{2}\tau_{\text{S}}^{2}}\cos{\big{[}\Delta\omega_{jk}\tau_{\text{S}}-\theta_{jk}\big{]}}$
$\displaystyle-|\braket{\phi_{j}}{\phi_{k}}|^{2}(e^{-\sigma_{\text{S}}^{2}\tau_{\text{S}}^{2}}+\cos{\theta_{jk}})\Bigg{)}.$
(48)

We perform yet another asymptotic behaviour analysis. For distant heralding
frequencies, we find that the limit is similar to Eq.(43) with an additional
dephasing:

$\displaystyle
P_{jk}(\tau_{\text{S}},\tau_{\text{I}})=\frac{1+e^{-\sigma_{\text{S}}^{2}\tau_{\text{S}}^{2}}\cos{\big{[}(\omega_{j}-\omega_{k})(\tau_{\text{S}}-\tau_{\text{I}}^{\prime})\big{]}}}{2},$
(49)

where we used (36) to factorize by the difference of heralded frequencies and
we defined
$\tau_{\text{I}}^{\prime}=\tau_{\text{I}}/2\alpha\sigma_{\text{S}}^{2}$. We
can see that the value of $\tau_{\text{I}}$ has a more noticeable effect
compared to the heralded JSI. The fringes are no longer synchronized to the
envelope and a phase shift occurs when the delay between the idlers photons is
nonzero. This effect is quite important in our case since a delay as small as
$100$ fs between the idler photons is sufficient to cause a dephasing of
$\pi$, due to relatively large spectral bandwidth of the JSA.

Figure 10: Simulated (left) and experimental (right) probability to get a
coincidence heralded by a BSM at frequencies $\Omega_{j}$ and $\Omega_{k}$.
The simulation utilizes the Gaussian, pure state approximation from Eq.(49)
with the same parameters as in the experimental case. The experimental plots
are acquired over 900s for different positions of the idler stage (on the
right). The frequency bins used are separated by 8nm (or about 10 THz). The
red curve represents a fit to the theoretical model. This is a different
representation of the data shown in Fig.11 in the bin labelled $(2,-2)$.

Note that using the more realistic mixed state model, the envelope of the
fringes is also affected by this delay and causes a reduction in visibility.
In the pure state approximation, the envelope in the idler’s direction is
infinite, but it is sufficient to show the most noticeable effect of the
dephasing. Note also that it is possible to find a value for both delays such
that $\tau_{\text{S}}=\tau_{\text{I}}^{\prime}$, in which case there is no
dephasing, but the visibility of the fringes would be decreased. This is
illustrated in Fig.10 which compares a simulation (left) that uses the
approximated model with the same parameters as the experimental results
(right). The entanglement verification utilizes the previous experimental
protocol for 7 different values of $\tau_{\text{I}}$. While we see again very
good agreement between the experimental and simulated fringes, the envelope in
the model has no dependency on $\tau_{\text{I}}$ whereas it is clearly the
case experimentally.

For near-degenerate heralding frequencies, we have the following limit:

$\displaystyle P_{j\rightarrow
k}(\tau_{\text{S}},\tau_{\text{I}})=\frac{1}{2}-\frac{1}{2}\cdot\frac{2\sigma_{\text{S}}^{2}(\tau_{\text{S}}-\tau_{\text{I}}^{\prime})^{2}-1}{1+4(\tau_{\text{I}}^{\prime})^{2}}e^{-\sigma_{\text{S}}^{2}\tau_{\text{S}}^{2}}.$
(50)

which is equal to Eq.(44) when setting $\tau_{\text{I}}=0$. This represents
again a Gaussian peak centered in a HOM dip at $\tau_{\text{S}}=0$. A change
of variable delay between the idler photons then decreases its visibility.

Putting both limits together, we find that a scan over $\tau_{\text{S}}$ and
$\tau_{\text{I}}$ of the verification signal $P_{jk}$ from Eq.(48) look very
similar to Fig.8a), except that the fringes will be offset as a function of
$\tau_{\text{I}}$ while the peak close to the diagonal remains centered. In
both cases, the visibility is decreased. This behaviour was very useful
experimentally to verify that the delay at both FBS was as close to zero as
possible. Moreover, the decrease in contrast for non-zero values of
$\tau_{\text{I}}$ is yet another proof that the visibility of the oscillations
in $P_{jk}$ is a marker of entanglement.

Figure 11: Waterfall plot of $P_{jk}$ plotted against $\tau_{\text{S}}$
(horizontal) for different values of $\tau_{\text{I}}$ (vertical) where each
plot represents a frequency bin $\Omega_{j},\Omega_{k}$ labelled by $j$ and
$k$. The plot range and units for $\tau_{\text{S}}$ and $\tau_{\text{I}}$ are
the same as shown in Fig. 10

.

In Fig.11 we show an experimental waterfall plot representing each $P_{jk}$
for values of $\tau_{\text{I}}$ ranging from -300 to +300 fs which originates
from the same dataset as Fig.10. Note that the spectrally resolved heralding
is done at half the resolution than used earlier (shown in Fig.3) to have
sufficient statistics. These waterfall plots confirm the expected behaviour as
we can see how the fringes are slanted in the
$(\tau_{\text{S}},\tau_{\text{I}})$ space for distant $j,k$ frequency bins as
predicted by (49), while in the degenerate $j=k$ case, the fringes collapse to
a single peak as described by Eq.(50). Note that the angle of the fringes is
defined by the proportionality factor between $\tau_{\text{I}}$ and
$\tau_{\text{I}}^{\prime}$, which depends on the amount of entanglement and on
the spectral bandwidth of the signal photons. These structures show that
entanglement swapping is still achieved for this specific range of delay
mismatch in the BSM, inducing distinguishability.

Finally, by repeating the same experiment either without spectral resolution
of the herald or by summing the individual $P_{jk}$ according to Eq.(30), we
find that the expression of the total verification probability is given by:

Figure 12: Top: Simulation of (51) with experimental parameters. Bottom:
experimental acquisition obtained by monitoring four-fold coincidences over 15
hours while scanning both $\tau_{\text{S}}$ and $\tau_{\text{I}}$ without
spectral resolution of the herald, fitted to the function
$P(\tau_{\text{S}},\tau_{\text{I}})$.

$\displaystyle P(\tau_{\text{S}},\tau_{\text{I}})=$
$\displaystyle\frac{1}{4}\Bigg{(}1+\left|\int\textrm{d}\omega\textrm{d}\Omega\,\left|f(\omega,\Omega)\right|^{2}e^{i(\omega\tau_{\text{S}}+\Omega\tau_{\text{I}})}\right|^{2}$
$\displaystyle-\int\textrm{d}^{2}\Omega\,\left|\rho_{\mathrm{I}}(\Omega,\Omega^{\prime})\right|^{2}e^{i(\Omega-\Omega^{\prime})\tau_{\text{I}}}$
$\displaystyle-\int\textrm{d}^{2}\omega\,\left|\rho_{\mathrm{S}}(\Omega,\Omega^{\prime})\right|^{2}e^{i(\omega-\omega^{\prime})\tau_{\text{S}}}\Bigg{)},$
(51)

which is similar to Eq.(45) with an additional dependency on
$\tau_{\text{I}}$. We simulated this expression in Fig.12 and show an
experimental acquisition as well. It is straightforward to identify the three
nonconstant terms as familiar quantities. The second term is the cross-
correlation between both JSA as a function of both delays. It can be written
as the product of Fourier transform of the joint spectrum with respect to
$\omega$ and $\Omega$, thus reducing to the product of envelope functions
centered at $\tau_{\text{S}}$ and $\tau_{\text{I}}$. This term is responsible
for the slanted peak in Fig.12. The last two terms correspond respectively to
the overlap integrals between the density matrices of the idlers and of the
signals. They show the effect of interferences between uncorrelated photons
and hence are visible as HOM dips along either $\tau_{\text{S}}$ or
$\tau_{\text{I}}$, as shown again in Fig.12. This figure allows us to identify
more clearly the range over which entanglement swapping can be verified, which
is essentially the area over which the slanted peak appears since the
structure outside corresponds to quantum interferences between unentangled
single photon states.

By exploring the range of both time delays in this configuration, we
explicitly see the ultimate equivalence between signals and the idlers in what
is effectively a four-photon interferometer. As such, it is adequate to
consider $P_{jk}(\tau_{\text{S}},\tau_{\text{I}})$ and
$P(\tau_{\text{S}},\tau_{\text{I}})$ as the most general representation of the
phase-insensitive part of four-photon interference where time-frequency
entanglement is present.

##  V Conclusion

In this work, we have undertaken a more thorough and more general analysis of
the results we report on in Ref. Merkouche  _et al._ (2021). In summary, we
demonstrated and analyzed a novel scheme for time-frequency entanglement
swapping, using a multimode, spectrally-resolved Bell-state measurement as the
heralding mechanism. The most salient feature of our method is the heralding
of several, mutually-orthogonal Bell states derived from identical multimode
entangled photon pairs. We further generalized our result to consider the case
of non-zero time delay in the heralding Bell-state measurement, giving rise to
Bell states with a varying amount of phase. Our setup is the first known to
the authors to incorporate four simultaneous time-of-flight spectrometers, and
thus points towards a promising venue of study of time-frequency entangled
four-photon interferometry. .

## Acknowledgements

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under Grant Agreement No. 665148, the United
Kingdom Defense Science and Technology Laboratory (DSTL) under contract No.
DSTLX-100092545, and the National Science Foundation under Grant No. 1620822.

## Appendix

##  Appendix A Mixed state model

###  A.1 Heralded state and JSI

In the realistic case, the idlers BSM is not performed with perfect
resolution, but rather with a finite spectral window. In our case, this is due
to the resolution of the time-of-flight spectrometer, which is a convolution
of multiple response functions in the frequency-to-time conversion. It is
dominated by the timing jitter ($\simeq 20$ ps) of the superconducting
nanowires.When this resolution is not perfect, then we can show that the
signal photons are heralded into a mixed state.

We begin by redefining the BSM operator as:

$\begin{gathered}\hat{\Pi}_{lm}^{\mathrm{BSM}}=\int\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}|t_{l}(\Omega_{j})|^{2}|t_{m}(\Omega_{k})|^{2}\
\hat{\Pi}^{\text{BSM}}_{jk}\\\
=\int\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}|t_{l}(\Omega_{j})|^{2}|t_{m}(\Omega_{k})|^{2}\ket{\Omega_{j},\Omega_{k}}\bra{\Omega_{j},\Omega_{k}},\\\
\ket{\Omega_{j},\Omega_{k}}=\hat{c}^{\dagger}(\Omega_{j})\hat{d}^{\dagger}(\Omega_{k})\ket{\text{vac}},\end{gathered}$
(52)

where $t_{l(m)}(\Omega_{j(k)})$ are transmission amplitudes centered at
$\Omega_{l(m)}$, and satisfy
$\sum_{lm}|t_{l}(\Omega_{j})|^{2}|t_{m}(\Omega_{k})|^{2}=1$. We have
introduced new indices $l$, $m$, so that we may incorporate the
$j,k$-dependent quantities from the main text into our analysis. Moreover, it
is straightforward to show that the POVM (6) is obtained by setting
$t_{l(m)}(\Omega)\rightarrow\delta(\Omega-\Omega_{l(m)})$. To simplify our
notation henceforth, we will use the following shorthand

$\iint_{lm}\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}:=\int\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}|t_{l}(\Omega_{j})|^{2}|t_{m}(\Omega_{k})|^{2}.$
(53)

It can be seen from Fig.5 that the finite resolution of the idler detection
bin reduces the purity of the heralded signal state, due to the entanglement
of the two-photon state. Because of this, the POVM element (52) requires that
we describe the heralded state as a mixed state

$\displaystyle\hat{\rho}_{lm}=\frac{\operatorname{Tr}_{\hat{b}}\Big{[}\hat{\Pi}_{lm}^{\mathrm{BSM}}\ket{\psi_{12}}\bra{\psi_{12}}\Big{]}}{p_{lm}},$
(54)

where $\operatorname{Tr}_{\hat{b}}$ is the partial trace over the subspace
defined by operators $\hat{b}_{1}$ and $\hat{b}_{2}$. Analogously to the pure
state case, the probability $p_{lm}$ are defined as

$\displaystyle p_{lm}$
$\displaystyle=\operatorname{Tr}\Big{[}\hat{\Pi}_{lm}^{\mathrm{BSM}}\ket{\psi_{12}}\bra{\psi_{12}}\Big{]}$
$\displaystyle=\frac{1}{2}\iint_{lm}\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}\Big{(}\rho_{\mathrm{I}}(\Omega_{j},\Omega_{j})\rho_{\mathrm{I}}(\Omega_{k},\Omega_{k})$
$\displaystyle\quad-\rho_{\mathrm{I}}(\Omega_{j},\Omega_{k})\rho_{\mathrm{I}}(\Omega_{k},\Omega_{j})\Big{)}$
$\displaystyle=\iint_{lm}\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}\ p_{jk}.$
(55)

where the idler density matrix is defined as Eq.(9), and we again obtain
Eq.(10) by setting the filters $t_{l(m)}$ as $\delta$ functions.

We may now compute the heralded state density matrix:

$\hat{\rho}_{lm}=\frac{1}{p_{lm}}\iint_{lm}\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}\
p_{jk}\ \ket{\Psi_{jk}}\bra{\Psi_{jk}}$ (56)

while recalling the definition of $\phi_{j(k)}(\omega)$ in $\ket{\Psi_{jk}}$
as

$\phi_{j(k)}(\omega)=\frac{f(\omega,\Omega_{j(k)})}{\sqrt{\rho_{\mathrm{I}}(\Omega_{j(k)},\Omega_{j(k)})}}.$
(57)

This has the intuitive interpretation of a mixed state as an incoherent sum of
pure states over the detection bandwidth of $t_{l(m)}$. All of the measured
quantities follow in a straightforward manner. In the absence of frequency
resolution, we herald again the mixed state

$\hat{\rho}=\sum_{lm}\ p_{lm}\ \hat{\rho}_{lm},$ (58)

just as with the pure state model.

The heralded JSI is given by:

$\displaystyle
F_{lm}(\omega_{1},\omega_{2})=\bra{\omega_{1},\omega_{2}}\hat{\rho}_{lm}\ket{\omega_{1},\omega_{2}},$
(59)

which can be expressed in terms of $F_{jk}$ as

$\displaystyle
F_{lm}(\omega_{1},\omega_{2})=\frac{1}{p_{lm}}\iint_{lm}\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}\
p_{jk}\ F_{jk}(\omega_{1},\omega_{2}).$ (60)

The integrated JSI corresponding to the state $\hat{\rho}$ is again given by

$\displaystyle
F(\omega_{1},\omega_{2})=\sum_{lm}p_{lm}F_{lm}(\omega_{1},\omega_{2}).$ (61)

When the signal photons in the state $\hat{\rho}_{lm}$ are incident on a 50:50
beamsplitter, the coincidence fringes at the output are given by

$\displaystyle
P_{lm}(\tau_{\text{S}},\tau_{\text{I}})=\text{Tr}\left(\hat{\Pi}_{\text{verif}}\,\hat{\rho}_{lm}\right)=\int\textrm{d}^{2}\omega\bra{\omega,\omega^{\prime}}\hat{\rho}_{lm}\ket{\omega,\omega^{\prime}},$
(62)

where
$\ket{\omega,\omega^{\prime}}=\hat{x}^{\dagger}(\omega)\hat{y}^{\dagger}(\omega^{\prime})\ket{\text{vac}}$
as before. Evaluating this in terms of $P_{jk}$, we

$\displaystyle
P_{lm}(\tau_{\text{S}},\tau_{\text{I}})=\frac{1}{p_{lm}}\iint_{lm}\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}\
p_{jk}\ P_{jk}(\tau_{\text{S}},\tau_{\text{I}}).$ (63)

Finally the integrated interference peak for the state $\hat{\rho}$ is
recovered by taking

$P(\tau_{\text{S}},\tau_{\text{I}})=\sum_{lm}\ p_{lm}\
P_{lm}(\tau_{\text{S}},\tau_{\text{I}}).$ (64)

We will make a few comments regarding the most interesting feature of
comparing this model with the pure state model. To this end, we shall
represent the pure state $\ket{\Psi_{jk}}\bra{\Psi_{jk}}$ as a density matrix
in the $\\{\ket{\phi_{j}}\ket{\phi_{k}},\ket{\phi_{k}}\ket{\phi_{j}}\\}$ basis
as follows:

$\ket{\Psi_{jk}}\bra{\Psi_{jk}}=\frac{1}{2\mathcal{C}_{jk}}\begin{pmatrix}1&e^{i\theta_{jk}}\\\
-e^{-i\theta_{jk}}&1\end{pmatrix}.$ (65)

Meanwhile, the mixed state $\hat{\rho}_{lm}$ has the representation:

$\hat{\rho}_{lm}=\frac{1}{p_{lm}}\iint_{lm}\textrm{d}\Omega_{j}\textrm{d}\Omega_{k}\
\frac{p_{jk}}{2\mathcal{C}_{jk}}\ \begin{pmatrix}1&e^{i\theta_{jk}}\\\
-e^{-i\theta_{jk}}&1\end{pmatrix}.$ (66)

Note that the off-diagonal terms $e^{\pm i\theta_{jk}}=e^{\pm
i(\Omega_{j}-\Omega_{k})\tau_{\text{I}}}$, the coherence terms, are the
hallmark of the bipartite entanglement in these states, and are responsible
for the interference we observe in $P_{jk}$. It seems reasonable then to ask
if, and to what extent, the averaging over the $l,m$ bands in the mixed state
is expected to reduce the entanglement. Take first the special case of
$\tau_{\text{I}}=0$. In this case, $\theta_{jk}=0\ \forall\ j,k$, and the off-
diagonal terms are equal to unity, and hence, no phase-averaging occurs for
the mixed state $\hat{\rho}_{lm}$. Indeed, the off-diagonals survive even for
the fully-averaged state $\hat{\rho}$, and one can interpret this as the
reason why the interference peak survives at full visibility for
$\tau_{\text{I}}=0$.

More generally, however, for $\tau_{\text{I}}\neq 0$, the disagreement between
the two models becomes more salient. The pure-state model assumes that the
idler photons are detected at monochromatic frequency modes at
$(\Omega_{j},\Omega_{k})$, which are necessarily infinite in extent in the
time domain. This means that the idler photons remain indistinguishable at the
output of the BSM beamsplitter, even for arbitrary $\tau_{\text{I}}$ delays at
the input. On the other hand, taking into account the averaging over the $l,m$
bands introduces distinguishability, and the BSM is no longer ideal in this
case. Furthermore, the phase-averaging becomes worse the larger
$\tau_{\text{I}}$ is, since the argument of the phase factor is
$(\Omega_{j}-\Omega_{k})\tau_{\text{I}}$. This feature results in a reduction
of the off-diagonal terms, and of the visibility of the $P_{lm}$ interference
for large enough $\tau_{\text{I}}$. This behavior is most clearly seen in the
profile of the fully-integrated two-dimensional peak, which results from
averaging over all phases with a weight $p_{jk}$. The peak vanishes in the
$\tau_{\text{I}}$ direction over a delay timescale comparable to the inverse
bandwidth of the idler photons. The intuitive interpretation is that the idler
photons become distinguishable at the output of the beamsplitter when the
relative delay is at least as long as their pulse durations.

##  Appendix B Source distinguishability

The entanglement verification protocol we use, that is, the two-photon
interference of the state $\ket{\Psi_{jk}}$, ultimately relies on the
indistinguishability of the two source states. To see this, we relabel the
source JSA’s as $f_{1}(\omega,\Omega)$ and $f_{2}(\omega,\Omega)$, and for
simplicity, we assume that they are identical up to a translation in frequency
space. Note now that this leads to a heralded state

$\displaystyle\ket{\Psi_{jk}}\propto\ket{\phi^{1}_{j}}_{1}\ket{\phi^{2}_{k}}_{2}-\ket{\phi^{1}_{k}}_{1}\ket{\phi^{2}_{j}}_{2},$
(67)

where

$\displaystyle\ket{\phi^{1(2)}_{j(k)}}=\int\textrm{d}\omega\phi^{1(2)}_{j(k)}(\omega)\hat{a}^{\dagger}_{1(2)}(\omega)\ket{\text{vac}},$
(68)

and

$\displaystyle\phi^{1(2)}_{j(k)}(\omega)=\frac{f_{1(2)}(\omega,\Omega_{j(k)})}{\rho_{1(2)}(\Omega_{j(k)},\Omega_{j(k)})}.$
(69)

Although this state is still entangled, the verification method using
coincidence fringes in $P_{jk}(\tau)$ will suffer from a reduction in
visibility due to the distinguishability of $f_{1}$ and $f_{2}$. To see this,
we recalculate $P_{jk}(\tau)$ in its approximate form (49), and find

$\displaystyle
P_{jk}(\tau_{S})\approx\frac{1}{2}\left(1-V_{jk}\,e^{\sigma_{\text{S}}^{2}\tau_{S}^{2}}\cos{\left[\left(\frac{\omega^{1}_{j}+\omega^{2}_{j}}{2}-\frac{\omega^{1}_{k}+\omega^{2}_{k}}{2}\right)\tau_{S}\right]}\right)$
(70)

where the visibility $V_{jk}$ is given by

$\displaystyle
V_{jk}=\left(\int\textrm{d}\omega\phi^{1*}_{j}(\omega)\phi^{2}_{j}(\omega)\right)\left(\int\textrm{d}\omega\phi^{1*}_{k}(\omega)\phi^{2}_{k}(\omega)\right).$
(71)

We can maximize this visibility by maximizing the overlap $f_{1}$ and $f_{2}$.
We see that this latter provides a lower bound on $V_{jk}$ by writing

$\displaystyle\int\textrm{d}\omega\phi^{1*}_{j}(\omega)\phi^{2}_{j}(\omega)=\frac{\int\textrm{d}\omega
f_{1}^{*}(\omega,\Omega_{j})f_{2}(\omega,\Omega_{j})}{\sqrt{\rho_{1}(\Omega_{j},\Omega_{j})\rho_{2}(\Omega_{j},\Omega_{j})}}$
$\displaystyle\geq\int\textrm{d}\omega\textrm{d}\Omega
f_{1}^{*}(\omega,\Omega)f_{2}(\omega,\Omega),$ (72)

and likewise for $k$.

It is relatively straightforward to maximize the quantity on the left by
tuning experimental parameters, namely pump wavelength, phasematching angle,
and transverse optical fiber position (due to residual spatial chirp), and
observing two-fold coincidences resulting from first order interference of the
sources. Because both sources are pumped with the same pulse, the two-photon
term of the state is given by

$\displaystyle\ket{\psi}\propto\int\textrm{d}\omega\textrm{d}\Omega\Big{(}f_{1}(\omega,\Omega)\hat{a}^{\dagger}_{1}(\omega)\hat{b}^{\dagger}_{1}(\Omega)+$
$\displaystyle
f_{2}(\omega,\Omega)\hat{a}^{\dagger}_{2}(\omega)\hat{b}^{\dagger}_{2}(\Omega)\Big{)}\ket{\text{vac}}.$
(73)

A straightforward calculation shows that the probability of a two-fold
coincidence between ports $\hat{c}$ (or $\hat{d})$ and $\hat{x}$ (or
$\hat{y}$) is given by

$\displaystyle
P_{cc}=\frac{1}{4}\int\textrm{d}^{2}\omega\Big{|}f_{1}(\omega,\Omega)\pm
f_{2}(\omega,\Omega)\Big{|}^{2}$
$\displaystyle=\frac{1}{2}\Big{(}1\pm\mathrm{Re}\int\textrm{d}^{2}\omega
f_{1}^{*}(\omega,\Omega)f_{2}(\omega,\Omega)\Big{)}$ (74)

Figure 13: Measured coincidence fringes $P_{\textit{cc}(\tau)}$ with a
contrast of $80\%$.

In the following, we will outline additional measurements to quantify the
source indistinguishability. In our case, our dual-pass geometry implies that
we need to match the JSI of both sources, which is achieved when both signals
and idlers from both sources have maximum overlap. We opted for a bulk crystal
source in Type II to enable pumping in both directions while being able to
separate our four photons into different paths. We used a BiBO crystal due to
its relatively high non linearity.

First, we measured the JSI by directing the two daughter photons from either
source into the fiber spools, since their large spectral bandwidth would be
cropped with the CFBG. The JSI from each source is depicted in Fig4. They show
that both sources are nearly indistinguishable; a singular value decomposition
yields a Schmidt number of $K_{1}=2.9\pm.1$ and $K_{2}=2.9\pm.1$. These values
are lower than the theoretical expectation ($K\sim 5$) because of the timing
jitter of our detectors that result in a broader distribution. This was
confirmed by measuring the JSI with the CFBG’s, which have a better resolution
but are limited in range. The correlation width was found to be lower and
therefore the Schmidt number can be expected to be at least $K=4$.

Figure 14: Left: HOM dip between the signal photons heralded by a coincidence
between the idler photons. Right: same measurement but with spectral
resolution of the heralding photons, labelled $j,k$ for
$\Omega_{k},\Omega_{k}$, where index $j,k=0$ corresponds to the center
frequency $\omega_{0}$.

Note that this method is insensitive to any spectral phase difference, such as
dispersion from the pump, since the second pump is slightly more dispersed
than the first due to propagation. This has been shown to increase the
entanglement and the Schmidt number Davis  _et al._ (2020); Ansari  _et al._
(2018). However, this difference should be negligible, and the method
presented that relies on Eq.(74) allows for a more accurate estimation of the
overlap. Nevertheless, the JSI measurement showed near-perfect correspondence
between the intensity of the two sources which is a critical step to ensure
indistinguishably between the uncorrelated photon pairs.

To further characterize the indistinguishability of the sources, we measure
their heralded $g^{(2)}$ by splitting their signal photon into a beamsplitter.
This yields a value of $g_{1}^{(2)}=0.16\pm 0.003$ and $g_{2}^{(2)}=0.14\pm
0.003$. These values are consistent with the relatively high optical power
that is utilized to pump the sources in order to maximize the probability of
four-fold coincidences. The lower value of $g^{(2)}$ for source 2 is
consistent with the fact that it also has a higher heralding efficiency than
source 1. The reason is not entirely clear, but it is likely that the previous
interaction with the PDC crystal on the first pass results in an additional
filtering on the pump as well as a slight reduction in optical power because
of absorption.

Finally, in Fig13 we measured the coincidences between ports $\hat{c}$ and
$\hat{x}$ (see Fig1) while scanning the relative phase between the two pump
fields with a piezoelectric stack, which is related to the probability from
Eq.(74). We scanned using a slow voltage ramp resulting in a few micrometers
of displacement over a few seconds. The visibility of those fringes is $80\%$,
which is a direct measurement of the overlap between the two sources, and
therefore a quantification of distinguishability.

Note that we also performed this measurement with spectral resolution,
essentially measuring those interferences in narrower spectral bins using our
time-of-flight spectrometer. This results in interferences over a much
narrower bandwidth, therefore restricting the degrees of freedom of the single
photon mode-function. Notably this method has the advantage of being less
sensitive to higher order phase mismatch between both sources, such as
dispersion. The measured contrast across all spectral bins was found to be
$90\%$.

##  Appendix C Purity of the heralded states

Since the state $\ket{\psi_{12}}$ from the sources is assumed to be a pure
state, the purity of the heralded states $\ket{\Psi_{jk}}$ is ultimately
dependent on the amount of spectral filtering in the heralding BSM. To assess
this purity, we measure HOM interference between the heralded signal photons
when there is no beamsplitter in the idler arms. In this case, upon a
coincidence detection of the idler photons at $(\Omega_{j},\Omega_{k})$, the
reduced state of the signal photons is separable, and given by

$\displaystyle\hat{\rho}_{j}\otimes\hat{\rho}_{k}=\left(\int\textrm{d}^{2}\omega\rho_{j}(\omega,\omega^{\prime})\right)\left(\int\textrm{d}^{2}\tilde{\omega}\rho_{k}(\tilde{\omega},\tilde{\omega}^{\prime})\right)$
$\displaystyle\hat{a}^{\dagger}_{1}(\omega)\hat{a}^{\dagger}_{2}(\tilde{\omega})\ket{\text{vac}}\bra{\text{vac}}\hat{a}_{1}(\omega^{\prime})\hat{a}_{2}(\tilde{\omega}^{\prime}),$
(75)

where

Figure 15: Left: $P(\tau)$ without removing the constant two photon
contribution from source 1 (dot) and source 2 (square). $P(\tau)$ resembles
Eq.(45) and the fit is obtained by summing the individual fits of $P_{jk}$ as
given by (30). Right: distribution of these background terms as a function of
the heralding frequencies $\Omega_{j},\Omega_{k}$

$\displaystyle\rho_{j(k)}(\omega,\omega^{\prime})=\int\textrm{d}\Omega|t_{j(k)}(\Omega)|^{2}f(\omega,\Omega)f^{*}(\omega^{\prime},\Omega).$
(76)

When the signal photons in this state are incident on a 50:50 beamsplitter,
the expected visibility of the HOM interference is given by Mosley  _et al._
(2008)

$\displaystyle V=\mathrm{Tr}(\hat{\rho}_{j}\hat{\rho}_{k}),$ (77)

and when the idlers are detected in identical frequency bins $(j=k)$, this
becomes

$\displaystyle
V=\mathrm{Tr}(\hat{\rho}_{j(k)}^{2})=\mathcal{P}(\hat{\rho}_{j(k)}),$ (78)

where $\mathcal{P}(\cdot)$ denotes the purity of a state. Thus, for $(j=k)$
the visibility of the HOM dip gives a lower bound on the purity of the state
$\hat{\rho}_{j(k)}$, and by extension, the state $\hat{\rho}_{jk}$. Our
measurements, shown in Fig. 14, indicate that purity of the heralded states is
at least 70%, as evidenced by the HOM visibility along the $j=k$ line. By
comparison, a direct calculation of the expected purity using our experimental
parameters gives $\sim 78\%$. The purity of our heralded state seems to be
dominated by the spectral resolution of our spectrometer. Without spectral
resolution, the purity of the heralded state is about 20% as shown in Fig. 14.

##  Appendix D Background signal

As shown by Eq.(3), the full four photon state in the interferometer (see Fig.
1) contains a contribution from photon pairs emitted by individual sources due
to the stochastic nature of parametric down conversion. These terms contribute
to $P(\tau_{S})$ in the form of interferences that get averaged over the
course of a measurement. It is therefore possible to remove that contribution
from the signal subsequently to the measurement by blocking a source and
recording the rate of four-fold coincidences.

We therefore repeated the measurement of $P_{jk}(\tau_{S})$ with either source
blocked to obtain the constant background signal for each $j,k$ frequencies,
as shown in Fig. 15. This shown that the background terms are similar between
both sources, therefore the two sources are similar. Summing over all the
bins, we can plot on the same scale the contribution of all term in Fig. 15.
The peak corresponds to interferences from $\ket{\psi_{12}}$ while the flat
terms represent $\ket{\psi_{11}}$ and $\ket{\psi_{22}}$. As expected from the
theory, both source contribute to $1/4$ of the full signal. Removing those
backgrounds at $\Omega_{j},\Omega_{k}$ from $P_{jk}$, we obtain the fringes
from the main paper with optimal visibility.

As stated in the main text, we have assumed that the three terms in the full
state Eq. (3) are mutually incoherent. This is because our measurements are
taken over a long timescale of a few hours where the optical phase drifts
significantly and any phase-sensitive interference can be neglected. Over a
shorter time-scale, we can measure this optical phase in real time by
measuring fourfold coincidences with the two beamsplitters present, while
scanning the PZT between both sources. A straightforward calculation taking
into account the full state (3) shows that there is a term that oscillates at
the sum frequency $\omega+\Omega\approx\omega_{p}$, where $\omega_{p}$ is the
pump frequency, corresponding to about $415$ nm in wavelength. In Fig. 16, we
plot the measured interference of the two-fold (red) coincidences against the
four-fold fringes (blue), where the latter can be seen to modulate at twice
the frequency of the two-fold modulation.

Figure 16: Interference fringes in the two-fold (red) and fourfold (blue)
coincidences, obtained in ”real-time”, while scanning a PZT as described in
the text. The four-fold fringes can be seen to occur at twice the frequency of
the two-fold fringes.

##  Appendix E Orthogonal modes

From Eq.(12), we see that the heralded state $\ket{\Psi_{jk}}$ is dependent on
the modes $\ket{\phi_{j}}$ and $\ket{\phi_{k}}$, which, in the pure state
case, results in a heralded joint spectrum (20) dependent on the outer
products $\phi_{j}(\omega_{1})\phi_{k}(\omega_{2})$. For each heralding bin
$j$ and $k$, we label the heralded JSI from (12) as
$F_{n}(\omega_{1},\omega_{2})$, where $n$ indexes a pair $(j,k)$. These are
normalized as $\int\textrm{d}^{2}\omega
F_{n}(\omega_{1},\omega_{2})=1~{}\forall~{}n$ but are not orthogonal, even in
the pure state case, i.e $\int\textrm{d}^{2}\omega
F_{n}(\omega_{1},\omega_{2})F_{m}(\omega_{1},\omega_{2})\neq\delta_{nm}$.
Orthogonality is usually a corner stone in any quantum protocol, and it is
therefore necessary to select the heralded states from our measurement that
are orthogonal. To do so, we utilize our measurement of $F_{jk}$, obtained by
measuring the spectral coincidences between the signal’s photon heralded by a
BSM on the idlers. We then obtain a figure similar to Fig. 7 albeit without
perfect spectral resolution, putting us in the mixed state configuration, but
the strategy to select orthogonal modes within this set is similar to the pure
state model.

Figure 17: Sets of orthogonal modes $F_{jk}$ that have less than 15 % of
mutual overlap within each part of the grid. The insets label $j$ and $k$.

First, it is important to notice the symmetry in (12), where $F_{jk}=F_{kj}$
for $j\neq k$. Since our TOFS are well-calibrated, it is reasonable to
symmetrize our measured heralded JSI by averaging the experimentally obtained
$F_{jk}$ and $F_{kj}$ (for $j\neq k$) thus defining the $F_{n}$ functions.
Then we compute the mutual overlaps $\int\textrm{d}^{2}\omega
F_{n}(\omega_{1},\omega_{2})F_{m}(\omega_{1},\omega_{2})$ and use an algorithm
to select a set of modes $\left\\{F_{n}\right\\}$ which all have an overlap
below a certain threshold of $15\%$. We represented a few of these JSI in Fig.
17. Since the spectral range of our high resolution TOFS is limited, so is the
range over which we can compute overlap, as can be seen from the modes that
are labelled with a large $j,k$. Nevertheless, there is a sufficient amount of
spectral coincidence in those cases to infer orthogonality with the other JSI.

Note that while this overlap is computed between the joint spectral
intensities and not between the states, it can be shown that if the overlap in
intensity is zero, then the states are necessarily orthogonal, hence the
strategy is valid to select which $\ket{\Psi_{jk}}$ are mutually orthogonal.
Therefore, it is fair to say that the JSI $F_{jk}$ from Fig. 17 correspond to
the heralded states $\ket{\Psi_{jk}}$ that are all mutually orthogonal.

##  Appendix F Miscellaneous functions and relations

Our theoretical derivation relies on the definition of the
$\phi_{j(k)}(\omega)$ functions which renders computation easier thanks to the
Gaussian approximation. These functions can be also written in a density
matrix formalism.
Using the definitions from Sec.II, the $\phi_{j(k)}$ functions are defined
from the JSA by:

$\displaystyle f(\omega,\Omega_{j(k)}=\sqrt{N_{j}}\phi_{j(k)}(\omega),$ (79)

where $N_{j}$ is a function that depends on the heralding frequency
$\Omega_{j}$, the JSA bandwidth over the idler axis $\sigma_{\text{I}}$ and
the amount of entanglement $\alpha$. The reduced density matrix of the idler,
given by Eq.(34), can then be written in the following manner:

$\displaystyle\rho_{\mathrm{I}}(\Omega_{j},\Omega_{k})$
$\displaystyle=\int\textrm{d}\omega\
f(\omega,\Omega_{j})f^{\ast}(\omega,\Omega_{k})$ (80)
$\displaystyle=\sqrt{N_{j}N_{k}}\int\textrm{d}\omega\
\phi_{j}(\omega)\phi_{k}^{\ast}(\omega)$ (81)
$\displaystyle=\sqrt{N_{j}N_{k}}\braket{\phi_{j}}{\phi_{k}},$ (82)

where we see that the idler density matrix can be linked to the overlap
integral between the heralded signal states. The $N_{j(k)}$ functions are then
found to be equal to the diagonal elements of the idlers density matrix:

$\displaystyle
N_{j(k)}=\rho_{\mathrm{I}}\left(\Omega_{j(k)},\Omega_{j(k)}\right),$ (83)

since the $\phi_{j(k)}$ are $\ell^{2}$ normalized.

The signal density matrices then follow a similar derivation, with

$\displaystyle\rho_{\mathrm{S}}(\omega,\omega^{\prime})$
$\displaystyle=\int\textrm{d}\Omega_{j(k)}\
f(\omega,\Omega_{j(k)})f^{\ast}(\omega^{\prime},\Omega_{j(k)})$ (84)
$\displaystyle=\int\textrm{d}\Omega_{j(k)}\
N_{j(k)}\phi_{j(k)}(\omega)\phi_{j(k)}^{\ast}(\omega),$ (85)

and the diagonal elements are given by

$\displaystyle\rho_{\mathrm{S}}(\omega,\omega)=\int\textrm{d}\Omega_{j(k)}\
N_{j(k)}\left|\phi_{j(k)}(\omega)\right|^{2}.$ (86)

It is easy to show that in the case of similar sources, we also have the
identity

$\displaystyle\rho_{\mathrm{S}}(\omega,\omega^{\prime})=\rho_{\mathrm{S}}^{\ast}(\omega^{\prime},\omega).$
(87)

Finally, the overlap between the modes $\phi_{j}$ and $\phi_{k}$ can also be
written in term of the density matrices, as:

$\displaystyle\braket{\phi_{j}}{\phi_{k}}=\frac{\rho_{\mathrm{I}}(\Omega_{j},\Omega_{k})}{\sqrt{N_{j}N_{k}}}$
(88)

## References

  * Kimble (2008) H. J. Kimble, Nature 453, 1023 (2008).
  * Brecht  _et al._ (2015) B. Brecht, D. V. Reddy, C. Silberhorn, and M. G. Raymer, Physical Review X 5, 041017 (2015).
  * Jozsa and Linden (2003) R. Jozsa and N. Linden, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 459, 2011 (2003).
  * Ekert (1991) A. K. Ekert, Phys. Rev. Lett. 67, 661 (1991).
  * de Riedmatten  _et al._ (2005) H. de Riedmatten, I. Marcikic, J. A. W. van Houwelingen, W. Tittel, H. Zbinden, and N. Gisin, Phys. Rev. A 71, 050302 (2005).
  * Raymer  _et al._ (2013) M. G. Raymer, A. H. Marcus, J. R. Widom, and D. L. Vitullo, Journal of Physical Chemistry B 117, 15559 (2013).
  * Ansari  _et al._ (2018) V. Ansari, J. M. Donohue, B. Brecht, and C. Silberhorn, Optica 5, 534 (2018), 1803.04316 .
  * Vértesi and Navascués (2011) T. Vértesi and M. Navascués, Phys. Rev. A 83, 062112 (2011).
  * Bennett  _et al._ (1993) C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, and W. K. Wootters, Phys. Rev. Lett. 70, 1895 (1993).
  * Żukowski  _et al._ (1993) M. Żukowski, A. Zeilinger, M. A. Horne, and A. K. Ekert, Phys. Rev. Lett. 71, 4287 (1993).
  * Pan  _et al._ (1998) J.-W. Pan, D. Bouwmeester, H. Weinfurter, and A. Zeilinger, Physical Review Letters 80, 3891 (1998).
  * Wagenknecht  _et al._ (2010) C. Wagenknecht, C.-M. Li, A. Reingruber, X.-H. Bao, A. Goebel, Y.-A. Chen, Q. Zhang, K. Chen, and J.-W. Pan, Nature Photonics 4, 549 (2010).
  * Blow  _et al._ (1990) K. J. Blow, R. Loudon, S. J. D. Phoenix, and T. J. Shepherd, Phys. Rev. A 42, 4102 (1990).
  * Ramelow  _et al._ (2009) S. Ramelow, L. Ratschbacher, A. Fedrizzi, N. K. Langford, and A. Zeilinger, Physical Review Letters 103, 253601 (2009).
  * Graffitti  _et al._ (2020) F. Graffitti, P. Barrow, A. Pickston, A. M. Brańczyk, and A. Fedrizzi, Phys. Rev. Lett. 124, 053603 (2020).
  * Fedrizzi  _et al._ (2009) A. Fedrizzi, T. Herbst, M. Aspelmeyer, M. Barbieri, T. Jennewein, and A. Zeilinger, New Journal of Physics 11, 103052 (2009).
  * Treps  _et al._ (2005) N. Treps, V. Delaubert, A. Maître, J. M. Courty, and C. Fabre, Phys. Rev. A 71, 013820 (2005).
  * Roslund  _et al._ (2014) J. Roslund, R. M. De Araujo, S. Jiang, C. Fabre, and N. Treps, Nature Photonics 8, 109 (2014).
  * Thiel  _et al._ (2017) V. Thiel, J. Roslund, P. Jian, C. Fabre, and N. Treps,  2, 034008 (2017).
  * Torres-Company  _et al._ (2011) V. Torres-Company, J. Lancis, and P. Andres, Progress in Optics 56, 1 (2011).
  * Goda and Jalali (2013) K. Goda and B. Jalali, Nature Photonics 7, 102 (2013).
  * Davis  _et al._ (2017) A. O. C. Davis, P. M. Saulnier, M. Karpiński, and B. J. Smith, Opt. Express 25, 12804 (2017).
  * Kim  _et al._ (2015) H. Kim, S. M. Lee, and H. S. Moon, Scientific reports 5, 1 (2015).
  * Merkouche  _et al._ (2021) S. Merkouche, V. Thiel, A. Davis, , and B. Smith, (2021), in preparation.
  * Davis  _et al._ (2020) A. O. C. Davis, V. Thiel, and B. J. Smith, Optica 7, 1317 (2020).
  * Mosley  _et al._ (2008) P. J. Mosley, J. S. Lundeen, B. J. Smith, P. Wasylczyk, A. B. U’Ren, C. Silberhorn, and I. A. Walmsley, Phys. Rev. Lett. 100, 133601 (2008).