* (85) S. Chakdar and P. Q. Hung, “Prospect of the Electroweak Scale Right-handed neutrino model in the Lifetime Frontier,” in International Conference on Neutrinos and Dark Matter.  5, 2020.  arXiv:2006.00381 [hep-ph].
  * (86) K. Hartling, K. Kumar, and H. E. Logan, “The decoupling limit in the Georgi-Machacek model,” Phys. Rev. D 90 no. 1, (2014) 015007, arXiv:1404.2640 [hep-ph].


$\texttt{r}_{k}$ falls below $1/M$. We do observe that the average fidelity is
slightly greater in the case $2$, especially for $M=10^{3}$. This could be due
to the large variance in the estimator for such a small number of samples. We
also expect to observe the converse for large $M$, and indeed the average
fidelity increases for $M=10^{6}$ and $M=10^{7}$ as $\mathcal{L}$ is lowered
to zero. The last column shows the observed change in the log-likelihood in
the last iteration step before termination for case $2$. It can be seen that
for $M=10^{3}$ the change is significantly larger than $10^{-8}$ and,
therefore, it is not surprising that the fidelity increases to $0.980$ when r
is set to zero, in which case the iterations terminate when the change in log-
likelihood is less than $10^{-8}$ for all values of $\epsilon$. This
measurable increase in fidelity when r is set to zero is also observable for
$M=10^{4}$. When $M=10^{5}$ we observe that there is no change in the average
fidelity when r is set to zero, and this is explained by the fact that the
change in log-likelihood at the point where r falls below $1/M$ is of the same
order as $\Delta\mathcal{L}=10^{-8}$. For these smaller values of $M$ we
notice that decreasing $\Delta\mathcal{L}$ to $10^{-12}$, which is the lowest
we can go given the time constraint, does not affect the average fidelity as
compared to case $3$. This makes us conclude that continuing the iterations
after the change in the log-likelihood reaches below $10^{-8}$ does not
increase the performance of the algorithm for $M\leq 10^{5}$. These results
make sense statistically, since an increase in the number of samples is
accompanied by a decrease in statistical noise and a decrease in the variance
of the maximum likelihood estimate.

When $M\geq 10^{6}$ we find that on average the best reconstruction occurs in
case $2$, but the reconstruction accuracy is only slightly worse in case $4$.
This is likely due to the number of trials being small, as the case $4$ is
designed to iterate for a long time, and is thus expected to achieve the best
reconstruction accuracy on average. We conclude that, for $M\geq 10^{6}$,
using the regulating conditions of case $2$ results in reconstruction accuracy
that is indistinguishable from the accuracy that results from letting the
iterations run for a long time. The results of the reconstructions in cases
$1$ and $3$ are also of similar accuracy, which we ascribe to the fact that
the iterations always terminate without satisfying the stopping condition
$\texttt{r}_{k}\leq 1/M$. Another observation is that when allowing the
simulations to run until $\texttt{r}_{k}$ falls below $1/M$ (case $2$), the
change in the log-likelihood calculated at the step of termination decreases
significantly when $M$ is increased by orders of magnitude. In particular, for
$M=10^{7}$, when the largest average fidelity is obtained, the iterations
increase the likelihood only by $5.3*10^{-13}$, which is only half of the
$\Delta\mathcal{L}=10^{-12}$ used in case $4$. The two simulations in case $2$
for this value of $M$ took more than an hour, and since the reconstruction
accuracy is not significantly different than in case $4$, we conclude that for
near-optimal performance, given the time constraint, we can use
$\Delta\mathcal{L}=10^{-12}$ together with $\texttt{r}=1/M$ in our further
simulations with $M\geq 10^{6}$.

So far we have fixed the input state to observe the effect of the regulating
parameters of the algorithm on the quality of the reconstructed state. We have
only done such studies with two low-magnitude coherent states due to the slow
speed of the reconstruction. But we do not have a reason to expect that the
conclusions we have drawn above do not extend to arbitrary low-magnitude
coherent states since our observations are in agreement with statistical
intuition. Thus, here we only use $\Delta\mathcal{L}=10^{-12}$ and
$\texttt{r}=1/M$, and aim to observe the reconstruction accuracy for different
coherent states with the same set of probe amplitudes. We draw $M=10^{7}$
samples once for each coherent state, perform the reconstruction and calculate
$\mathcal{F}(\hat{\rho},\rho_{true})$. We show the results for three different
$\Gamma$, each having $N_{p}=6$ probe amplitudes. Fig. 5.1 shows the results
of the simulations with
$\Gamma=\\{0.3,0.6e^{i\pi/6},0.9e^{i2\pi/6},1.2e^{3\pi/6},1.5e^{i4\pi/6},1.8e^{i5\pi/6}\\}$,
while in Fig. 5.2 we used probe amplitudes randomly drawn from the $2D$ grid
spanned by magnitudes from the range $[0.3,3]$ and phases from the range
$[0,\pi]$. It can be seen that in both cases the fidelity is greater than
$0.9988$ for input states with magnitudes ranging from $0.1$ to $0.7$ and
phases spanning $[0,\pi]$. It turns out that for $N=5$ there exist IC sets of
POVMs where the magnitudes of all probes are the same. We used such an IC set
of POVMs composed of $N_{p}=6$ probes with magnitudes $0.5$ and randomly
chosen phases from $[0,\pi]$. The results are shown in Fig. 5.3, and it can be
seen that fixing the magnitude of the probe does not seem to observably affect
the accuracy of the reconstruction, and the lowest fidelity we observe is
$\approx 0.9989$. We performed these simulations for several other randomly
chosen sets of probe amplitudes and observed similar values for the fidelities
in the same range of input state magnitudes and phases.

Figure 5.1: Fidelity for different single-mode coherent input states. The
horizontal (vertical) axis denotes the magnitude (phase) of the input state.
The set of probe amplitudes composing the measurement context is
$\Gamma=\\{0.3,0.6e^{i\pi/6},0.9e^{i2\pi/6},1.2e^{3\pi/6},$
$1.5e^{i4\pi/6},1.8e^{i5\pi/6}\\}$. $M=10^{7}$ samples were drawn for each
estimate.  Figure 5.2: Fidelity for different single-mode coherent input
states. The horizontal (vertical) axis denotes the magnitude (phase) of the
input state. The set of probe amplitudes composing the measurement context is
$\Gamma=\\{3e^{i0.79\pi},1.9e^{i0.5\pi},1.4e^{i0.6\pi},$
$1.7e^{i0.55\pi},1.5e^{i0.31\pi},2.9e^{i0.63\pi}\\}$. $M=10^{7}$ samples were
drawn for each estimate.  Figure 5.3: Fidelity for different single-mode
coherent input states. The horizontal (vertical) axis denotes the magnitude
(phase) of the input state. There are $N_{p}=6$ probe states composing the
measurement context, each having magnitude $0.5$ and a phase randomly drawn
from $[0,\pi]$. $M=10^{7}$ samples were drawn for each estimate.

####  5.3.3 Analysis With Other Input States and Non-Ideal Counters

The results of the simulations with single mode coherent states give us
confidence that we have correctly implemented the numerical model and the
reconstruction algorithm. Further, the high values of state fidelity achieved
for different coherent states when the sample size is large ($\geq 10^{6}$),
and with different sets of probe amplitudes, indicate that the measurement
configuration generally allows for constructing good statistical models. By
“good” we mean that the relationship between the parameter space and the
statistical distribution of the data is such that one is very likely to end up
with a likelihood function that has high Fisher information and is maximized
near the true state for a given dataset. Similarly, the fact that these high
values of state fidelity are achieved under regulating conditions that allow
for large number of iterations during the likelihood maximization process
implies that the specific algorithm we chose to use performs well in finding
the maximum of the likelihood with our parametrization of the state space and
of the POVMs. One can argue that these simulations, for a given set of probe
amplitudes and sample size, were done with a small subset of the corresponding
statistical model that is characterized by only two separate parameters (the
magnitude and the phase of the coherent state), and therefore do not say much
about how “good” the statistical model is or how well the reconstruction
algorithm will perform for a dataset chosen at random from the statistical
model. But notice that we did not encode any assumptions about the parameter
space during the reconstruction process - the parameter space over which the
likelihood maximization was performed was the whole of $\mathcal{R}_{2,N}$.
Nevertheless, to increase our confidence in our approach, we performed
simulations with different kinds of input states described in Subsec. 5.3.1.

We considered the truncated TMSV state with $r=0.5$ and $N=8$, such that the
fidelity with the corresponding TMSV state is $\approx 0.9996$. In this case
there exist IC sets of POVMs with $N_{P}=10$ probe amplitudes. We obtained one
such set of POVMs by randomly drawing from the magnitudes in the range
$[0.3,3]$ and from the phases in the range $[0,\pi]$. More specifically,

$\displaystyle\Gamma$
$\displaystyle=\\{2.4e^{i0.31\pi},1.5e^{i0.57\pi},1.9e^{i0.03\pi},1.1e^{i0.07\pi},2.8e^{i0.59\pi},2.6e^{i0.98\pi},$
$\displaystyle
0.9e^{i0.02\pi},2.3e^{i0.78\pi},0.9e^{i0.86\pi},1.3e^{i0.35\pi}\\}.$

We drew $M=10^{7}$ samples and used the regulating parameters
$\Delta\mathcal{L}=10^{-12}$ and $\texttt{r}=1/M$, and obtained a fidelity of
$\approx 0.996$ between the true and the estimated state. The reconstruction
took more than $4$ hours and terminated when $\texttt{r}_{k}$ reached below r.
The slow reconstruction speed and the lower value of fidelity achieved than in
the case of single-mode coherent states is most likely due to the fact that
the vast majority of the entries of the true twirled state are zero. Another
reason is the larger value of $N$ which significantly increases the parameter
space, decreases the Fisher information for the same $M$, and significantly
increases the computation time for a single iteration.

We performed another set of simulations with the truncated (at $N=5$) twirled
cat state introduced in Eq. 5.3.7. We used the minimum number of different
probe states that seem to be necessary for an IC set of POVMs, $N_{p}=6$. We
again used $\texttt{r}=1/M$ and $M=10^{7}$. For $\alpha=0.5e^{i\pi/4}$, chosen
for no particular reason, and for
$\Gamma=\\{2.8e^{i0.71\pi},2.3e^{i0.48\pi},1.3e^{i0.25\pi},2e^{i0.06\pi},2.3e^{i0.23\pi},1.6e^{i0.33\pi}\\}$
we achieved a fidelity of $\approx 0.997$ when $\Delta\mathcal{L}=0$. For the
same $\alpha$, using

$\Gamma=\\{2.5e^{i0.89\pi},0.6e^{i0.23\pi},1.8e^{i0.2\pi},0.7e^{i0.21\pi},0.5e^{i0.19\pi},e^{i0.85\pi}\\}$

and $\Delta\mathcal{L}=10^{-12}$ we achieved a fidelity of $\approx 0.9990$.
Both simulations terminated with $\texttt{r}_{k}$ falling below $1/M$. The
difference in the achieved fidelities for the two different $\Gamma$ implies
that there is an observable sensitivity to the choice of $\Gamma$. To see how
the size of $\Gamma$ affects the reconstruction accuracy, we simulated for the
same $\alpha$ but used $N_{p}=10$ and $N_{p}=15$ probe states, each drawn
randomly from the magnitudes in the range $[0.3,3]$ and from the phases in the
range $[0,\pi]$. The fidelities obtained were $\approx 0.9996$ for $N_{p}=10$
and $\approx 0.9997$ for $N_{p}=15$. In both cases the iterations terminated
when the $\epsilon$ reached $10^{-30}$. The implication seems to be that for
the same r the neighborhood of the likelihood function around its maximum is
closer to (the twirled version of) $\rho_{true}$ when more probe amplitudes
are used. We performed more simulations with several other values of $\alpha$
using $N_{p}=6$ randomly selected probe states and observed similar
reconstruction accuracy, with the fidelities averaging around $0.998$.

