Sample-efficient adaptive calibration of quantum networks
using Bayesian optimization

We envision using a GP Bayesian optimization algorithm rasmussen2006 to maximize the indistinguishability of photons generated at a quantum node and, after travelling through deployed fiber optics cable, arriving at a two-photon measurement station, as sketched on the right hand side of Fig. 1. In addition to being a powerful optimization framework, the final GP surrogate model allows for a deep analysis of a high-dimensional parameter space which, as described below, corresponds to degrees of freedom of the photons in this work. The calibration problem consists of optimizing an expensive objective function, $f(\mathbf{x})$, corresponding to a measure of the change in coincidence detections due to two-photon interference, as described in the next section. GP Bayesian optimization is based on the assumption that the likelihood and prior distributions correspond to multivariate Gaussian distributions, where the objective function is assumed to be a random variable sampled from a Gaussian Process,

$$f(\mathbf{x})\sim\mathcal{GP}[\mu(\mathbf{x}),K(\mathbf{x},\mathbf{x}^{\prime})],$$

(1)

which is defined by the mean function, $\mu(\mathbf{x})=\mathbb{E}[f(\mathbf{x})]$, and covariance (or kernel) function, $K(\mathbf{x},\mathbf{x}^{\prime})=\mathbb{E}[(f(\mathbf{x})-\mu(\mathbf{x})(f(\mathbf{x^{\prime}})-\mu(\mathbf{x^{\prime}}))]$ respectively. The Bayesian optimization algorithm proceeds by feeding the optimizer, a single function call which is equivalent to a single measurement result, $f(\mathbf{x}^{(i)})$, during the $i$th iteration of the algorithm,. The variable $\mathbf{x}^{(i)}=(x_{1}^{(i)},x_{2}^{(i)},\cdots,x_{D}^{(i)})$ is a $D$-dimensional vector corresponding to $D$ experimental parameters. For two-photon interference, these parameters control the photonic degrees of freedom that are adjusted to compensate for channel transformations. For our proof-of-principle laboratory demonstration described in the next section, these correspond to a time delay controlling the mutual arrival time of the photons, half-waveplate angle controlling the polarization of the photons, and spectral filter pass-band determining the frequency detuning (offset) of the photons. Once the $i$th measurement result has been fed into the Bayesian optimizer, it builds a Gaussian Process surrogate model by performing an optimization with respect to kernel hyperparameters (see Appendix .1). An acquisition function, such as the lower confidence bound, which we found worked best in this work (see Appendix .5), is then used to suggest the next measurement, $\mathbf{x}^{(i+1)}$. Various initial sampling strategies, kernel functions, and hyper-parameter optimization strategies are also tested with details of the algorithmic fine-tuning provided in Appendices .4-.6.

To test the algorithm, we performed a HOM interference measurement of two photons. When two indistinguishable independent photons are sent to different input ports of a beam splitter, they will bunch at one of the output ports Hong et al. (1987). This results in a minimum (zero without imperfections) of the zero-delay-time normalized second-order correlation function Beck (2007),

$\displaystyle g^{(2)}(0)=\frac{C_{12}T}{S_{1}S_{2}\Delta\tau},$

(2)

where $S_{1}$ and $S_{2}$ ($C_{12}$) denote the total number of photons detected at each (in coincidence at both) of the outputs of the beamsplitter within time interval $\Delta\tau$ over a total amount of time $T$. This quantity is related to the commonly used HOM interference visibility $V=(C_{12}^{\mathrm{max}}-C_{12}^{\mathrm{min}})/C_{12}^{\mathrm{max}}$, where $C_{12}^{\mathrm{max}}$ ($C_{12}^{\mathrm{min}}$) denotes the maximum (minimum) number of photons detected in coincidence as a degree of freedom of a photon is varied (e.g. its time of arrival) Hong et al. (1987); Valivarthi et al. (2020). Either $V$ or $g^{(2)}(0)$ are used to quantify the impact of imperfections in HOM interference, which ultimately impacts Bell-state measurement fidelities, and hence the fidelity of qubit distribution in quantum networks Nielsen and Chuang (2002); Metcalf et al. (2014); Valivarthi et al. (2020).

In our experiment, which is depicted in Fig. 2,a continuous-wave laser emits near-visible-wavelength light that pumps a periodically-poled lithium niobate crystal waveguide to create a photon pair at telecommunication wavelength through spontaneous parametric down conversion (SPDC) Burnham and Weinberg (1970). The leftover pump light is removed with a 50 GHz bandpass filter, and the photon pair is sent to an initial beam-splitter to probabilistically separate the photons. Here we employ photons originating from the same source to demonstrate our principle even though photons would originate from independent sources in a quantum network. To avoid first-order interference between correlated photons generated by the SPDC process, one of them passes through a fiber-based phase-modulator used to randomize the phase difference between the photon pair Jin et al. (2013). The second photon goes through a free-space optics setup in which a set of half- and quarter-waveplates adjust the photon polarization while mirrors on a translation stage vary the path length difference (corresponding to a mutual time delay) between each of the photons. Each photon then passes an optional, independent, and tunable 12 GHz-bandwidth bandpass filter. The waveplates, translation stage, and tunable filters, surrounded by a dashed outline in Fig. 2, are adjusted according to the algorithm. Next, each photon passes through fiber-based polarizing beam-splitters (PBSs), which ensure that the polarizations of both photons are identical and that any polarization rotations are converted into intensity variations. Then, each photon is directed into a different input port of a second fiber-based beamsplitter at which two-photon interference occurs. Finally, the output photons are guided to two cryogenically cooled superconducting nanowire threshold single-photon detectors. A time-to-digital converter records the time-of-arrival of the photons at the detectors.

In our first measurement, we remove the tunable band-pass filters such that the path-length difference and the half-waveplate angle correspond to the experimental parameters $\mathbf{x}$ that are fed into the Bayesian optimization algorithm. It is important to point out that there are several choices for possible objective functions $f(\mathbf{x})$ one may use for the optimization algorithm. Perhaps unsurprisingly, we find that

$\displaystyle f(\mathbf{x})\equiv g^{(2)}(0)[\mathbf{x}]$

(3)

provides the best performance when compared to others and is what we employ for our scheme. For example, $C_{12}$ can vary for experimental settings that do not yield quantum interference, or $V$ which requires adjustment of experimental parameters to assess both $C_{12}^{\mathrm{min}}$ and $C_{12}^{\mathrm{max}}$. Note here that the zero in $g^{(2)}(0)$ indicates that we always measure about the same relative time-of-arrival of the photons (about $\Delta\tau=2$ ns) when we vary the degrees of freedom, specifically a time that corresponds to a minimum $g^{(2)}(0)$ when photons are rendered indistinguishable.

Before discussing the results of our optimization, note that we develop a theoretical model for $g^{(2)}(0)$ to validate our experimental measurements and the predictions of the Bayesian optimization algorithm (see Appendix Fig. 4). The degenerate output of our Type-0 SPDC crystal is described by a squeezed vacuum state, which is a Gaussian state, i.e. it is completely characterized by a displacement vector and a covariance matrix. Since all optical operations in the experiment, including the photon detection, can be described as Gaussian operations, that is, they map Gaussian states to other Gaussian states, we can apply a characteristic function formalism LaukInPrep ; Takeoka et al. (2015) to determine the final displacement vector and the final covariance matrix, which allow us to predict $S_{1}$, $S_{2}$ and $C_{12}$. In this model, the experimental degrees of freedom are modeled as virtual beam splitters with variable transmittances $\eta_{u},\eta_{d}$, and $\zeta$, where $\eta_{u/d}$ corresponds to overall photon coupling efficiency in the upper/lower path (see Fig. 2) and $\zeta=\exp(-x^{2})$ corresponds to the mode overlap parameterization between the two photons. Details of the model are provided in the Appendix .9 along with the theoretical plot of $g^{(2)}(0)$ in Appendix Fig. 5a.

Our optimization using two degrees of freedom is depicted in Fig. 3a, presenting a two-dimensional map of the measured $g^{(2)}(0)[\mathbf{x}]$ as a function of the full range of settings $\mathbf{x}$ of the translation stage and half-waveplate angle - using 60 s for changing the parameter settings and $T=60$ s of data-collection per setting. Our result shows good agreement with that predicted by the theoretical model as plotted in Appendix Fig. 5a. In particular, $g^{(2)}(0)[\mathbf{x}]$ displays $V\leq 0.5$, as expected for a squeezed vacuum input state as the input to the initial beam splitter, which splits the pair probabilistically. Unlike the $g^{(2)}(0)[\mathbf{x}]$ predicted by the theoretical model, the experimental data exhibits random variations due to sampling noise.

To benchmark the Bayesian optimization algorithm, we used the baseline data in Fig. 3a as the source of values for our objective function $g^{(2)}(0)[\mathbf{x}]$ which are fed to the algorithm. Starting with a set of twelve measurement settings $[\mathbf{x}^{(1)},\mathbf{x}^{(2)},\cdots,\mathbf{x}^{(12)}]$ selected based on Latin hypercube sampling (see Appendix .6 for details), the algorithm proceeds to select the next setting $\mathbf{x}^{(13)}$, and then the next and so on, for a total number of measurements, $n$, that it deems optimal for establishing the minimum of $g^{(2)}(0)[\mathbf{x}]$. As shown in Fig. 3b, we find that the Bayesian optimization algorithm reliably converges to the minimum $g^{(2)}[\mathbf{x}]=20.2$ in less than $n=30$ measurements, on average, (dark blue line). The light blue lines correspond to single instances of the Bayesian optimization algorithm, which converge more slowly or more quickly depending on the random instance of initially sampled points. By benchmarking the algorithm with over one-hundred independent trials, we are able to remove the effect of randomness, resulting in the average (expected) convergence shown by the dark blue line.

While this benchmark confirms the validity of the approach, it is important to test the efficacy of our algorithm in a live setting where every measurement is performed in real-time, as opposed to using the tabulated baseline data from Fig. 3a. The results for four different real-time trials are shown in Fig. 3c. The plots for trials 1-3, for which data was acquired for $T=60$ s (same as for the baseline data), show that the minimum value of $g^{(2)}(0)[\mathbf{x}]$ is found within $n=15$ in all three cases. The maps on the right of Fig. 3c show the final Gaussian Process surrogate model prediction made after $n=30$ measurements. It is worth noting that the surrogate models do not, and indeed are not intended to, accurately reproduce $g^{(2)}(0)[\mathbf{x}]$ over the entire parameter range of the input variables. The surrogate models need only accurately predict the location of the minimum of $g^{(2)}(0)[\mathbf{x}]$. In fact, this pursuit of incomplete but sufficient information is the underlying reason why GP Bayesian optimization is more efficient than brute-force mapping of $g^{(2)}(0)[\mathbf{x}]$. In the fourth trial we reduce $T$ to just 1 s, which significantly increases the level of shot noise. Yet, the GP optimization still finds the minimum value after just $n=23$, albeit the polarization setting is slightly off the optimum value. As the PBS transmission follows a cosine of the half-waveplate angle $g^{(2)}(0)[\mathbf{x}]$ becomes less sensitive to small variations and, thus, renders the optimization more difficult. However, the consequence of a slight offset of the polarization angle is also less significant due to the cosine dependence of the loss. On the other hand, without the PBS the distinguishability would be directly proportional to the offset of the half-waveplate angle.

Our results using two degrees of freedom validate the effectiveness of the algorithm for a realistic signal that is subject to various sources of noise including shot noise, while demonstrating fast convergence. In comparison, the full experimental map (shown in Fig. 3a) consists of 350 data points acquired over 13 h.

Our next measurements test the Bayesian optimization algorithm in a higher-dimensional setting. Specifically, we add the 12 GHz-bandwidth tunable filters to the experimental setup so as to include the frequency offset (between 0 and 6 GHz) of the photons as a third degree of freedom. At 0 GHz offset, both filters are resonant with the degenerate spectral mode of the photon pair, while at the maximum detuning of 6 GHz between the two photons, the spectral filters transmit correlated but partially non-degenerate photons. Compared to the previous setup, a broader region of interference with respect to varied path-length difference is expected due to the spectral filtering. Since the parameter space of settings in three degrees of freedom is too large to acquire a baseline data-set as that in Fig. 3a, we proceed directly to a live demonstration, with results shown in Fig. 4a. We find excellent convergence towards the minimum value within $n=25$ measurements. Two-dimensional partial-dependence plots of the thirty-point surrogate model are shown in Fig. 4b. The partial dependence plots illustrate the surrogate model prediction with respect to two degrees of freedom, with the remaining degree of freedom averaged out. We note that the first two degrees of freedom, corresponding to the half-waveplate angle and path length difference, display a near minimum at around 85 degrees and 7.5 mm respectively, as expected. The final degree of freedom, corresponding to the frequency offset, displays a near-constant dependence with a small slope predicting a minimum at $X_{2}=0$ (zero-frequency offset).

Using our theoretical model, we generate a map of the predicted $g^{(2)}(0)$ for all possible parameter settings. Note that our model captures the spectral offset of the two photons as a transition between degenerate squeezed vacuum to non-degenerate two-mode squeezed vacuum states as inputs to the final beamsplitter, as shown in Appendix Fig. 5a,b. Thus, our theoretical model, instead of an experimentally acquired baseline map, is used as the source of the $g^{(2)}(0)[\mathbf{x}]$ values employed for systematically benchmarking the Bayesian optimization algorithm for the three-dimensional parameter set. In Fig. 4c we plot the value of the predicted $g^{(2)}(0)[\mathbf{x}]$ as a function of the total number of measurements $n$ averaged over a hundred trials. A set of fifteen measurement settings based on Latin hypercube sampling are used as the initial sampling points. The results show excellent convergence towards the predicted minimum within $n=40$ measurement calls, showing the efficacy of the result in the higher dimensional setting while also demonstrating the utility of the theoretical model.

For deployed quantum networks, where photon sources are separated, one will not interfere correlated photons as in our demonstration. Instead, two-photon measurements are typically performed on two uncorrelated coherent photons, as in the case of MDI-QKD Lo et al. (2012), or two photons, each entangled with another photon not partaking in the measurement, as in the case of quantum repeaters Briegel et al. (1998). For the latter case, photons from separate SPDC sources will be independently calibrated in order to maximize the degree of photon indistinguishability. When one photon from an SPDC source is observed without its correlated partner, it will follow a thermal photon number distribution Burnham and Weinberg (1970). Hence, we provide additional benchmarks with thermal states at the input of the beamsplitter – corresponding to two independent members of separate photon pairs – using our theoretical model for the predicted $g^{(2)}(0)$, see Appendix Fig. 5c. We simulate the effect of finite sampling by performing Poisson sampling of the theoretical model, which we define as the baseline data. The independent degrees of freedom correspond to $\zeta=\exp(-x^{2})$ (due to path length difference) and $\eta_{u}$ (loss due to polarization misalignment), which in an experimental setting would be proportional to the cosine of the polarization, see the Appendix .9 for more details. The result is shown in Fig. 5a. The effect of finite sampling is evident by the random noise, which reflects what a real acquired signal would look like, e.g. as in Fig. 3a. The baseline objective function is, thus, used to test the Bayesian optimizer under experimental conditions for which the photon detection rate can be very limited. The results shown in Fig. 5 provide the final benchmark results for Bayesian optimization using thermal input states, averaged with 150 different instances. As before, a set of twelve measurement settings based on Latin hypercube sampling are used as the initial sampling points. The results shown in Fig. 5c show that the Bayesian optimizer finds the expected minimum within $n=30$ measurements (located near $\eta_{u}=0.2$ and $x=0$). These results are promising for in situ quantum network calibration.

Bayesian optimization is a sample-efficient technique that performs sequential optimization of time-consuming black-box functions. There are two main steps to any Bayesian optimization algorithm: (1) the construction of a conditional probabilistic model for the objective function based on a set of observations, and (2) the construction of an acquisition function that uses this model to predict future observations that optimize the objective function of interest. The first step requires an understanding of Bayesian inference, which quantifies how we can update our belief of a particular hypothesis of the objective function based on the current set of observations. The second step requires defining an effective criterion that may be used to predict new observation points that will (most likely) be close to the optimum. In the following section, we will introduce basic concepts related to Bayesian inference. The second section will define Gaussian Processes (GP) which are often used as computationally tractable surrogate models for Bayesian optimization. The third section will present various kernel functions relevant to Bayesian optimization. The fourth section will define and discuss various acquisition functions which are often used in practice. The fifth section will define the quantum network black-box objective used for this study.

Bayesian inference aims to construct a conditional probability distribution $p(\mathcal{H}|\mathcal{D})$ for a hypothesis $\mathcal{H}$ based on observed data $\mathcal{D}$. For our purposes, the hypothesis represents the real values of an objective function $f(\mathbf{x})$ where $\mathbf{x}$ is a set of experimental knobs/parameters within the quantum optical experiment. In the context of quantum networks, we are interested in maximizing photon indistinguishability in order to maximize the performance of a wide variety quantum network operations and protocols. A single-measurement objective function which is able to quantify the degree of photon indistinguishability is the normalized second order photon correlation function, $f(\mathbf{x})\equiv g^{(2)}(0)[\mathbf{x}]$, measured with the Hong-Ou-Mandel experimental set-up shown in the main text. Here, the goal is to minimize the second order correlation function in order to maximize indistinguishability. For general Bayesian optimization, the experimental knobs can be categorical, integer, or real-valued quantities but for our purposes, we will only consider real-valued quantities such as the polarization rotation angle, position delay stage, and frequency filter bandwidth.

In the following, we define the basic concepts of Bayesian inference without reference to the photon correlation function, however, we will later refine the likelihood functions and kernel functions for the Gaussian Process in order to refine the Bayesian optimization algorithm to the quantum network problem.

For arbitrary likelihood and prior distributions, the exact calculation of the posterior distribution is not tractable. In particular, the denominator, which acts as a normalization constant for the posterior to remain a valid probability distribution, requires the calculation of a high-dimensional integral that is generally intractable. It is possible to use Monte-Carlo-based methods, however, these approaches are generally computationally expensive and ultimately become intractable for high dimensions. By assuming that the likelihood and prior correspond to multivariate Gaussian distribution functions, it is possible to calculate the posterior and predictive posterior exactly. In the Gaussian Process framework, we assume that the training and test data have Gaussian noise and may be written in the form,

$$y^{(i)}=f(\mathbf{x}^{(i)})+\epsilon^{(i)},$$

(A.3)

where $\epsilon_{i}$ corresponds to Gaussian-distributed noise with zero mean, $p(\epsilon)=\mathcal{N}(0,\sigma^{2})$. The likelihood of the observed data $\mathcal{D}$ is then given by:

$$p(\mathbf{y}|X,\mathbf{f})=\mathcal{N}(\mathbf{f},\sigma^{2}I).$$

(A.4)

The “noiseless” signal $\mathbf{f}=[f(\mathbf{x}^{(1)},f(\mathbf{x}^{(2)}),\cdots,f(\mathbf{x}^{(n)})]^{T}$ plays the role of the hypothesis $\mathcal{H}$ which corresponds to the parameters we wish to estimate. In this regard, the Gaussian Process approach is model-agnostic, which may be advantageous in scenarios where parametric models are not available, or perhaps undesirable to avoid the introduction of model bias. Since the noiseless vector $\mathbf{f}$ plays the role of the effective parameters, we require defining a prior distribution on these parameters. The Gaussian Process framework assumes the prior distribution is a Gaussian Process written as, $f(\mathbf{x})\sim\mathcal{GP}(\mu(\mathbf{x}),k(\mathbf{x},\mathbf{x}^{\prime}))$, where the symbol “$\sim$” is read as: “is sampled from” or “is distributed as.” This implies that we treat the measured outputs $f_{i}$ as random variables that are sampled from the Gaussian process distribution function. Note that a Gaussian process is a distribution over functions completely defined by the mean function $\mu(\mathbf{x})=\mathbb{E}[f(\mathbf{x})]$ and covariance function $k(\mathbf{x},\mathbf{x}^{\prime})=\mathbb{E}[(f(\mathbf{x})-\mu(\mathbf{x}))(f(\mathbf{x}^{\prime})-\mu(\mathbf{x}^{\prime}))]$. The prior distribution for $\mathbf{f}$ can therefore be written as a multivariate Gaussian,

$$p(\mathbf{f})=\mathcal{N}(\mu(\mathbf{x}),K)$$

(A.5)

where $K$ is the kernel covariance matrix,

	$\displaystyle K$	$\displaystyle\equiv K(X,X)$		(A.6)
		$\displaystyle=\begin{bmatrix}k(\mathbf{x}^{(1)},\mathbf{x}^{(1)})&k(\mathbf{x}^{(1)},\mathbf{x}^{(2)})&\cdots&k(\mathbf{x}^{(1)},\mathbf{x}^{(n)})\\ k(\mathbf{x}^{(2)},\mathbf{x}^{(1)})&k(\mathbf{x}^{(2)},\mathbf{x}^{(2)})&\cdots&k(\mathbf{x}^{(2)},\mathbf{x}^{(n)})\\ \cdots&\cdots&\cdots&\cdots\\ k(\mathbf{x}^{(n)},\mathbf{x}^{(1)})&k(\mathbf{x}^{(n)},\mathbf{x}^{(2)})&\cdots&k(\mathbf{x}^{(n)},\mathbf{x}^{(n)})\end{bmatrix},$

and must be a positive semi-definite matrix. Given the likelihood (A.4) and prior distribution (A.5), it is possible to calculate the posterior as:

	$\displaystyle p(\mathbf{f}|X,\mathbf{y})$	$\displaystyle=\frac{p(\mathbf{y}|X,\mathbf{f})p(\mathbf{f})}{p(\mathbf{y}|X)}$		(A.7)
		$\displaystyle=\mathcal{N}(K(K+\sigma^{2}I)^{-1}\mathbf{y},\sigma^{2}(K+\sigma^{2}I)^{-1}K)$

where, as an aside, it is worth noting that $\sigma^{2}(K+\sigma^{2}I)^{-1}K=K-K(K+\sigma^{2}I)^{-1}K$, and the latter quantity is often found in textbooks. Let us now calculate the predictive distribution for unobserved data, $\mathcal{D}^{*}=\{X^{*},\mathbf{y}^{*}\}$. The joint Gaussian process prior for the observed and unobserved data is written as:

$$p\begin{pmatrix}\mathbf{f}\\ \mathbf{f}^{*}\end{pmatrix}=\mathcal{N}\left[\begin{pmatrix}\mu(\mathbf{x})\\ \mu(\mathbf{x}^{*})\end{pmatrix},\begin{pmatrix}K_{x,x}+\sigma^{2}I&K_{x,x^{*}}\\ K_{x^{*},x}&K_{x^{*},x^{*}}\end{pmatrix}\right].$$

(A.8)

The predictive distribution, $p(\mathbf{f}^{*}|\mathcal{D})=\int p(\mathbf{f}^{*}|\mathbf{f})p(\mathbf{f}|\mathcal{D})d\mathbf{f}$, is then found to be:

$$p(\mathbf{f}^{*}|\mathcal{D})=\mathcal{N}(\mathbf{f}^{*}|\mu_{p}(\mathbf{x}^{*}),\Sigma_{p}(\mathbf{x}^{*}))$$

(A.9)

with $\mu_{p}(\mathbf{x}^{*})$ and $\Sigma_{p}(\mathbf{x}^{*})$ corresponding to the predicted mean and variance of the model at point $\mathbf{x}^{*}$,

	$\displaystyle\mu_{p}(\mathbf{x}^{*})$	$\displaystyle=\mu(\mathbf{x}^{*})+k_{*}^{T}(K+\sigma^{2}I)^{-1}(\mathbf{f}-\mu(\mathbf{x}))$		(A.10)
	$\displaystyle\Sigma_{p}(\mathbf{x}^{*})$	$\displaystyle=k(\mathbf{x}^{*},\mathbf{x}^{*})-k_{*}^{T}(K+\sigma^{2}I)^{-1}k_{*}$		(A.11)

where $k_{*}=k(\mathbf{x},\mathbf{x}^{*})$ is an $(n\times m)$ matrix corresponding to the $n$ training points $\mathbf{x}$ and $m$ test points $\mathbf{x}^{*}$. Equations (A.10) and (A.11) are the key equations which describe Gaussian-Process-based Bayesian inference.

The covariance function encodes information about the shape and structure the objective function. The kernel ultimately affects the convergence rate of the Bayesian optimization algorithm. Below, we write several well-known kernel functions from the literature. For example, the squared exponential kernel is written as:

$$k_{(}\mathbf{x},\mathbf{x}^{\prime})=\exp\left(-\frac{||\mathbf{x}-\mathbf{x}^{\prime}||^{2}}{2\ell^{2}}\right),$$

(A.12)

which ensures that nearby points have similar function values within the length scale given by $\ell$. The periodic kernel is given by:

$$k(\mathbf{x},\mathbf{x}^{\prime})=\exp\left(-\frac{2\sin^{2}(\pi||\mathbf{x}-\mathbf{x}^{\prime}||/p)}{\ell^{2}}\right)$$

(A.13)

where $\ell$ is a length scale parameter and $p$ is the periodicity of the kernel. Finally, we also consider the Matern kernel,

$\displaystyle k(\mathbf{x},\mathbf{x}^{\prime})=\frac{1}{\Gamma(\nu)2^{\nu-1}}\left(\tfrac{\sqrt{2\nu}}{l}||\mathbf{x}-\mathbf{x}^{\prime}||^{2}\right)^{\nu}\!\!K_{\nu}(\tfrac{\sqrt{2\nu}}{l}||\mathbf{x}-\mathbf{x}^{\prime}||^{2}),$

(A.14)

where $K_{\nu}$ is a modified Bessel function and $\Gamma$ is the gamma function. The parameter $\nu$ controls the smoothness of the function. The smaller the $\nu$, the less smooth the function will be. A comparison of various kernels is shown in Appendix Fig. 1, demonstrating that the Matern kernel provides the best performance compared to all other kernels.

An acquisition function $a(x)$ aims to evaluate the expected loss associated with evaluating $f(\mathbf{x})$ at point $\mathbf{x}$, and selects the point with the lowest expected loss. In the following, we compare three different acquisition functions which are often used in the literature.

While Bayesian optimization can start with zero knowledge of the objective function (i.e. the size of training data is null), it is possible to accelerate the convergence of the Bayesian optimization algorithm by using a set of judiciously chosen initial sampling points. Conventionally, there are many techniques available for sampling such as: uniform or random sampling, Sobol sampling, Halton sampling, and Latin hypercube sampling (lhs). By using the skopt initial sampling package, we performed additional benchmarks with respect to the number of initial sampling points as well as the type of sampling used. In Appendix Fig. 3, we provide a comparison between Sobol, Halton, random, conventional Latin hypercube sampling and Maximin Latin hypercube sampling showing that the maximin and conventional Lain hypercube sampling provided the best performance. For live testing purposes, we opted to use the conventional Latin Hypercube Sampling scheme, however, further work will be required to determine optimal initial sampling methods for other response functions of interest for quantum network calibration experiments.

The proposed Bayesian optimization algorithm can be summarized by the following steps:

1. 1.

   Perform $n$ initial measurements $\mathbf{y}=(y_{1},\cdots,y_{n})$ with respect to the experimental degrees of freedom, $X=[\mathbf{x}_{1},\cdots,\mathbf{x}_{n}]$, using Maximin latin hypercube sampling.

2. 2.

   Build Gaussian Process surrogate model based on existing measurements.

3. 3.

   Suggest new measurement location $\mathbf{x}_{i+1}$ based on the acquisition function prediction.

4. 4.

   Repeat 2-3 until stopping criterion is met.

The Gaussian Process model built in step 2 of the Bayesian optimization algorithm above can be further decomposed into the following steps:

1. 1.

   Build gram matrix, equation (6), for Matern Kernel ($\nu=5/2$) based on equation (14).

2. 2.

   Optimize the Kernel hyperparameters based on the marginal likelihood, equation (17).

3. 3.

   Provide predictive mean and variances based on equations (A.10) and (A.11).

The numerical implementation and benchmarking of the Bayesian optimization algorithm was done with in-house code written in Python with standard numpy and scipy linear algebra and optimization packages. However, there are a wide variety of excellent open-source software implementations that are available through github. We recommend skopt (Bayesian optimization) and GPytorch (Gaussian Process regression) Gpytorch which provide an excellent starting point for initial testing. It is worth noting that we found the customization of the open-source software difficult (for example, adding custom acquisition functions, kernel functions, as well as other functionalities) which is why we opted to perform the benchmarking (shown in this Appendix) with our in-house code. The live testing was performed with a modified version of the skopt Bayesian optimization algorithm.

To model the experiment we use methods of phase space quantum optics, in particular the characteristic function formalism LaukInPrep ; Takeoka2015 . This method allows us to describe many experimental imperfections, such as coupling losses, non-perfect detector efficiencies, high photon number contributions, photon number statistics, etc., and is best suited to deal with the Gaussian quantum optical states. Gaussian states are the states whose characteristic function, or in fact any other phase space representation, is given by a multidimensional Gaussian function

$\displaystyle\chi(\xi)=\exp(-\frac{1}{4}\xi\gamma\xi-id\xi),$

(A.23)

with $d$ and $\gamma$ being displacement vector and covariance matrix of the system respectively, $\xi\in\mathbb{R}^{2n}$ with $n$ being the number of independent bosonic mode of the system. The examples of the Gaussian states include the vacuum state, coherent states, as well as squeezed vacuum and two-mode squeezed vacuum states. The later two are directly relevant to our experimental studies since they describe the output states of degenerate and non-degenerate spontaneous parametric down conversion process respectively. It is known that linear optic operations, such as beam splitters, phase shifters etc., preserve Gaussian states Weedbrook2012 ; Braunstein2005 , i.e. they map Gaussian states onto Gaussian states. For the characteristic function it means that for any linear optic operation there exists a symplectic matrix $S$ that transforms the initial displacement vector and the covariance matrix to the output form, $\gamma^{\prime}=S^{T}\gamma S$ and $d^{\prime}=Sd$.

The model representation of the experimental setup is shown in Appendix Fig. 4. The output of the SPDC crystal is mixed with a vacuum input on the first beam splitter BS. The overall coupling efficiencies in the upper and the lower arm can be modelled by the virtue of a virtual beam splitter with the transmittance $\eta_{u}$ and $\eta_{d}$ respectively. The polarization degree of freedom can be modeled by a variable transmittance, e.g. $\eta_{u}$, that can take values from zero, corresponding to crossed polarizations, to some finite value $\eta_{u}^{f}$, corresponding to parallel polarizations. Another virtual beam splitter with a transmittance $\zeta$ is used to model the path length or the time of arrival degree of freedom. $\zeta$ corresponds to the mode overlap of the arriving photons and can take values from 0, corresponding to the case of completely distinguishable photon wave packets with zero overlap, to 1, corresponding to the case of completely indistinguishable photons. To resemble the experimental data more closely we parametrize the $\zeta$ parameter as a Gaussian function $\zeta=\exp(-x^{2})$ reflecting the fact that the overlap of two Gaussian pulses is again a Gaussian function. As we can see from the figure only transmitted parts interfere at the following beam splitter whereas reflected parts do not interfere but are mixed with a vacuum input instead. To model the phase averaging we introduce a phase shifter in the lower arm that shifts the relative phase between the two arms by an amount $\phi$. To obtain the experimental values of the (phase dependent) quantities we average over $\phi$ from $0$ to $2\pi$. In the experimental setup this averaging is achieved by the phase-modulator in one of the arms continuously sweeping over a $2\pi$ phase-shift.

After determining individual symplectic transformations and combining them together we obtain the covariance matrix that completely describes the final state. Using the characteristic function description we can now calculate the probabilities for a single and coincidence photon detection as follows:

	$\displaystyle P_{D1/D2}$	$\displaystyle=\mathrm{Tr}\{\hat{\rho}(1-\ket{0}\bra{0}_{1/2})\}=1-\mathrm{Tr}\{\hat{\rho}(\ket{0}\bra{0})_{1/2}\}$		(A.24)
	$\displaystyle P_{D1D2}$	$\displaystyle=\mathrm{Tr}\{\hat{\rho}(1-\ket{0}\bra{0}_{1})(1-\ket{0}\bra{0}_{2})\}$
		$\displaystyle=1-\mathrm{Tr}\{\hat{\rho}(\ket{0}\bra{0})_{1}\}-\mathrm{Tr}\{\hat{\rho}(\ket{0}\bra{0})_{2}\}$
		$\displaystyle+\mathrm{Tr}\{\hat{\rho}(\ket{0}\bra{0})_{1}(\ket{0}\bra{0})_{2}\},$		(A.25)

where $\ket{0}_{1/2}$ corresponds to a vacuum state of all modes impinging on the detector $D1/D2$. From these we can now easily determine the normalized $g^{(2)}(0)$ correlation function as:

$\displaystyle g^{(2)}(0)=\frac{P_{D1D2}}{P_{D1}P_{D2}}.$

(A.26)

We studied the $g^{(2)}(0)$ correlation function for different input states corresponding to different physical situations. If the underlying nonlinear process is a degenerate down conversion, i.e. two identical photons are created, the input state is described by a squeezed vacuum state with the corresponding covariance matrix

$\displaystyle\gamma_{SV}=\begin{pmatrix}\mathrm{e}^{-2r}&0\\ 0&\mathrm{e}^{2r}\end{pmatrix},$

(A.27)

where $r$ is the squeezing parameter. This input corresponds to all of the experimental two dimensional parameter space scans, i.e., the data in Fig. 3 of the main text. In Appendix Fig. 5a the $g^{(2)}$ is plotted as a function of $\zeta=\exp(-x^{2})$ (path length difference) and $\eta_{u}$ (loss due to polarization misalignment of upper arm).

If the down conversion process is non-degenerate, i.e. two different photons are created for example at different frequencies, the input state is a two-mode squeezed vacuum state with the corresponding covariance matrix

$\displaystyle\gamma_{TMSV}=\begin{pmatrix}\cosh^{2}(r)+\sinh^{2}(r)&0&2\cosh(r)\sinh(r)&0\\ 0&\cosh^{2}(r)+\sinh^{2}(r)&0&-2\cosh(r)\sinh(r)\\ 2\cosh(r)\sinh(r)&0&\cosh^{2}(r)+\sinh^{2}(r)&0\\ 0&-2\cosh(r)\sinh(r)&0&\cosh^{2}(r)+\sinh^{2}(r)\end{pmatrix}$

(A.28)

This input corresponds to the case in the three-dimensional parameter space for which the spectral filters in the two arms have been shifted with respect to each other by an amount much greater than the bandwidth of the filters. However, in this case the photons in each arm are clearly distinguishable, so the final HOM interference will always be between distinguishable modes. We can tweak the model to accommodate this feature by introducing another virtual beamsplitter with reflectivity $\zeta^{\prime}$. Essentially, the same result can be achieved by fixing $\zeta=0$. The predicted $g^{(2)}(0)$ is plotted as a function of $\zeta=\exp(-x^{2})$ (path length difference and $\eta_{u}$ (loss due to polarization misalignment of upper arm in Appendix Fig. 5b. The cross-correlation function does not depend on $\zeta$, but only on $\eta_{u}$, since the latter induces actual loss in the PBS, i.e., changes the number of photons reaching the detector. In reality we do not fully achieve the scenario in which the filters are detuned by more than bandwidth as the filters are maximally shifted by 6 GHz with respect to each other while their bandwidth is about 12 GHz. Hence, the actual expected $g^{(2)}$ will be a combination of the squeezed vacuum input case (Appendix Fig. 5a) and the non-degenerate two-mode squeezed state (Appendix Fig. 5b). The relative weight of the two contributions can be found by calculating the spectral overlap of the photons after the filter, which, for Gaussian filter pass-bands, will again yeild a Gaussian shaped weight parameter as a function of the spectral separation of the two filters.

Finally, in a real world implementation of a quantum network the Bell state measurements will occur between two different entangled photon sources. This means that the two partaking photons will originate from different SPDC sources and, thus, not be correlated. The input state in this case will correspond to two thermal states with the corresponding covariance matrix

$\displaystyle\gamma_{TMSV}=\begin{pmatrix}1+2\mu&0&0&0\\ 0&1+2\mu&0&0\\ 0&0&1+2\mu^{\prime}&0\\ 0&0&0&1+2\mu^{\prime}\end{pmatrix},$

(A.29)

where $\mu$ and $\mu^{\prime}$ are the mean photon numbers at the two inputs. Appendix Fig. 5c shows the $g^{(2)}(0)$ function for the thermal input state scenario.