Finally, we briefly considered simulations with lossy counters described by
the POVM elements in Eq. 5.2.2. The BS is again balanced, and both counters
are utilized. To estimate the effect of the efficiency $\nu$ (which is the
same for both counters) on the reconstruction accuracy, we performed
simulations for different values of $\nu$ in the range $[0.5,1]$ for the same
truncated ($N=5$) single-mode coherent state as an input state, and for the
same $\Gamma$ using $M=10^{6}$ samples. The regulating parameters
$\Delta\mathcal{L}=10^{-12}$ and $\texttt{r}=1/M$ were used. The fidelities
observed are shown in Table 5.3 for $4$ different single-mode coherent states
and $\Gamma$ composed of $N_{p}=6$ probe states. The members of the latter
were randomly drawn from the magnitudes in the range $[0.3,3]$ and from the
phases in the range $[0,\pi]$. The simulations where the observed fidelities
are much closer to $1$ compared to the rest (in particular, for $\nu=0.9$ for
the input state $\alpha=0.5e^{i\pi/4}$, where the fidelity has $5$ nines after
the decimal point) terminated when $\epsilon$ reached $10^{-30}$. But these
are not the only cases where r was not reached (though in general the
fidelities observed are higher when the algorithm runs until $\epsilon$
reaches $10^{-30}$), so we cannot explain the significantly higher fidelities
for these few cases by the fact that the number of iterations before
termination was several times higher. The most plausible explanation is that
the likelihood function for the particular sample peaks at a $\rho_{max}$ that
is much closer to $\rho_{true}$. Otherwise, it seems the efficiency of the
counters does not impact the accuracy of the reconstruction significantly.

Efficiency |  $\alpha=0.5e^{i\pi/4}$ and $\Gamma=$ $\displaystyle\\{2.8e^{i0.31\pi},2.8e^{i0.06\pi},0.4e^{i0.76\pi},$ $\displaystyle 1.9e^{i0.02\pi},2.1e^{i0.77\pi},0.3e^{i0.48\pi}\\}$ |  $\alpha=0.9e^{i\pi/10}$ and $\Gamma=$ $\displaystyle\\{2.6e^{i0.76\pi},2.5e^{i0.83\pi},2.1e^{i0.33\pi},$ $\displaystyle 0.5e^{i0.34\pi},2.4e^{i0.3\pi},1.1e^{i0.89\pi}\\}$
---|---|---
1 |  0.998890 |  0.998511
0.9 |  0.999997 |  0.999826
0.8 |  0.999294 |  0.999000
0.7 |  0.999621 |  0.998930
0.6 |  0.999304 |  0.998369
0.5 |  0.999417 |  0.998334
Table 5.3: The fidelity between the estimate and the true state for different values of the efficiency $\nu$, and for different single-mode coherent states and probe amplitudes. $M=10^{7}$, $\Delta\mathcal{L}=10^{-12}$ and $\texttt{r}=1/M$. The fidelities reported are from a single estimate. Efficiency |  $\alpha=-0.7$ and $\Gamma=$ $\displaystyle\\{1.9e^{i0.97\pi},2.6e^{i0.55\pi},2.2e^{i0.87\pi},$ $\displaystyle 2.4e^{i0.31\pi},1.7e^{i0.77\pi},2.8e^{i0.53\pi}\\}$ |  $\alpha=0.5e^{i5\pi/6}$ and $\Gamma=$ $\displaystyle\\{2.6e^{i0.87\pi},1.3e^{i0.32\pi},1.0e^{i0.45\pi},$ $\displaystyle 1.8e^{i0.51\pi},1.6e^{i0.02\pi},2.1e^{i0.8\pi}\\}$
---|---|---
1 |  0.998296 |  0.999362
0.9 |  0.997922 |  0.999749
0.8 |  0.999325 |  0.999957
0.7 |  0.998457 |  0.998770
0.6 |  0.998462 |  0.999625
0.5 |  0.998282 |  0.999751
Table 5.3: The fidelity between the estimate and the true state for different
values of the efficiency $\nu$, and for different single-mode coherent states
and probe amplitudes. $M=10^{7}$, $\Delta\mathcal{L}=10^{-12}$ and
$\texttt{r}=1/M$. The fidelities reported are from a single estimate.
(continued)

To test the last hypothesis, we performed $50$ simulations for the single-mode
coherent state with $\alpha=0.5e^{i\pi/4}$ and with
$\Gamma=\\{2.8e^{i0.31\pi},2.8e^{i0.06\pi},0.4e^{i0.76\pi},1.9e^{i0.02\pi},2.1e^{i0.77\pi},0.3e^{i0.48\pi}\\}$
at the efficiency $\nu=0.9$. We want to observe the variations in the fidelity
depending on the drawn sample, in order to determine whether the very high
fidelity observed above is due to the value of $\nu$ or due to the particular
sample that was drawn in that trial. We obtained an average fidelity of
$0.999718$ with a standard deviation of $0.000100$. The histogram is shown in
Fig. 5.4. It can be seen that the estimates are clustered around the fidelity
of $0.9997$, which is near the same order of accuracy as the results for the
other input states in Table 5.3.

Figure 5.4: Histogram of the fidelity for the single-mode coherent state
$\alpha=0.5e^{i\pi/4}$ and for
$\Gamma=\\{2.8e^{i0.31\pi},2.8e^{i0.06\pi},0.4e^{i0.76\pi},1.9e^{i0.02\pi},2.1e^{i0.77\pi},0.3e^{i0.48\pi}\\}$.
Efficiency of the counters is $\nu=0.9$. $50$ simulated reconstructions are
shown, where $M=10^{7}$ for each simulation.

##  Chapter 6 Analysis of the Experiment

In this chapter we describe the experiment that was performed by David S.
Phillips, Thomas Gerrits and Michael Mazurek at NIST, Boulder, and our
subsequent data analysis. The data analysis was abandoned midway after some
systemic discrepancies were discovered in the data that we could not trace
back to the conditions of the experiment. Nevertheless, these systemic
discrepancies are not large enough to prevent what one might consider a good
agreement between the reconstructed state and the (assumed) true state.

The experiment was initially meant to model the WFH configuration where the BS
is characterized by a trivial partition and thus the input state could be
described as a two-mode twirled state in $\mathcal{R}_{2,N}$. A pulsed laser
was used to simultaneously produce the LO and input states. Only coherent
input states were generated. The mode of the LO was set to a fixed
polarization. The polarization of the input state was controlled with respect
to the polarization of the LO. The mode of the input state matching the LO
(mode $1$) is defined by having the same polarization and temporal structure
as the LO. All occupied modes in the orthogonal polarization, as well as the
modes in the same polarization as the LO but with zero temporal overlap with
it, are orthogonal to mode $1$. According to our results at the end of Sec.
3.3 we could assume only one of these modes is occupied and call it mode $2$.
However, after the experiment was conducted, it was found that the BS acts
differently depending on the polarization. The WFH configuration that models
the experiment belongs to the family with a BS characterized by a partition
size $K=2$, and, therefore, we need to model the input state as living in
$\mathcal{R}_{3,N}$. Mode $1$ is still the same, but the set of modes in the
same polarization as the LO but with orthogonal temporal shapes are
distinguishable from the set of modes in the orthogonal polarization to the
LO, as the BS acts differently on these two sets of modes. According to our
theoretical results, we can assume only one mode in the former set of modes is
occupied, and similarly for the latter set of modes. These two modes are
denoted as mode $2$ and mode $3$, respectively.

We start by giving a summary of the experiment in Sec. 6.1 (a more detailed
description can be found in David S. Phillips’ thesis [35, Chap. 7]). The
following two sections are devoted to data analysis. Sec. 6.2 describes how
the data for reconstruction is generated - namely, how the relative
frequencies corresponding to the outcomes of a chosen measurement context are
obtained. It also discusses how the “best guess” of the input state is
obtained. Sec. 6.3 describes the results from performing maximum likelihood
reconstruction with this data.

###  6.1 Description of the Experiment

We divide our description of the experimental setup into four parts: (1)
generation of the input and the LO states, (2) manipulation of the phase and
the magnitude of the LO, (3) implementation of the BS action, and (4)
measurement by the counters. The schematic of the optical setup is shown in
Fig. 6.1. A pulsed laser with a wavelength of $1550\textrm{ nm}$ and pulse
length of $100\textrm{ ps}$ was used to generate both the LO and the input
state. A function generator was used to trigger the laser pulses at
$100\textrm{ kHz}$, resulting in $10^{5}$ measurement shots per second. This
number is chosen to give enough time for the photon counters to reset after
each shot. The laser light was first passed through a polarization-maintaining
(PM) and single-mode (SM) fiber-optic cable to shape it into a well-defined
mode. A variable attenuator was used to decrease the magnitude of the coherent
pulse to having an average of at most several photons. This caused some de-
polarization, and thus a fiber polarizer (FP) was used before a $90:10$ fiber
beam splitter (FBS). The $90\%$ arm of the FBS was used to generate the LO
state (the upper arm in Fig. 6.1), while the $10\%$ arm was used for
constructing the input state (the lower arm in Fig. 6.1).

Figure 6.1: Schematic of the experimental setup (taken from David S. Phillips’
thesis [35, Fig. 7.2]). Some of the abbreviations of the components, starting
from the left, are as follows: TESs stand for transition-edge sensors, FBS
stands for a fiber-optic BS, HWP stands for a half-wave plate, PBS stands for
a polarizing BS, and FP stands for a fiber-optic polarizer. The description is
in the main text, and further information can be found in Ref. [35].

To generate the input state, the corresponding output of the FBS was first
passed through a formation of four mirrors to control the path length and the
beam alignment. Afterwards the pulse was passed through a polarizing beam
splitter (PBS), which separates the pulse into orthogonal polarization
components. We refer to the polarization that is passed through to the next
stage of the setup by “vertical”, and this is aligned with the polarization of
the FP. The orthogonal “horizontal” component is discarded. The vertically
polarized input is then passed through a half-wave plate (HWP) to set its
polarization with respect to the vertical (different polarizations correspond
to different input states). To construct the LO state, the $90\%$ arm of the
FBS is passed through an arrangement of mirrors, two of which are attached to
a stage controlled by a piezoelectric device, which is in turn attached to a
more coarsely controllable translation stage. This arrangement is used to
alter the path length of the LO with a precision that is of the order of a nm,
and thus allows for an adjustment of the relative phase between the input and
the LO states. We describe the intricacies of this arrangement in more detail
in the next section. After phase adjustment, the LO state is passed through a
sequence composed of a PBS, a HWP and a polarizer to adjust the magnitude of
the pulse. The polarizer is set to vertical polarization so that the
horizontal polarization component of the LO is blocked from passing through.
The magnitude of the pulse in the vertical polarization, and thus the
magnitude of the pulse after the polarizer, can be controlled by tuning the
angle of the HWP. The input state as well as the LO state in the model are
assumed to be the states in the input fibers of the FBS right before
interaction.

The next part of the setup is the BS, which is realized by a PM FBS that was
designed to act as a $50:50$ BS on the vertically polarized modes.
Measurements by power meters at the outputs of the FBS without using the
variable attenuator and directing the light in a vertical polarization to one
of the input ports showed that the splitting ratio of the BS on the vertical
modes was $49.46:50.54$. The same measurement with horizontally polarized
light gave a splitting ratio of $73.58:26.42$. We used these numbers in our
model, but the error in them is estimated to be up to $1\%$. The pulses at the
outputs of the FBS are then sent to two transition-edge sensors (TESs), which
we describe below. The fibers leading to the TESs were not PM, and the
experimenters used mating couplers to connect them with the outputs of the
FBS. The resulting losses from the couplings were measured and incorporated
into the POVM elements corresponding to the outcomes of the TESs (described in
Subsec. 6.2.1).

The TESs operate by keeping a film of superconducting material at the critical
temperature and bias current - that is, the temperature and current at which
the transition occurs from the superconducting to the conducting regime. This
transition is accompanied by a rapid increase of the resistance from a value
of zero to a finite value, which then changes very little upon further
increase of the temperature. Due to the shape of the dependence of the
resistance on the temperature around the critical point, very small increases
in the temperature can be detected by monitoring the resistance across the
superconducting film. This is precisely what happens when one or several
photons are absorbed by the film. The critical temperature of a TES is
typically around $100\textrm{ mK}$, and an adiabatic demagnetisation
refrigerator is used to cool it to such a low temperature and provide good
temperature control. The change in the resistance across the superconducting
film is detected by measuring the change in the magnetic field generated by an
inductor that is in series with the TES and through which a constant current
is maintained. The measurement of the magnetic field is performed by an array
of superconducting quantum interference devices (SQUIDs), which are held at a
temperature of about $10\textrm{ K}$. The design of a TES is its own field of
study, and we will not go into further detail here. More information about the
design of the circuits incorporating the TES, the signal readout and its
analysis can be found in [14]. For the purposes of this thesis, it suffices to
mention that the experimenters provided us with a measurement of a quantity
proportional to the energy absorbed by the TES from the light pulse. We
describe our analysis with these values in the next section.

The experiment lasted five days, with a different input coherent state
$\ket{\alpha}$ prepared and measured for each day. An effort was made to keep
the magnitudes of the input states the same at
$\absolutevalue{\alpha}\approx\sqrt{0.7}$, and the states prepared on
different days differed by the angle of the polarization with respect to the
vertical polarization as set by the HWP. The HWP angles used were
$54.8^{\circ}$, $64.80^{\circ}$, $76.74^{\circ}$, $84.59^{\circ}$ and
$90.66^{\circ}$. The relationship between the HWP angle ($\theta_{HWP}$) and
the polarization angle with respect to the vertical ($\theta_{p}$) was
measured to be $\theta_{p}=2\theta_{HWP}-109.6^{\circ}$ by using power meters
instead of TESs and bypassing the attenuator. The chosen HWP angles were meant
to test the whole polarization range, with $\theta_{HWP}=54.8^{\circ}$
corresponding to vertical polarization and $\theta_{HWP}=90.66^{\circ}$
corresponding to a polarization angle of $71.72^{\circ}$ with respect to the
LO polarization. At the beginning of each day, the input state was prepared
and measured by the experimental setup by blocking the LO arm, thus providing
us with TES data to obtain a more accurate estimate of the input state to
compare the reconstructed state with. We refer to these datasets as “input
state only” datasets.

Three different LO magnitudes were prepared for each input state, and
individual measurements of these (by blocking the input state from entering
the FBS) were made. We thus obtained TES data to estimate the magnitudes of
the LO states. We refer to these datasets as “LO state only” datasets, while
the measurements performed corresponding to the WFH configuration are referred
to as “interference experiments”. In the interference experiment, the
piezoelectric stage was used to sweep the position of the holding mirrors
(shown in Fig. 6.1) repeatedly. More precisely, the position of the
piezoelectric stage is linearly related to the voltage supplied to the
piezoelectric device, and a saw-tooth voltage function was fed to the device.
The period of the function was set to $\approx 0.8$ seconds, and the height of
the saw-tooth with respect to the beginning point of the sweep was chosen such
that the resulting voltage difference roughly corresponded to $8.4$ radians of
phase change of the LO pulse (thus giving us data from the full range of the
relative phase). TES data was collected in $0.8$ second intervals, which
corresponded to $83856$ shots, with each interval followed by $0.2$ seconds of
no measurements, until around $3\cdot 10^{8}$ measurement instances were
collected for each LO magnitude. We refer to a dataset taken within a single
interval as a single “piezo ramp” dataset, even though that is not strictly
accurate, as the period of the saw-tooth function is slightly different from
the acquisition window. The fast speed of the ramp of the piezoelectric stage
was chosen to mitigate the effect of the drift of the stage. In particular,
the experimenters found that the stage can drift significantly (in terms of
its impact on the relative phase) within dozens of seconds at the same
voltage, and this made it impossible to rely on the voltage input to estimate
the relative phase. Therefore, the data was divided into consecutive pairs of
piezo ramps, with the $83856$ data points from the first ramp used to obtain
the relationship between the voltage and phase, and the second set of $83856$
data points used for reconstruction. Therefore, the phase stability required
was of the order of $2$ seconds, which is achievable by the measurement setup.

###  6.2 Generating Data for Reconstruction

In this section we describe how the “true” input and the “true” LO states for
the different interference experiments were estimated from the “input state
only” and “LO state only” datasets, respectively. We also describe how the
different POVMs for the different TES measurement outcomes at different LO
amplitudes were constructed, as well as how the relative frequencies
corresponding to the measurements by different POVM elements were obtained
from the interference experiments. The relative frequencies were then used for
state reconstruction analysis (described in the next section).

####  6.2.1 Construction of the POVMs

The construction of the POVMs corresponding to the measurement outcomes of the
TESs as part of the measurement configuration of the experiment is a multistep
process. The first stage was performing tomography of the TES detectors. This
was done by the experimenters according to the method described in Ref. [6] on
a different day before the start of the experiment. Their method uses a set of
known coherent states as inputs to the TES. As mentioned in the last section,
the analysis of the TES output for a given shot results in a single number,
which we refer to as the matched filter value (MFV), and which contains
information about the number of photons absorbed by the TES from the light
pulse. As an example, a histogram of MFVs of one of the TESs for one of the
“LO state only” datasets is shown in Fig. 6.2. It can be clearly seen that the
MFVs are naturally grouped according to photon number, and that there are
overlaps between the groupings. Broadly speaking, there are two kinds of noise
that affect the estimation of the photon number measurement of the state from
the MFV outcome. The first kind of noise has to do with all sorts of processes
occurring in the complicated circuitry that links the act of absorption of the
light pulse to the final readout. Using the example above, if one were to
assume that the state absorbed is a Fock state with a well-defined frequency,
the corresponding distribution of the MFVs would be one of the Gaussian-
looking profiles in Fig. 6.2. Because of the overlap between the different
profiles, this source of noise would produce an uncertainty in the photon
number of the Fock state. The second kind of noise has to do with losses and
background light present in the absorption process. This includes losses in
the coupling with the fiber leading to the TES. In this case, even if the
Gaussian profiles had near-zero overlap, so that the number of photons
absorbed could be deduced with high certainty, there would still be
uncertainty in the photon number measurement on the input state. Thus, the
tomography of the TESs needs to take these two kinds of sources of noise into
account and link the MFVs to the projectors onto photon numbers on the input
state.

Figure 6.2: Example of a histogram of MFVs for many shots. The dataset shown
is obtained from one of the TESs during one of the “LO state only”
measurements for the input state with $\theta_{HWP}=54.8^{\circ}$.

The method of detector tomography that was used is based on assigning photon
numbers to MFVs. In particular, the method produces a set of MFV intervals
corresponding to each photon number. Using Fig. 6.2 for illustration again,
the edges of the intervals (we refer to these as “bounds”) would be near the
troughs in the histogram. While more advanced data analysis can be imagined
given the knowledge of the TES design, this grouping strategy requires minimal
processing and is based on the statistically sound assumption that the
absorption of a given number of photons produces an MFV with Gaussian error.
The TES detector tomography also produces a $10\times 26$ matrix which encodes
the conditional probabilities $P(n|m)$ of detecting $n=0,\ldots,8$ or $n>8$
photons when $m=0,\ldots,24$ or $m>24$ photons are present in the fiber
leading to the TES. We refer to these matrices as $T_{1}$ and $T_{2}$ for the
two TESs. There are additional losses in the couplings between the outputs of
the FBS and the fibers leading to the TESs, which were estimated by the
experimenters. The losses are modeled as an additional Markov process as
described in Subsec. 5.2.2, to obtain the conditional probabilities for each
TES of detecting a particular photon number when a certain number of photons
is present at the corresponding output port of the FBS. More specifically, let
us denote the updated conditional probabilities of detecting $n=0,\ldots,8$ or
$n>8$ photons when $m=0,\ldots,24$ or $m>24$ photons are present at the
outputs of the FBS by $T^{\prime}_{1}$ and $T^{\prime}_{2}$ for the two TESs.
These were approximated as

$\displaystyle
T^{\prime}_{i,nm}=\sum_{j=0}^{m}T_{i,nj}\binom{m}{j}\nu_{i}^{j}(1-\nu_{i})^{m-j},$
(6.2.1)

where $i=1$ or $i=2$, $\nu_{i}$ is the transmission coefficient for the $i$’th
TES, and the approximation assumes that there are no more than $25$ photons in
the outputs of the FBS. This approximation is justified given the low
magnitudes of the LO and input states. The transmission coefficients were
measured by the experimenters to be $\nu_{1}=0.9586$ and $\nu_{2}=0.9413$.

Then, for the $i$’th TES, the positive operator on the state at the matching
output port of the FBS corresponding to the measurement of $k$ photons is
given by

$\displaystyle\hat{D^{\prime\prime}}_{i,k}=\sum_{m=0}^{25}T^{\prime}_{i,km}\hat{D}_{i,m}.$
(6.2.2)

The final POVM elements corresponding to $k$ and $l$ photon number
measurements by the two TESs for a given LO amplitude are calculated similarly
to Eq. 5.2.2. Let us denote these by $\Pi^{\prime\prime}_{kl}$. Skipping some
of the steps in the derivation that are already familiar from Subsec. 5.2.2,
we obtain

$\displaystyle\Pi^{\prime\prime}_{kl}$
$\displaystyle=\tr_{\textsl{p}}\left(\left(I\otimes\ket{\gamma}\bra{\gamma}\otimes\ket{\vec{0}}\bra{\vec{0}}\right)U^{\dagger}(\hat{D^{\prime\prime}}_{1,k}\otimes\hat{D^{\prime\prime}}_{2,l})U\right)$
$\displaystyle=\sum_{m=0}^{25}\sum_{n=0}^{25}T^{\prime}_{1,km}T^{\prime}_{2,ln}\Pi_{mn},$
(6.2.3)

where
$\Pi_{mn}=\tr_{\textsl{p}}((I\otimes\ket{\gamma}\bra{\gamma}\otimes\ket{\vec{0}}\bra{\vec{0}})U^{\dagger}\hat{D}_{kl}U)$.
As a reminder, since our BS is characterized by the partition size $K=2$, the
unitary $U$ is acting on $3$ modes. We encode its action in mode space as

$\displaystyle B=\begin{pmatrix}\sqrt{0.4946}&\sqrt{0.5054}\\\
-\sqrt{0.5054}&\sqrt{0.4946}\end{pmatrix}\oplus\begin{pmatrix}\sqrt{0.4946}&\sqrt{0.5054}\\\
-\sqrt{0.5054}&\sqrt{0.4946}\end{pmatrix}\oplus\begin{pmatrix}\sqrt{0.2642}&\sqrt{0.7358}\\\
-\sqrt{0.7358}&\sqrt{0.2642}\end{pmatrix}$ (6.2.4)

using the experimentally measured splitting ratios for the “vertical” and
“horizontal” polarizations. For a given LO amplitude $\gamma$, these
$\Pi_{kl}$ are calculated similarly to Eq. 5.2.1 and are encoded as elements
of $\mathcal{R}_{3,N}^{s}$. We do not present this calculation here.

The circuit by which the resistance across the TES is measured utilizes a
SQUID array, which is cooled anew each day. This reset of the cooling chamber
causes a slight shift in the MFVs by an overall scaling factor. More
precisely, the output of the SQUID indirectly measures the change in the
resistance of the TES, and we assume that the effect of the reset on this
measurement is small enough that it can be assumed to be linear to first
order. This necessitates a correction in the bounds for the MFVs by an overall
scaling factor. We find that a correction in the bounds throughout the day is
also desirable. In particular, the interference experiment with a particular
LO magnitude is performed right after the corresponding “LO state only”
measurement. The latter can be used to estimate the correct scaling factor, as
described in the next paragraph, which is then used to assign photon numbers
to the TES data for the corresponding interference experiment. The “input
state only” dataset is also used to estimate the correct scaling factor at the
time of its measurement, which is only used to gain an estimate of the true
state to compare the reconstruction of the twirled state with.

We perform model selection to obtain the estimates of the scaling factors for
the bounds of both TESs. More specifically, since we know that both the LO and
the input state are prepared in coherent states, the photon number
distributions at the outputs of the FBS are Poissonian with a fixed Poisson
parameter (given by the square of the coherent state magnitude) for each
output for each of the four datasets (that is, the three “LO state only” and
the “input state only” datasets). The probability of measuring $m$ photons in
output port $i$ is determined by $\hat{D}_{i,m}$. We have access to the
$\hat{D^{\prime\prime}}_{i,k}$, which we use to relate the relative
frequencies for the output of each TES to the corresponding Poisson parameter.
We refer to this statistical model as the “Poissonian model”. The optimization
is based on calculating likelihood ratios for different scaling factors. In
particular, a given scaling factor gives different photon number frequencies
because the scaling factor influences the MFV intervals corresponding to each
photon number. We use the Poissonian model to estimate the most likely Poisson
parameter for these frequencies and calculate the value of the likelihood
function for this parameter. We do this for each of the four datasets. We then
consider the ratio of this value of the likelihood function over the value of
the likelihood function of the unconstrained model (that is, when the
frequencies are taken as the probabilities). We maximize this likelihood ratio
over a range of scaling factors, and choose the scaling factor that gives the
largest likelihood ratio for each dataset.

To estimate the accuracy of the Poissonian model when the optimal scaling
factors are used, we performed parametric bootstrapping on the “LO state only”
datasets corresponding to the input state with $\theta_{HWP}=76.74^{\circ}$.
More specifically, for each of the three datasets, we computed the coherent
state magnitude that maximizes the likelihood function and calculated the
corresponding probabilities for each TES. We then simulated random samples
from these probabilities the same number of times as the original number of
experimental data points (which was of the order of $10^{6}$). We did this
$1000$ times for each LO state. After generating the likelihood ratios for the
bootstrapped distributions (that is, the ratio of the maximized likelihood
function according to the Poissionian model over the unconstrained likelihood
function) we looked at the deviation of the original likelihood ratio from the
bootstrapped likelihood ratio distribution. We found that overall the model
mismatch is not too bad, with the largest deviation being about $7\sigma$ for
the LO state with the largest magnitude. While this deviation might seem to be
very big, the likelihood ratio test in this case is very sensitive given the
number of data points. It is possible that this deviation can be partly
explained by the residual memory effects of the TESs. The “LO state only” data
is also taken at $100\textrm{ kHz}$, and we found that $10\textrm{
}\mu\textrm{s}$ is not enough time to completely eliminate the effects from
the preceding absorption event.

The magnitude-squares of the LO states that were incorporated into the POVM
elements were calculated from the values that maximized the likelihoods of the
Poissonian models using the corresponding “LO state only” dataset for each
TES. More specifically, the magnitude-square of the LO state is assumed to be
equal to the sum of the magnitude-squares at the output ports of the FBS, and
the latter are estimated from the corresponding TES data by the procedure
described above. In the next subsection we describe how the relative phases
between the LO and input states were estimated from the “piezo ramp” datasets.

####  6.2.2 Estimation of the Relative Phase and Generation of Frequencies

A segment of the voltage function fed to the piezoelectric device is shown in
Fig. 6.2.2. The vertical axis is the voltage measured across the device, while
the horizontal axis is time (more specifically, the “experimental time” taken
by the clock that stops ticking when no measurements are made). One can
observe the existence of fast transients (a zoomed-in segment is shown in part
(b) of the figure), which occur when the voltage fed by the function generator
starts a new saw-tooth cycle by making a sudden jump from the largest value to
the lowest. The piezo cannot change its position instantaneously, and this
results in the recorded transient response. One can also notice that there are
instantaneous jumps in the voltage, occurring at roughly the same intervals as
the transients. These are due to the $0.2$ seconds of no data taking within
each second during which the piezoelectric device continuous its cycle. More
specifically, data of the voltage readout was taken continuously, but the part
of this data when no measurements were made was edited out. The data taken
between two such jumps is what we denote as a single “piezo ramp” dataset. We
use the first “piezo ramp” dataset to estimate the relationship between the
voltage and the phase of the LO, which is assumed to be linear during the non-
transient part of the time trace. We use this estimate to assign phases to the
voltage values in the consecutive “piezo ramp” dataset, and we repeat this for
all the data we have.

Figure 6.3: (a) The trace of the voltage measured across the piezoelectric
device over several seconds. The second figure (b) shows a zoomed-in segment
of (a) for a better view of the transient response of the device to an
instantaneous shift in the voltage input.

.

The relationship between the relative phase $\phi$ and the voltage readout $V$
is assumed to be linear and modeled by the equation $\phi=a\cdot V+b$. The
estimation of the constants $a$ and $b$ for a given “piezo ramp” dataset is
done as follows. Since both inputs to the FBS are coherent states, we expect
the states at the outputs of the FBS to be coherent states with magnitudes
that have sinusoidal dependence on the relative phase (more explicit
description is provided in the next subsection). The algorithm first removes
the part of the data corresponding to the transient voltage response described
above, and assigns photon numbers to the MFVs of both TESs according to the
estimated bounds for each voltage value. Once that is done, it orders the data
of the TESs according to the corresponding voltage values, and divides the
remaining data into $241$ consecutive segments. The number of data points in
each segment averages to about $348$, which corresponds to almost $2^{\circ}$
of relative phase change. The number of data points in each segment is not
always the same, since the procedure that removes the transient part of the
piezo trace does not leave the same number of data points each time. We then
make the assumption that the magnitude of the coherent state at the outputs of
the FBS can be considered constant throughout the aforementioned $2^{\circ}$
window, so that we can use the Poissonian model described in the previous
subsection to estimate the magnitude of the coherent state at each output
using maximum likelihood estimation with each of the $348$ data points. We
also calculate the Fisher Information for each estimate, which is used to
estimate the corresponding variance by assuming the variance is equal to the
inverse of the Fisher Information. We then use the estimated coherent state
magnitudes for each segment, and perform nonlinear weighed least squares
regression, with the weights given by the inverses of the variances, to a
general sinusoid $c\cdot\sin(a\cdot V+b)+d$ for each TES. This allows us to
obtain two estimates for each of the parameters $a$ and $b$, and we take the
averages of the two for the final estimates of these parameters. An example of
a single fit for one of the TESs for the input state with
$\theta_{HWP}=76.74^{\circ}$ is shown in Fig. 6.4.

Figure 6.4: Fitting a sinusoid to the estimated coherent state magnitudes at
one of the outputs of the FBS at a range of relative phase values. The HWP
angle of the input state is $\theta_{HWP}=76.74^{\circ}$.

There are several potential issues with this procedure. First, the nonlinear
least-squares fitting of a sinusoid is not a convex optimization problem and
the algorithm can return a bad estimate. Second, the coherent state magnitudes
are not constant over a $2^{\circ}$ phase change period, which affects the
accuracy of the Poissonian model. Third, the inverse Fisher Information may
not be a very accurate estimate of the variance when the number of data points
is only about $348$, which affects the accuracy of the weights assigned to the
points in the least-squares procedure. We performed some simulations to
investigate the accuracy of the last assumption. In particular, we simulated
the maximum likelihood estimation procedure for $348$ datapoints $1000$ times
for a set of different magnitude-squares using the Poissonian model to
generate the datapoints. We found that the difference between the sample
variance (based on the $1000$ estimates) and the inverse Fisher information is
generally within $2\%$.

We keep track of the sum of the squares of the weighed residuals, which should
approximately follow a chi-squared distribution. More specifically, this sum
is a chi-squared distribution when the distributions of each of the $241$
estimates of the magnitudes are Gaussians with variances equal to the
corresponding weights. If the model were exact, we expect the average of the
sum of the squares of the weighed residuals to be equal to the number of
degrees of freedom of the model. The latter is calculated to be $241-4=237$,
where we subtracted the number $4$ to account for the $4$ parameters that
define the sinusoid. We find that, very rarely, the optimization algorithm
converges to a very bad fit. We therefore use a cutoff value of $310$ for the
sum of the squared weighed residuals to reject the fits associated with a
higher value of this quantity. The corresponding p-value is $\approx 0.1\%$.
For such fits, the next “piezo ramp” dataset is also discarded. This seems
like an unnecessary waste of data, but we find that the fits that are rejected
are very few. We found that the observed average chi-squared value had a
p-value of around $25\%$ instead of the $50\%$ one would expect from a perfect
model. We originally assumed that this was consistent with our expectations,
since there is a clear model mismatch as argued above. But, as described in
the next section, there is a more profound model mismatch that cannot be
explained with these arguments, and which eventually convinced us to abandon
the further analysis of the experiment.

At the end of the procedure described above, for a given input and LO state,
we obtain a collection of photon-number outcomes associated with phase values
for each TES. We divide the phase range into $2^{\circ}$ consecutive intervals
starting from $[0^{\circ},2^{\circ})$, and group the data of the photon-number
outcomes into these intervals according to where their associated phase values
fall. We then count the number of outcomes in each phase interval for each
TES, and divide these numbers by the total number of outcomes to obtain the
corresponding relative frequencies. For the interference experiment with the
weakest LO magnitude, we first divide the data into two parts (described in
the next subsection) and then generate the relative frequencies for each part.
This is necessary to have a separate set of data of photon-number outcomes for
estimating the true state (again, described in the next subsection). We make
the assumption that the relative LO phase at which the frequencies in a
particular phase interval were generated is equal to the value of the phase in
the middle of that interval. We do this for several reasons. First, it is
expensive to generate the POVM elements, to store them and to use them to
construct the functions in the $R\rho R$ algorithm, and it would be
practically impossible to use a separate POVM element for each phase value as
the number of samples is of the order of $10^{7}$. Second, this should not
have a large effect on the reconstruction, as we expect any two POVM elements
with the same photon-number outcomes but differing up to $2^{\circ}$ in phase
to not be too different from each other. This is confirmed by parametric
bootstrapping studies described in the next section. Perhaps the best argument
for choosing this phase interval is that we want to make use of the order of
several dozen LO amplitudes out of which to construct our measurement context,
and we want to have a reasonably large number of samples for our POVMs. In
this case, we can obtain around $2.5\cdot 10^{7}$ samples for each of the
measurement contexts we decide to construct.

####  6.2.3 Determining the “True” Input State

Here we describe how we use part of the data of photon-number outcomes from
the interference experiment with the weakest LO state, as well as the “input
state only” dataset to obtain a “best guess” of the true input state. For this
purpose, we assume that the input state is in a coherent state, and we denote
it by $\ket{\psi_{\textrm{in}}}$. The corresponding truncated (at $N=4$
maximum photon number), normalized and then twirled three-mode state is
denoted as $\rho_{true}$. The latter is needed for comparing the best guess
with the state obtained by the maximum likelihood reconstruction procedure
(described in the next section) using state fidelity as a measure.

Before we describe our estimation procedure, a set of observations and
definitions are in order that also apply to the maximum likelihood
reconstruction procedure. The first observation is that the LO and the input
state may undergo individual phase changes before their interaction in the
middle of the FBS. We thus define the input state as the state within its
corresponding input branch of the FBS right before the action of the BS. The
absolute value of the phase of the input state is not observable, and for
convenience we set it to $0$. That is, we choose our reference frame such that
the absolute phase is $0$. We similarly define the LO state as the state in
its corresponding input branch of the FBS right before the interaction with
the input state. With these definitions and choice of the reference frame, and
assuming the unitary in Eq. 6.2.4 describes the BS action, the relative phase
$\phi$ between the LO and the input state is equal to the absolute phase of
the LO.

The magnitude of the input state is estimated using the same procedure as for
the LO states. Namely, the Poissonian model is used to estimate the magnitudes
of the pulses at the outputs of the FBS from the “input state only” dataset,
and the magnitude-square of the input state is computed as the sum of the
estimated magnitude-squares at the outputs. It is left to describe how the
polarization angle and the temporal overlap with the LO were estimated.
Remember that we denoted the polarization axis of the LO by “vertical”. The
input state with magnitude $\absolutevalue{\alpha}$ can be written as

$\displaystyle\ket{\psi_{\textrm{in}}}=\ket{\absolutevalue{\alpha}\cos(\theta_{p})\cos(\theta_{t})}_{1}\otimes\ket{\absolutevalue{\alpha}\cos(\theta_{p})\sin(\theta_{t})}_{2}\otimes\ket{\absolutevalue{\alpha}\sin(\theta_{p})}_{3},$
(6.2.5)

where the subscripts of the kets denote the corresponding mode number,
$\theta_{p}$ is the polarization angle with respect to the LO polarization,
and $\theta_{t}$ signifies the temporal overlap between the modes of the LO
and the modes of the input state in the vertical polarization. With this
parametrization, the LO state is written as
$\ket{\gamma}_{1}\otimes\ket{0}_{2}\otimes\ket{0}_{3}$. We omit the tensor
products in our descriptions in the rest of this section for brevity.

After the action of the FBS, the states at the output ports of the FBS are
given by

$\displaystyle\ket{\psi_{\textrm{out},1}}=\ket{B_{1,11}\absolutevalue{\alpha}\cos(\theta_{p})\cos(\theta_{t})+B_{1,12}\gamma}_{1}\ket{B_{1,11}\absolutevalue{\alpha}\cos(\theta_{p})\sin(\theta_{t})}_{2}\ket{B_{2,11}\absolutevalue{\alpha}\sin(\theta_{p})}_{3},$
(6.2.6)

and

$\displaystyle\ket{\psi_{\textrm{out},2}}=\ket{B_{1,21}\absolutevalue{\alpha}\cos(\theta_{p})\cos(\theta_{t})+B_{1,22}\gamma}_{1}\ket{B_{1,21}\absolutevalue{\alpha}\cos(\theta_{p})\sin(\theta_{t})}_{2}\ket{B_{2,21}\absolutevalue{\alpha}\sin(\theta_{p})}_{3},$
(6.2.7)

where $B_{1}=\begin{pmatrix}\sqrt{0.4946}&\sqrt{0.5054}\\\
-\sqrt{0.5054}&\sqrt{0.4946}\end{pmatrix}$ and
$B_{2}=\begin{pmatrix}\sqrt{0.2642}&\sqrt{0.7358}\\\
-\sqrt{0.7358}&\sqrt{0.2642}\end{pmatrix}$ according to Eq. 6.2.4. For the
“input state only” datasets, the LO amplitude is zero, and the magnitudes of
$\ket{\psi_{\textrm{out},1}}$ and $\ket{\psi_{\textrm{out},2}}$ are constant.
We estimate $\theta_{p}$ from the estimates of these magnitudes (which are
obtained using the Poissonian model). In particular, we obtain a single
equation with single unknown ($\theta_{p}$) from the data of each TES, which
we solve and take their average as the estimate of our phase. These equations
are

$\displaystyle\textrm{Mag}^{2}_{1}$
$\displaystyle=\absolutevalue{B_{1,11}\absolutevalue{\alpha}\cos(\theta_{p})}^{2}+\absolutevalue{B_{2,11}\absolutevalue{\alpha}\sin(\theta_{p})}^{2}$
$\displaystyle=\absolutevalue{\alpha}^{2}(B_{1,11}^{2}\cos^{2}(\theta_{p})+B_{2,11}^{2}\sin^{2}(\theta_{p}))$
(6.2.8)

and

$\displaystyle\textrm{Mag}^{2}_{2}$
$\displaystyle=\absolutevalue{B_{1,21}\absolutevalue{\alpha}\cos(\theta_{p})}^{2}+\absolutevalue{B_{2,21}\absolutevalue{\alpha}\sin(\theta_{p})}^{2}$
$\displaystyle=\absolutevalue{\alpha}^{2}(B_{1,21}^{2}\cos^{2}(\theta_{p})+B_{2,21}^{2}\sin^{2}(\theta_{p})),$
(6.2.9)

where $\textrm{Mag}^{2}_{1}$ and $\textrm{Mag}^{2}_{2}$ stand for the squares
of the magnitudes for counters 1 and 2, respectively. Note that we also have
an estimate of $\theta_{p}$ from the measured HWP angle, but we find that
these do not match too closely, which was another factor in the decision to
terminate the data analysis. We detail this at the end of the next section.

To estimate the temporal overlap $\absolutevalue{\cos(\theta_{t})}^{2}$
between the LO and the part of the input state in the vertical polarization we
use the first $1/12$’th of the photon-number data generated from the
interference experiment with the weakest LO state to estimate $\theta_{t}$.
The rest of the data generated from all three interference experiments are
used for reconstructing the three-mode twirled state using the maximum
likelihood estimation procedure. Now, the squares of the magnitudes at the
outputs of the FBS are functions of $\theta_{t}$. In particular,

$\displaystyle\textrm{Mag}^{2}_{1}$
$\displaystyle=\absolutevalue{B_{1,11}\absolutevalue{\alpha}\cos(\theta_{p})\cos(\theta_{t})+B_{1,12}\gamma}^{2}+\absolutevalue{B_{1,11}\absolutevalue{\alpha}\cos(\theta_{p})\sin(\theta_{t})}^{2}+\absolutevalue{B_{2,11}\absolutevalue{\alpha}\sin(\theta_{p})}^{2},$
(6.2.10)

and

$\displaystyle\textrm{Mag}^{2}_{2}$
$\displaystyle=\absolutevalue{B_{1,21}\absolutevalue{\alpha}\cos(\theta_{p})\cos(\theta_{t})+B_{1,22}\gamma}^{2}+\absolutevalue{B_{1,21}\absolutevalue{\alpha}\cos(\theta_{p})\sin(\theta_{t})}^{2}+\absolutevalue{B_{2,21}\absolutevalue{\alpha}\sin(\theta_{p})}^{2}.$
(6.2.11)

Notice that the dependence on the phase of the LO is due to the appearance of
$\gamma$ in the first summands. We remind the reader that we use roughly half
of the data from the intereference experiment for each LO state to generate a
set of samples of photon-number for both TESs associated with $180$ equally
spaced different values of the relative phase. We calculate the relative
frequencies from the photon-number data obtained with this procedure from the
first $1/12$’th duration of the experiment. We use the Poissonian model on
each set of frequencies to obtain an estimate of the magnitude of the coherent
state at both outputs of the FBS for all phase points. We then perform
nonlinear least squares optimization to fit the equations above to the squares
of the magnitudes, and obtain two estimates of $\theta_{t}$. We take their
average as our final estimate.

We should note that the procedure for estimating $\rho_{true}$ is fairly
involved, uses multiple datasets, relies on the accuracy of the assumed BS
coefficients and on detector tomography, and on the procedure for the
estimation of the relative phase corresponding to the photon-number outcomes.
Thus, our “best guess” of the input state cannot be assumed to be a good guess
unless we are sure that the steps leading to it have high accuracy. As we show
later, we have good evidence to believe that this is not so. Therefore, the
analysis comparing the twirled state estimate from the maximum likelihood
reconstruction with $\rho_{true}$, shown in the next section, should not be
treated as tests of the performance of the reconstruction process.

Our analysis was mostly restricted to the input state associated with the HWP
setting $\theta_{HWP}=76.74^{\circ}$. The corresponding
$\ket{\psi_{\textrm{in}}}$ estimate has magnitude
$\absolutevalue{\alpha}=0.8324$, polarization angle $\cos(\theta_{p})=0.7388$
(or $\theta_{p}=42.37^{\circ}$) and overlap angle $\cos(\theta_{t})=0.7138$.
For comparison, the polarization angle $\theta_{p}$ extracted from
$\theta_{HWP}$ is $43.88^{\circ}$, which differs from the estimate that we
trust more by almost $1.5^{\circ}$. Of course, there is significant
uncertainty in these numbers, which is mostly due to the error propagation
from the uncertainties in the BS coefficients, in the operators in Eq. 6.2.1
that are used to estimate the magnitude in each output, in the LO magnitude
and phase, which are all used to obtain $\ket{\psi_{\textrm{in}}}$. We
truncate $\ket{\psi_{\textrm{in}}}$ at $N=4$ maximal photon-number, followed
by normalization. The fidelity between the truncated state $\rho_{true}$ and
$\ket{\psi_{\textrm{in}}}$ is given by Eq. 5.3.1, which in our case equals
$0.9992$. We describe the “true state” estimates of the other $4$ input states
at the end of the next section in the context of the discussion of the
discrepancies that caused us to terminate the analysis.

###  6.3 State Reconstruction Analysis

In this section we present the results of reconstruction of the three-mode
twirled state $\tilde{\rho}_{3}$ using maximum likelihood estimation. We first
describe how we constructed IC measurement contexts, and then describe and
compare the results of the reconstruction with the relative frequencies
associated with the different measurement contexts. We spend significant time
investigating the extent and causes of model mismatch, and argue towards the
end why we decided to abandon the further analysis of the experiment.

####  6.3.1 Construction of Measurement Contexts

According to our strategy of binning the experimental outcomes according to
the corresponding phase of the LO, there are three LO magnitudes and $180$
phases for each magnitude to construct POVM elements from. For a WFH
configuration with a BS that is characterized by partition size $K=2$ and
$S_{1}>1$, according to Thm. 4.3.3 and Prop. 3, $(N+1)^{2}=25$ different LO
amplitudes are sufficient to construct an IC measurement context. We
numerically find (using the procedure described in Sec. 5.2.1) that even a
smaller number is sufficient, but when only $3$ LO magnitudes are available,
the number of phases associated with each LO magnitude must be at least $10$.
This is, of course, conditioned on our choice to use a maximum photon number
of four. Let us denote the three magnitudes by $R_{1}$, $R_{2}$ and $R_{3}$ by
the increasing order of their values. For the input state with
$\theta_{HWP}=76.74^{\circ}$ these were estimated to be $R_{1}=1.082$,
$R_{2}=1.286$ and $R_{3}=1.531$. To have an equal representation across the
magnitudes, we decided to use $10$ different phase values for each magnitude
to construct $18$ different measurement contexts that have mutually distinct
POVM elements. We perform reconstruction with each measurement context to
observe the sensitivity of the estimation procedure to different datasets.

The $18$ different IC measurement contexts were constructed as follows. Let us
refer to these measurement contexts by a number ranging from $1$ to $18$. The
measurement context $\\#1$ is constructed from the following LO amplitudes:
The phase values $\\{(1+36i)^{\circ}\mid i=0,1,\ldots,9\\}$ are used for the
magnitude $R_{1}$, while for the magnitudes $R_{2}$ and $R_{3}$ we use the
phase values $\\{(13+36i)^{\circ}\mid i=0,1,\ldots,9\\}$ and
$\\{(25+36i)^{\circ}\mid i=0,1,\ldots,9\\}$, respectively. So the phase values
for the different magnitudes are mutually different, and the set of phase
values for each magnitude are equidistantly distributed along the whole phase
range. The other $17$ measurement contexts are constructed similarly: the
phase values for $R_{1}$ are chosen to be $\\{(n+36i)^{\circ}\mid
i=0,1,\ldots,9\\}$, for $R_{2}$ they are $\\{(n+12+36i)^{\circ}\mid
i=0,1,\ldots,9\\}$, and for $R_{3}$ they are $\\{(n+24+36i)^{\circ}\mid
i=0,1,\ldots,9\\}$ where $n\in\\{3,5,\ldots,35\\}$ and corresponds to the the
measurement context $\\#(n+1)/2$. The reader can verify that these $18$
measurement contexts are built with mutually different LO amplitudes. Further,
this choice of measurement contexts insures that all samples that were
generated from the interference experiments are used in our analysis.

####  6.3.2 Reconstruction With Different Measurement Contexts

For the reconstruction we use the diluted $R\rho R$ algorithm described in the
last chapter. As a reminder, this algorithm starts iterating with the $R\rho
R$ algorithm until either a stopping condition is met or the increase in the
value of the logarithmic likelihood during the iteration is less than some
predetermined value $\Delta\mathcal{L}$. If the latter, the algorithm
transitions to a diluted $R\rho R$ regime, with the value of the initial
dilution parameter $\epsilon$ being $10^{30}$. If the stopping condition is
not met and the change in logarithmic likelihood during an iteration becomes
less than $\Delta\mathcal{L}$, the value of $\epsilon$ is halved. This
continues until either the stopping condition is met at some point of the
iterations, or the value of $\epsilon$ reaches $10^{-30}$ at which point we
terminate the iterations. The stopping condition is based on calculating the
maximum eigenvalue of $\hat{R}$ (Eq. 5.2.20) and terminating the algorithm if
it is less than $1+\texttt{r}$ where r is decided beforehand, as described in
the next paragraph. The number of samples for each measurement context
averages to around $26.5\cdot 10^{6}$ which are almost equally divided among
the $30$ different LO amplitudes. Our simulations with coherent states in the
previous chapter suggest that if the experiment perfectly matched the model,
we should expect such a sample size to result in a value of fidelity $\geq
0.999$ between the estimate and $\rho_{true}$ if the estimate is close enough
to the state that maximizes the likelihood. Note that $\rho_{true}$ is
truncated, but the state that generated the experimental data,
$\ket{\psi_{\textrm{in}}}$, was not truncated, which is different from the
simulations where we used the truncated state to simulate data. However, since
the fidelity between $\ket{\psi_{\textrm{in}}}$ and $\rho_{true}$ is $\geq
0.9992$, the fidelity between the estimate and $\rho_{true}$ should still be
near $0.999$ if there is no model mismatch.

In our analysis we decided to set $\Delta\mathcal{L}$ to zero, and used a
value of $\texttt{r}=200/M$, where $M$ is the number of samples, for the
stopping condition. Our simulations in the last chapter suggest that this
should result in as good of a performance as we can expect from our algorithm.
The value of r was chosen to terminate the reconstruction in a reasonable time
(around one hour). Some experimentation with smaller values of r showed that
the termination time of the algorithm would increase substantially with little
gain in the performance. In particular, we found that the value of
$\texttt{r}=1/M$ which we commonly employed with our simulations in the last
chapter did not result in a termination of the algorithm even after several
days. We did not observe such an issue with our simulations. One likely reason
is that all our simulations were with a WFH configuration with a BS
characterized by a trivial partition, which required us to encode our
variables and perform computations on matrices in $\mathcal{R}_{2,N}^{s}$. In
contrast, we encode our state and perform computation on matrices in
$\mathcal{R}_{3,N}^{s}$ which have an order $N$ times more parameters -
increasing the computation time and accumulated error in each iteration step.

The fidelity between each of the $18$ estimates and $\rho_{true}$ is shown in
Fig 6.5. It can be seen that the fidelities are consistently within $0.985\pm
0.002$. While these are not bad results, they fall short of the values
(consistently around $0.999$) we observed during our simulations with coherent
states in the last chapter. This is expected because, as we argued towards the
end of Subsec. 6.2.3, $\rho_{true}$ cannot be taken to be an accurate guess of
the true state. The fact that the fidelities observed for all $18$ estimates
are similar, even though the data they were reconstructed from are mutually
distinct, suggests that the estimates are close to each other. Indeed,
comparing the estimates with each other, we find that the fidelity between any
two estimates is above $0.9975$ and there are pairs which have a fidelity well
above $0.9995$ (Fig. 6.6). These facts suggest that we are not handicapped by
statistical noise or convergence issues in the reconstruction algorithm.
Rather, we are observing a model mismatch, such that the frequencies generated
from the data for all phase values and magnitudes of the LO are not
represented very well by the total photon statistics at the outputs of the FBS
resulting from the interference of $\rho_{true}$ with the LO by the FBS with
the coefficients we assume in our model. In the next subsection we aim to
investigate the degree of the model mismatch and where the model errs
significantly.

Figure 6.5: The fidelity between the estimate and $\rho_{true}$ for each of
the $18$ measurement contexts. The HWP angle of the input state is
$\theta_{HWP}=76.74^{\circ}$.  Figure 6.6: The fidelities between the $18$
estimates. The HWP angle of the input state is $\theta_{HWP}=76.74^{\circ}$.

####  6.3.3 Studies With Parametric Bootstrapping

One common method to inspect the validity of a statistical model, is to
perform parametric bootstrapping. Parametric bootstrapping uses the estimate
obtained from a particular measurement context and associated experimental
frequencies to calculate the probabilities given a model and employs a random
number generator to generate samples from these probabilities. The number of
generated samples is equal to the number of samples generated by the
experiment, and the bootstrapped estimate is obtained from these samples by
the same procedure that was used to obtain the original estimate from the
experimental samples. Ideally, this is done many times to obtain a number of
bootstrapped estimates and thus to get an idea about the distribution of the
estimator when the true state is the original estimate. Bootstrapping can be
used to generate confidence intervals for the estimate, but if there is
significant model mismatch, the confidence intervals may have coverage
probabilities that are different from their nominal confidence levels.

We use parametric bootstrapping to calculate and compare likelihood ratios
(LRs). In particular, for the original estimate and for each bootstrap we have
the final value of $\mathcal{L}$ as well as the logarithmic likelihood of the
unconstrained model (that is, where the probabilities are taken to be equal to
the frequencies). We compare the LRs of the bootstraps with the LR of the
estimate. If the model is accurate, the sample size is large enough, and the
estimate is close to the state that maximizes $\mathcal{L}$, we expect the
difference between the LRs of the bootstraps and the LR of the estimate to be
near zero. A histogram based on $60$ bootstraps generated for the estimate
obtained using the measurement context $\\#1$ is shown in Fig. 6.7, where the
logarithmic LR is defined as $-2(\mathcal{L}-\mathcal{L}_{u})$ with
$\mathcal{L}_{u}$ being the logarithmic likelihood of the unconstrained model.
The logarithmic LR of the estimate is $226$ times the sample standard
deviation of the bootstrap values away from the mean of the bootstraps.

Figure 6.7: The histogram of the logarithmic likelihood ratio (LR) of the
bootstraps with the logarithmic LR of the estimate shown as a dot for
illustration. $60$ bootstraps were generated. Further details are in the text.

We performed further bootstrapping studies to be more confident that we are
observing model mismatch. We simulated the experimental model and the
bootstrapping procedure described above, with the “true state” being the
estimate obtained using the measurement context $\\#1$ and the corresponding
frequencies. Namely, we used this “true state” to generate samples by a random
number generator according to the measurement context $\\#1$ and then ran the
reconstruction algorithm on these frequencies to generate an estimate. We then
used this estimate to generate $12$ bootstraps and calculated the logarithmic
LR for each of them. The logarithmic LR value for the estimate lies within the
range of the logarithmic LR values of the bootstraps - in sharp contrast with
the $226.43\sigma$ deviation observed in the analysis of the experimental
data.

We tried to investigate how much the various approximations and assumptions
used in the data generation procedure might contribute to the model mismatch.
In particular, one approximation we make is binning the outcomes of the TESs
into $2^{\circ}$ phase bins, and assigning the LO amplitude with the mean
value of the bin to the outcomes. To get an idea about the impact of this
binning procedure when the model matches the process of data generation, we
again started with the estimate obtained from the experimental data using the
measurement context $\\#1$. In contrast to the simulations described in the
previous paragraph, we used a measurement context consisting of $10$ times
more LO phases to generate the samples. In particular, we chose $10$ equally
spaced phase values falling within $1^{\circ}$ in each direction of each LO
phase appearing in the measurement context $\\#1$. We generated samples using
this measurement context followed by binning the samples into $2^{\circ}$
phase bins as in the generation of experimental data, and used the measurement
context $\\#1$ for reconstruction. We then produced $24$ bootstraps from the
estimate using the measurement context $\\#1$ for generation of frequencies
and looked at the logarithmic LR of the estimate and the bootstraps. We again
found that the logarithmic LR of the estimate falls in the range of those of
the bootstraps, indicating that the effect of phase binning on the estimator
is of the order of the statistical noise from an order of $10^{7}$ samples.

We report yet another investigation we did with phase binning. In the
reconstruction procedure from the experimental data, instead of using the
measurement context $\\#1$, we used slightly different POVM elements for each
LO amplitude. In particular, the POVM element corresponding to the observation
of $k$ and $l$ photons from a particular phase bin was constructed as an
equally weighed average of operators corresponding to the same photon number
observations but associated with different phases in the phase bin. To
clarify, the new POVM elements are the numerical integrals of the precise-
phase POVM elements in the corresponding phase intervals. We used $10$ equally
spaced phase values for each phase interval for the numerical integration. We
used this modified measurement context to generate the estimate from the
experimental data, followed by parametric bootstrapping with the same model
using $55$ bootstraps. We again observed $\approx 240\sigma$ deviation in the
logarithmic LR of the estimate from the mean of the bootstraps.

We conclude that the model mismatch we observe when comparing the LRs of the
bootstraps with the LR of the estimate obtained from the experimental data
cannot be due to assigning slightly different phase values to the samples and
using a measurement context associated with the corresponding LO phases. The
sensitivity of the statistical distribution of logarithmic LRs to such
deviations in the POVMs is negligible when the sample size is $\approx
2.7\cdot 10^{7}$.

####  6.3.4 Investigation of the Discrepancies

This is a good place to recount which of the parameters that we use in our
data analysis are given to us from measurements preceding the experiment (that
we take as being their true values), and which of the parameters that we use
in the reconstruction are estimated from the experimental data. The parameters
we assume are true are the reflection and transmission coefficients of the BS
for each polarization, the matrices $T_{1}$ and $T_{2}$, the transmission
coefficients in the couplings of the FBS outputs to the fibers leading to the
TESs which we incorporate into $T_{1}$ and $T_{2}$ producing $T^{\prime}_{1}$
and $T^{\prime}_{2}$ as in Eq. 6.2.1, and the bounds for the TESs (used for
assigning photon numbers to the outcomes) up to a scaling factor. The rest of
the parameters that are used in the data analysis - namely, the scaling
factors for the bounds, and the magnitudes and phases of the LO, are estimated
from the “LO state only” datasets and from a part of the datasets of the
interference experiments using the parameters above.

While it is not clear how/if the bounds for the TESs or the matrices
$T^{\prime}_{1}$ and $T^{\prime}_{2}$ can be double-checked using the
available data, the BS coefficients corresponding to the polarization of the
LO are straightforward to verify (assuming the bounds and the matrices
$T^{\prime}_{1}$ and $T^{\prime}_{2}$ are exact) from the “LO state only”
datasets using the Poissonian model for both TESs. There are $15$ such
datasets - one for each of the $3$ LO magnitudes for each of the $5$ different
input states. The estimates for the splitting ratio from these datasets are in
a tight range from $48.38:51.62$ to $48.61:51.39$, while the corresponding
splitting ratio used in the model is $49.46:50.54$ \- a $1.9\%$ difference
from the average of the estimates.

A stronger disagreement between the value of a quantity measured before the
experiment and our estimate of that quantity based on the experimental data
pertains to the temporal overlap $\theta_{t}$ in the input state. In
particular, before each interference experiment in each day, the temporal
overlap between the input state and the LO was maximized by measuring the
visibility. A temporal overlap of $\geq 95\%$, corresponding to
$\cos(\theta_{t})\geq 0.974$, for each input state was reported. In contrast,
we found by the procedure described in Subsec. 6.2.3 that the values of
$\cos(\theta_{t})$ were well below that number for all input states. These
numbers, along with the other parameters that describe the estimate of the
true state, are detailed in Table 6.1 for completeness. The value of the
temporal overlap reported is obtained according to the procedure described in
Subsec. 6.2.3. Another thing to note is that the values of the temporal
overlap according to the two counters were found to be somewhat different from
each other, and this difference became more pronounced for larger $\theta_{p}$
values. In particular, for the input state with the HWP setting
$\theta_{HWP}=90.66^{\circ}$, which has the highest value of $\theta_{p}$, the
temporal overlaps were $33.09\%$ and $38.63\%$ according to the two different
TESs.

Input state (HWP setting) | $\absolutevalue{\alpha}$ | $\theta_{p}$ |  $\theta_{p}$ according to $\theta_{HWP}$ | $\cos(\theta_{t})$ | Temporal overlap
---|---|---|---|---|---
$54.8^{\circ}$ | 0.8145 | $0^{\circ}$ |  $0^{\circ}$ | 0.8237 | $67.85\%$
$64.8^{\circ}$ | 0.8080 | $15.97^{\circ}$ |  $20^{\circ}$ | 0.7897 | $62.36\%$
$76.74^{\circ}$ | 0.8324 | $42.37^{\circ}$ |  $43.88^{\circ}$ | 0.7138 | $50.95\%$
$84.59^{\circ}$ | 0.8337 | $58.66^{\circ}$ |  $59.58^{\circ}$ | 0.6853 | $46.97\%$
$90.66^{\circ}$ | 0.8289 | $71.05^{\circ}$ |  $71.72^{\circ}$ | 0.5986 | $35.84\%$
Table 6.1: The estimates of the input states using the method described in
Subsec. 6.2.3. $\absolutevalue{\alpha}$, $\theta_{p}$ and $\theta_{t}$ form
the input state according to Eq. 6.2.5. The temporal overlap is defined as
$\cos^{2}(\theta_{t})\cdot 100\%$.

One useful test we performed to investigate the nature of the model mismatch
was to look at the sum of the magnitude-squares estimated at the outputs of
the FBS for different values of the relative phase during the interference
experiment. In particular, if our model of the experiment is correct, we
expect the interference fringes to combine to form a flat horizontal line.
Remember that we used the data of the TESs during one acquisition window to
estimate the relationship between the voltage across the piezo and the
relative phase between the LO and the input state. We then used this estimate
to assign a relative phase to the outcomes of the TESs during the consequent
acquisition window and binned the outcomes into $\approx 2^{\circ}$ phase bins
according to their corresponding phase value. We performed this for the whole
dataset of each interference experiment for the input state with
$\theta_{HWP}=76.74^{\circ}$. Thus, the magnitude-squares at the outputs of
the FBS for different values of the relative phase can be estimated using the
Poissonian model on the samples for each phase bin. We can do this for the
samples extracted from the whole dataset as well as for the samples extracted
from a single acquisition window. We show the scatter plots of the estimated
magnitude-squares in each phase bin for both TESs, as well as their sum at
each phase value, for a sample size coming from $8$ acquisition windows of the
interference experiment with the weakest LO magnitude in Fig. 6.8a. The next
subfigure (Fig. 6.8b) is based on the whole set of samples extracted from the
same interference experiment. It can be seen that the combination of the
magnitude=squares exhibits a clear sinusoidal pattern across the phase values.
The relative height of this sinusoid is $\approx 2\%$ of the mean of the sum
of the magnitude-squares. Similar results are observable for the other two LO
magnitudes.

Figure 6.8: (a) The estimates of the magnitude-squares at the outputs of the
FBS and their sum using samples from $8$ acquisition windows of an
interference experiment. The input state is $\theta_{HWP}=76.74^{\circ}$ and
the LO has magnitude $1.082$.

.

Figure 6.9: (b) The estimates of the magnitude-squares at the outputs of the
FBS and their sum using all the samples extracted from the data of an
interference experiment. The input state is $\theta_{HWP}=76.74^{\circ}$ and
the LO has magnitude $1.082$. (cont.)

One culprit could be the memory effect in the TESs. In particular, it was
established that the TES does not fully reset during the $10\textrm{
}\mu\textrm{s}$ time interval between the shots. This implies that the initial
conditions of the TES signal at the beginning of the absorption of the
following pulse are slightly higher than the fully relaxed conditions, which
causes the MFV estimate to be slightly higher than it would be otherwise. This
in turn can cause our estimates of the coherent state amplitudes using the
Poissonian model to be higher than they are. If the memory effect is stronger
for the larger magnitudes, it is plausible that this could cause an apparent
phase shift in the interference fringes. We worked on methods to filter out
this bias in the final MFV estimate, and the experimental group redid the MFV
extraction analysis from the signals of the TESs using this improved method
for a part of the experimental data for each input state. In fact, the
estimated parameters for the “true state” reported in Table 6.1 for all input
states besides the one with $\theta_{HWP}=74.74^{\circ}$ were extracted using
the new datasets. However, this did not appear to affect the dependence of the
sum of the magnitude-squares on relative phase, and a similar degree of
discrepancy was still observable.

The shift in the phase between the sinusoids corresponding to the two TESs
seems to be around $4^{\circ}$ shy of $180^{\circ}$. It is interesting to
explore how much the model mismatch would express itself in the value of the
LR if the phases assigned to POVMs are all shifted in the same direction from
the correct values by $4^{\circ}$. In particular, we perform a bootstrapping
simulation, starting with the estimate produced using the measurement context
$\\#1$ as our true state. We use the measurement context $\\#1$ to generate
random samples but use a modified version of the measurement context $\\#1$
with all LO phases shifted by $4^{\circ}$ in the same direction to perform
reconstruction. We generate $24$ bootstraps from the resulting estimate, but
using the correct LO phases. We find that the deviation of the logarithmic LR
of the estimate is $\approx 2.7\sigma$ from the mean of the bootstraps. So it
seems that the large model mismatch that we observe in the reconstruction from
the experimental data cannot be explained by assuming that for whatever reason
there is a systemic error in the assignment of the relative phase to the
outcomes of the TESs.

These discrepancies, and our inability to explain them and, ideally, to
incorporate them in our model, forced us to abandon the data analysis of the
experiment. Currently, there is an ongoing effort to set up and conduct a
similar experiment to demonstrate the theory developed in this thesis, and
care is taken to utilize the lessons from our failures to account for many
possible sources of uncertainty that affect the data analysis.

## References

  * [1] A. Allevi, M. Bina, S. Olivares, and M. Bondani, Homodyne-like detection scheme based on photon-number-resolving detectors, International Journal of Quantum Information, (2017), p. 1740016.
  * [2] H. A. Bachor and T. C. Ralph, A Guide to Experiments in Quantum Optics, Wiley-VCH, third ed., Oct 2019.
  * [3] K. M. Beck, M. Hosseini, Y. Duan, and V. Vuletić, Large conditional single-photon cross-phase modulation, Proceedings of the National Academy of Sciences, 113 (2016), pp. 9740–9744.
  * [4] M. Bina, A. Allevi, M. Bondani, and S. Olivares, Phase-reference monitoring in coherent-state discrimination assisted by a photon-number resolving detector, Scientific Reports, 6 (2016), p. 26025.
  * [5] M. Bina, A. Allevi, M. Bondani, and S. Olivares, Homodyne-like detection for coherent state-discrimination in the presence of phase noise, Opt. Express, 25 (2017), pp. 10685–10692.
  * [6] G. Brida, L. Ciavarella, I. P. Degiovanni, M. Genovese, L. Lolli, M. G. Mingolla, F. Piacentini, M. Rajteri, E. Taralli, and M. G. A. Paris, Quantum characterization of superconducting photon counters, New Journal of Physics, 14 (2012), p. 085001.
  * [7] P. Busch, M. Grabowski, and P. Lahti, Operational Quantum Physics, Lecture Notes in Physics Monographs, Springer Berlin, Heidelberg, 1995.
  * [8] G. Casella and R. L. Berger, Statistical Inference, Thomson Learning, second ed., 2002.
  * [9] G. M. D’Ariano, M. G. A. Paris, and M. F. Sacchi, Quantum tomography, 2003.
  * [10] J. Diestel and A. Spalsbury, The Joys of Haar Measure, Graduate studies in mathematics, American Mathematical Society, 2014.
  * [11] G. Donati, T. J. Bartley, X. min Jin, M. D. Vidrighin, A. Datta, M. Barbieri, and I. A. Walmsley, Observing optical coherence across Fock layers with weak-field homodyne detectors, Nature communications, 5 (2014), p. 5584.
  * [12] A. Einstein, Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen, Annalen der Physik, 322 (1905), pp. 549–560.
  * [13] W. Fulton and J. Harris, Representation Theory: A First Course, no. 129 in Graduate Texts in Mathematics, Springer, New York, 1991.
  * [14] T. Gerrits, A. Lita, B. Calkins, and S. W. Nam, Superconducting Transition Edge Sensors for Quantum Optics, Springer International Publishing, Cham, 2016, pp. 31–60.
  * [15] S. Glancy, E. Knill, and M. Girard, Gradient-based stopping rules for maximum-likelihood quantum-state tomography, New Journal of Physics, 14 (2012), p. 095017.
  * [16] R. J. Glauber, Coherent and incoherent states of the radiation field, Phys. Rev., 131 (1963), pp. 2766–2788.
  * [17] R. J. Glauber, The quantum theory of optical coherence, Phys. Rev., 130 (1963), pp. 2529–2539.
  * [18] R. Haag, Local Quantum Physics, Springer, 1992.
  * [19] F. Hernández and C. J. Riedel, Rapidly decaying Wigner functions are Schwartz functions, Journal of Mathematical Physics, 63 (2022), p. 022104.
  * [20] H. Hertz, Ueber einen Einfluss des ultravioletten Lichtes auf die electrische Entladung, Annalen der Physik, 267 (1887), pp. 983–1000.
  * [21] Z. Hradil, J. Rehacek, J. Fiurasek, and M. Jezek, Maximum-likelihood methods in quantum mechanics, in Quantum State Estimation, Springer-Verlag, New York, 2004, pp. 163–172.
  * [22] J. Řeháček, Z. Hradil, E. Knill, and A. I. Lvovsky, Diluted maximum-likelihood algorithm for quantum tomography, Phys. Rev. A, 75 (2007), p. 042108.
  * [23] M. Koashi, K. Kono, M. Matsuoka, and T. Hirano, Probing the two-photon phase coherence of parametrically down-converted photons by a local oscillator, Phys. Rev. A, 50 (1994), pp. R3605–R3608.
  * [24] A. Kuzmich, I. A. Walmsley, and L. Mandel, Violation of Bell’s inequality by a generalized Einstein-Podolsky-Rosen state using homodyne detection, Phys. Rev. Lett., 85 (2000), pp. 1349–1353.
  * [25] U. Leonhardt, Measuring the Quantum State of Light, Cambridge University Press, Cambridge, UK, 1997.
  * [26] U. Leonhardt, M. Munroe, T. Kiss, T. Richter, and M. Raymer, Sampling of photon statistics and density matrix using homodyne detection, Optics Communications, 127 (1996), pp. 144 – 160.
  * [27] R. A. Lorentz, Multivariate Birkhoff Interpolation, Lecture Notes in Mathematics, Springer Berlin, Heidelberg, 1992.
  * [28] A. I. Lvovsky and M. G. Raymer, Continuous-variable optical quantum-state tomography, Rev. Mod. Phys., 81 (2009), pp. 299–332.
  * [29] D. Maxwell, D. J. Szwer, D. Paredes-Barato, H. Busche, J. D. Pritchard, A. Gauguet, K. J. Weatherill, M. P. A. Jones, and C. S. Adams, Storage and control of optical photons using Rydberg polaritons, Phys. Rev. Lett., 110 (2013), p. 103001.
  * [30] A. Migdall, S. V. Polyakov, J. Fan, and J. C. Bienfang, Single-Photon Generation and Detection, vol. 45 of Experimental Methods in the Physical Sciences, Academic Press, 2013.
  * [31] C. R. Müller, M. A. Usuga, C. Wittmann, M. Takeoka, C. Marquardt, U. L. Andersen, and G. Leuchs, Quadrature phase shift keying coherent state discrimination via a hybrid receiver, New Journal of Physics, 14 (2012), p. 083009.
  * [32] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, Cambridge University Press, Cambridge, UK, 2001.
  * [33] S. Olivares, A. Allevi, G. Caiazzo, M. G. A. Paris, and M. Bondani, Quantum tomography of light states by photon-number-resolving detectors, New Journal of Physics, 21 (2019), p. 103045.
  * [34] T. Peyronel, O. Firstenberg, Q.-Y. Liang, S. Hofferberth, A. Gorshkov, T. Pohl, M. Lukin, and V. Vuletic, Quantum nonlinear optics with single photons enabled by strongly interacting atoms, Nature, 488 (2012), pp. 57–60.
  * [35] D. S. Phillips, Advanced measurements for quantum photonics and quantum technologies, PhD thesis, University of Oxford, 2020.
  * [36] G. Puentes, A. Datta, A. Feito, J. Eisert, M. B. Plenio, and I. A. Walmsley, Entanglement quantification from incomplete measurements: applications using photon-number-resolving weak homodyne detectors, New Journal of Physics, 12 (2010), p. 033042.
  * [37] G. Puentes, J. S. Lundeen, M. P. A. Branderhorst, H. B. Coldenstrodt-Ronge, B. J. Smith, and I. A. Walmsley, Bridging particle and wave sensitivity in a configurable detector of positive operator-valued measures, Phys. Rev. Lett., 102 (2009), p. 080404.
  * [38] K. J. Resch, J. S. Lundeen, and A. M. Steinberg, Quantum state preparation and conditional coherence, Phys. Rev. Lett., 88 (2002), p. 113601.
  * [39] K. Schmüdgen, The Moment Problem, Graduate Texts in Mathematics, Springer International Publishing, 2017.
  * [40] A. Serafini, Quantum Continuous Variables: A Primer of Theoretical Methods, CRC Press, 2017.
  * [41] S. Slussarenko and G. J. Pryde, Photonic quantum information processing: A concise review, Applied Physics Reviews, 6 (2019), p. 041303.
  * [42] D. Steck, Quantum and atom optics, 2007.
  * [43] A. Stuart, J. K. Ord, and S. Arnold, Advanced Theory of Statistics, Volume 2A: Classical Inference and the Linear Model, Oxford University Press, sixth ed., 1999.
  * [44] E. C. G. Sudarshan, Equivalence of semiclassical and quantum mechanical descriptions of statistical light beams, Phys. Rev. Lett., 10 (1963), pp. 277–279.
  * [45] G. S. Thekkadath, D. S. Phillips, J. F. F. Bulmer, W. R. Clements, A. Eckstein, B. A. Bell, J. Lugani, T. A. W. Wolterink, A. Lita, S. W. Nam, T. Gerrits, C. G. Wade, and I. A. Walmsley, Tuning between photon-number and quadrature measurements with weak-field homodyne detection, Phys. Rev. A, 101 (2020), p. 031801.
  * [46] J. D. Thompson, T. L. Nicholson, Q.-Y. Liang, S. H. Cantu, A. V. Venkatramani, S. Choi, I. A. Fedorov, D. Viscor, T. Pohl, M. D. Lukin, and V. Vuletic, Symmetry-protected collisions between strongly interacting photons, Nature, 542 (2017), pp. 206–209.
  * [47] J. Tiedau, V. S. Shchesnovich, D. Mogilevtsev, V. Ansari, G. Harder, T. J. Bartley, N. Korolkova, and C. Silberhorn, Quantum state and mode profile tomography by the overlap, New Journal of Physics, 20 (2018), p. 033003.
  * [48] S. Wallentowitz and W. Vogel, Unbalanced homodyning for quantum state measurements, Phys. Rev. A, 53 (1996), pp. 4528–4533.
  * [49] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Rev. Mod. Phys., 84 (2012), pp. 621–669.
  * [50] H. S. Wilf, Generatingfunctionology, A K Peters/CRC Press, third ed., 2005.
  * [51] H. Xu, F. Xu, T. Theurer, D. Egloff, Z.-W. Liu, N. Yu, M. B. Plenio, and L. Zhang, Experimental quantification of coherence of a tunable quantum detector, Phys. Rev. Lett., 125 (2020), p. 060404.
  * [52] L. Zhang, H. B. Coldenstrodt-Ronge, A. Datta, G. Puentes, J. Lundeen, X.-M. Jin, B. Smith, M. Plenio, and I. Walmsley, Mapping coherence in measurement via full quantum tomography of a hybrid optical detector, Nature Photonics, 6 (2012), p. 364.
  * [53] L. Zhang, A. Datta, H. B. Coldenstrodt-Ronge, X.-M. Jin, J. Eisert, M. B. Plenio, and I. A. Walmsley, Recursive quantum detector tomography, New Journal of Physics, 14 (2012), p. 115005.

##  Appendix A A Result About Anti-Normally Ordered Operators

Here we derive a result about the relationship of the anti-normally ordered
powers of the total number operator and its usual powers. We could not find
this result in the literature, and since we derive a similar result for the
normally-ordered case, it seems like a good opportunity to include this result
in this thesis as well.

We use the notation introduced in Sec. 2.5. Let us introduce the anti-normal
ordering transformation on expressions of mode variables as the map that
rearranges the terms in each monomial by moving the variables representing
creation operators to the right of the monomial. We use triple vertical dots
to represent anti-normal ordering - namely, if $P$ is a polynomial of mode
variables then $\mathrel{\vdots}\mathrel{P}\mathrel{\vdots}$ is the
corresponding anti-normally ordered operator. A calculation similar to
LABEL:eq:monexprcov shows that anti-normal ordering commutes with PLTs. In
particular:

######  Lemma A.0.1.

$\displaystyle\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{\mathrel{\vdots}\mathrel{\boldsymbol{n}^{k}}\mathrel{\vdots}}]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\mathrel{\vdots}\mathrel{\boldsymbol{n}^{k}}\mathrel{\vdots}}$
$\displaystyle=(\hat{n}+(S-1)+k)_{k},$ (A.0.1)

###### Proof.

The map
$P(\boldsymbol{a},\boldsymbol{a}^{\dagger})\mapsto\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{\mathrel{\vdots}\mathrel{P(\boldsymbol{a},\boldsymbol{a}^{\dagger})}\mathrel{\vdots}}]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\mathrel{\vdots}\mathrel{P(\boldsymbol{a},\boldsymbol{a}^{\dagger})}\mathrel{\vdots}}{\tmpbox}$
is an intertwiner of the group of PLTs, and $\boldsymbol{n}^{k}$ is invariant
under the PLT action on polynomial expressions. It follows that
$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{\mathrel{\vdots}\mathrel{\boldsymbol{n}^{k}}\mathrel{\vdots}}]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\mathrel{\vdots}\mathrel{\boldsymbol{n}^{k}}\mathrel{\vdots}}{\tmpbox}$
is also invariant under the action of PLTs, which in turn means (according to
Cor. 2.4.2) that it is a polynomial of the total number operator. We use
induction to obtain our proof. One can verify the formula in the lemma for
$S=1$ by direct calculation:

$\displaystyle a_{1}^{k}(a_{1}^{\dagger})^{k}$
$\displaystyle=\sum_{m=0}^{\infty}(m+1)\cdots(m+k)\ket{m}\bra{m}$
$\displaystyle=\left[\sum_{m_{1}=0}^{\infty}(m_{1}+1)\ket{m_{1}}\bra{m_{1}}\right]\left[\sum_{m_{2}=0}^{\infty}(m_{2}+2)\ket{m_{2}}\bra{m_{2}}\right]\cdots\left[\sum_{m_{k}=0}^{\infty}(m_{k}+k)\ket{m_{k}}\bra{m_{k}}\right]$
$\displaystyle=(a_{1}^{\dagger}a_{1}+1)(a_{1}^{\dagger}a_{1}+2)\cdots(a_{1}^{\dagger}a_{1}+k)=(a_{1}^{\dagger}a_{1}+k)_{k}.$
(A.0.2)

For general $S>1$, let $\boldsymbol{n^{\prime}}$ be the expression for the
total number operator on the first $S-1$ modes and $\boldsymbol{n}_{S}$ the
expression for the mode $S$ number operator. Then, assuming the formula holds
for the first $S-1$ modes, one obtains

(A.0.3)

We want to show that the last line in the above is equal to
$(\hat{n}+S-1+k)_{k}$. Since
$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{\mathrel{\vdots}\mathrel{\boldsymbol{n}^{k}}\mathrel{\vdots}}]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\mathrel{\vdots}\mathrel{\boldsymbol{n}^{k}}\mathrel{\vdots}}{\tmpbox}$
is a polynomial of $\hat{n}$, its expectations with all Fock state with total
photon number $m$ are equal. Thus, it suffices to calculate its expectation on
the state $\ket{m,0}$ with $m$ photons in the first $S-1$ modes and $0$ in the
last:

$\displaystyle\bra{m,0}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{\mathrel{\vdots}\mathrel{\boldsymbol{n}^{k}}\mathrel{\vdots}}]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\mathrel{\vdots}\mathrel{\boldsymbol{n}^{k}}\mathrel{\vdots}}{\tmpbox}$
$\displaystyle=\sum_{j=0}^{k}\binom{k}{j}(m+S-2+j)_{j}(k-j)_{k-j}$
$\displaystyle=k!\sum_{j=0}^{k}\binom{m+S-2+j}{j}$
$\displaystyle=k!\binom{m+S-1+k}{k}=(m+S-1+k)_{k}.$ (A.0.4)

∎