
\section{Optimization Algorithm}

Figure \ref{flowchart} shows a flowchart of the phase estimation algorithm implemented in real time in the FPGA to obtain the optimal values of the displacement field amplitude $|\beta_{opt}|$ and phase $\theta_{opt}$ for each adaptive step $k$, and to perform the Bayesian updating. This algorithm has two cases depending on the current adaptive step $k$. For $k\leq8$ the algorithm requires full calculation of the optimal values $|\beta_{opt}|$ and $\theta_{opt}$ using the complete prior phase distributions. For $k\geq9$ the algorithm approximates the prior phase distributions as Gaussian distributions allowing for a very efficient estimation of $|\beta_{opt}|$ and $\theta_{opt}$ in real time in the FPGA.

Starting at I in figure \ref{flowchart}, given the prior probability distribution $\mathrm{P}_{prior}(\phi)$ at adaptive step $k$, the FPGA calculates the optimal displacement amplitude $|\beta_{opt}|$ and phase $\theta_{opt}$ for the current measurement step (shown as I-III in Fig. \ref{flowchart}). For convenience, $\theta_{opt}$ is expressed as the sum of the current phase estimate $\hat{\phi}$ and an optimal phase offset $\Delta_{opt}$, $\theta_{opt} = \hat{\phi} + \Delta_{opt}$. Two cases can be identified depending on the measurement progression in $k$:

\textbf{Case I.--} For $k\leq8$ the optimal values $|\beta_{opt}|$ and $\Delta_{opt}$ are obtained from a lookup table (LUT) whose input is the complete detection record $\{d_{1}, d_{2}, ... d_{k-1}\}$ up to the adaptive step $k$ (see dashed blue rectangle). This LUT is pre-calculated for all possible detection records based on the maximization of the expected sharpness or the mutual information to provide the exact optimal values of $|\beta_{opt}|$ and $\Delta_{opt}$ using the full probability distribution $\mathrm{P}_{prior}(\phi)$ for every adaptive step up to $k$.

\textbf{Case II.--}  If $k\geq9$ the FPGA switches to a mode where the prior distribution  $\mathrm{P}_{prior}(\phi)$ for step $k$ is approximated as Gaussian (Gaussian approximation) in order to calculate approximately optimal values $|\beta_{opt}|$ and $\Delta_{opt}$ in an efficient way. This approximation is justified by observing that after $k \geq 9$ steps, the prior distribution $\mathrm{P}_{prior}(\phi)$ is close to Gaussian, so that it can be described by its first three moments: the mean, the variance $\sigma^{2}$, and skewness $\mu_{3}$ of the distribution, while its support over higher order moments is negligible. In particular, the skewness $\mu_{3}$ mainly describes the deviation of the prior $\mathrm{P}_{prior}(\phi)$ from a Gaussian distribution and provides information to correct the inferred LO phase offset $\delta$ under this approximation to obtain $\Delta_{opt}$ in step $k$.

To obtain the optimal values $|\beta_{opt}|$ and $\Delta_{opt}$ for $k\geq9$, the FPGA calculates the variance $\sigma^{2}$ and skewness $\mu_{3}$ of $\mathrm{P}_{prior}(\phi)$ (shown as II in Fig. \ref{flowchart}). The variance $\sigma^{2}$ is used to calculate the optimal amplitude $|\beta_{opt}|$ and phase $\delta$ under the Gaussian approximation using a second LUT with pre-calculated values (see dashed orange rectangle). Then, the phase $\delta$ obtained under this approximation is corrected by a function $f(\mu_{3})$ of skewness:

\begin{equation}
\Delta_{opt} = f(\mu_{3}) \times \delta \text{ , with  } f(\mu_{3})=-\frac{\textrm{sign}(\mu_{3})}{1+|\mu_{3}|}.
\label{skew}
\end{equation}

This function provides the correction from this approximation, and is empirically found through Monte Carlo simulations by comparing the optimal phases obtained from a Gaussian distribution and the full distribution with nonzero $\mu_{3}$. Note that calculating the optimal displacement values takes approximately 500 ns.

Given obtained $|\beta_{opt}|$ and $\Delta_{opt}$, the offset $\Delta_{opt}$ is added to the current phase estimate $\hat{\phi}$ at step $k$ to define the amplitude and phase $\theta_{opt}$ of the displacement operation for step $k$ (see dashed black rectangle). This displacement is implemented and a photon detection result $d_{k}$ for step $k$ is acquired (shown as V in Fig. \ref{flowchart}).

To obtain $\mathrm{P}_{post}(\phi)$ for step $k$, the FPGA obtains the likelihood function $\mathcal{L}(d_{k}|\phi, \beta)$ for the detected number of photons. Based on the photon detection result $d_{k}$ and the amplitude $|\beta|$ of local oscillator field, the FPGA retrieves from memory a segment of the likelihood function corresponding to the phase range of $\phi\in[0,\pi)$ (shown as VI in Fig. \ref{flowchart}). The FPGA then constructs the complete likelihood function using the fact that the likelihoods are symmetric about $\phi=\theta_{LO}$ and  periodic with period $2\pi$ (shown as VII in Fig. \ref{flowchart}). The posterior distribution $\mathrm{P}_{post}(\phi)$ for step $k$ is then calculated through Bayes' rule as the product of the likelihood $\mathcal{L}(d_{k}|\phi, \beta)$ and the prior probability distribution $\mathrm{P}_{prior}(\phi)$ (see dashed green rectangle). We note that the Bayesian updating uses the full prior probability distribution and the Gaussian prior approximation is used only for the LO optimization when $k\geq9$.

An updated phase estimate $\hat{\phi}$ is then obtained as the phase with maximum posterior probability. Next, the FPGA performs a quasi-normalization of the posterior (shown as VIII in Fig. \ref{flowchart}) by multiplying the posterior by $2^{p}$ where $p$ is an integer such that the most significant bit of the updated phase estimate (maximum of the posterior) is one. This prevents computational underflow in the Bayesian step of the algorithm so that the FPGA can maintain an accurate distribution for the phase for the entire measurement. The FPGA uses discretized probability distributions with 256 phase points and each probability value is represented in fixed point notation as a 10 bit number. Finally, given the posterior distribution $\mathrm{P}_{post}(\phi)$ for step $k$, the recursive Bayesian algorithm now defines the posterior $\mathrm{P}_{post}(\phi)$ in step $k$ as the prior $\mathrm{P}_{prior}(\phi)$ for the next step $k+1$ (shown as IX in Fig. \ref{flowchart}). The Bayesian updating takes approximately 500 ns, giving a feedback bandwidth of $\approx 1 \mathrm{MHz}$.

We note that as the input mean photon number $|\alpha|^{2}$ increases, the posterior phase distributions in general become very narrow (i.e. have small variance). With increasing $|\alpha|^{2}$, accurately representing the phase distributions in the simulations and in the experimental implementation requires increasing the number of discretization points for the phase. A larger number of discretization points will increase the computational complexity of the algorithm roughly linearly with the number of points. Within the FPGA, the increase in complexity is mainly due to an increase in the number of multiplications needed to calculate the posterior phase distribution, which increases the computation time and the memory required to store the likelihood functions. While the resources required for the experimental implementation increase with $|\alpha|^{2}$, in theory this does not affect the complexity of the proposed strategy for phase estimation. On the other hand, if the PNR used in the strategy is increased, the time required to compute the expected sharpness in Eq. (2), and mutual information in Eq. (3), will increase linearly with the number of terms in the sums in Eq. (2) and Eq. (3).

\begin{figure*}[!tp]
	\includegraphics[width = 0.9\textwidth]{VarianceLUTplots}
	\caption{Optimal LO amplitude and phase based on the Gaussian approximation for Gaussian prior distributions with different mean photon numbers $|\alpha|^2$ and variances $\sigma_{pr}^{2}$. (a,b) Optimal phase offset $\Delta_{opt}$ for estimation strategies maximizing (a) the expected sharpness and (b) the mutual information in each adaptive step. (c,d) Ratio of the optimal LO amplitude $|\beta_{opt}|$ to $|\alpha|$ on a log scale for strategies maximizing (c) the sharpness and (d) the mutual information. (e,f)  Ratio of the optimal LO amplitude $|\beta_{opt}|$ to the amplitude maximizing the Fisher information $|\beta_{FI}(\Delta_{opt})|$ from Eq. (\ref{optcond_delta}) for strategies maximizing (e) the sharpness and (f) the mutual information. Note the similarity for optimal phase offsets and amplitudes for both estimation strategies. Also note that the ratio $|\beta_{opt}|/|\beta_{FI}(\Delta_{opt})|$ approaches 1 as $\sigma_{pr}^{2}$ decreases for all mean photon numbers.
}
	\label{varlut}
\end{figure*}
\section{Gaussian Approximation}
In general, accurately describing the prior phase distribution $\mathrm{P}_{prior}(\phi)$ at step $k$ requires specifying many moments of the distribution. However, as mentioned above, in the adaptive estimation strategy based on recursive Bayesian updating, for $k \geq 9$ the prior $\mathrm{P}_{prior}(\phi)$ at step $k$ can be approximately described by its first three moments: mean, variance, and skewness. We find that the prior distribution can be well approximated as a Gaussian distribution with the mean equal to the current phase estimate $\hat{\phi}$ at the adaptive step $k$ and a variance $\sigma_{pr}^{2}$ equal to the variance of the prior distribution:
\begin{equation}
\sigma_{pr}^{2} = \int_{0}^{2\pi} (\phi - \hat{\phi})^{2} \mathrm{P}_{prior}(\phi),
\label{VarPrior}
\end{equation}
so that
\begin{equation}
\mathrm{P}_{prior}(\phi)\approx\tilde{\mathrm{P}}(\phi) = \mathcal{N}e^{-\frac{(\phi - \hat{\phi})^{2}}{2\sigma_{pr}^{2}}}
\label{Gaussian}
\end{equation}
where $\mathcal{N}$ is a normalization factor.

Given the calculated variance $\sigma_{pr}^{2}$ from Eq. (\ref{VarPrior}), it is possible to very efficiently obtain optimal values of the LO field for the amplitude $|\beta_{opt}|$ and the phase $\delta$ under the Gaussian approximation by optimizing the mutual information or sharpness for a Gaussian distribution. This is possible because the Gaussian distribution is uniquely defined by its mean and its variance, see Eq. (\ref{Gaussian}). Moreover, since we are looking for an optimal phase offset, we can find the optimal values for a zero-mean Gaussian, and then shift the phase offset $\Delta_{opt}$ by $\hat{\phi}$ to obtain the optimal phase $\theta_{opt}$ for the LO field. As a result, for a given  mean photon number $|\alpha|^{2}$, it is possible to construct a look up table (LUT) in the FPGA that associates the variance $\sigma_{pr}^{2}$ to the optimal values of the LO amplitude $|\beta_{opt}|$ and phase offset $\Delta_{opt}$, after the correction with the function $f(\mu_{3})$ in Eq. (\ref{skew}) (also see Fig. \ref{flowchart}).

Figure \ref{varlut} shows the calculated optimal values for the LO amplitude $|\beta_{opt}|$ and phase offset $\Delta_{opt}$ for different mean photon numbers $|\alpha|^{2}$ as a function of the variance of a Gaussian prior distribution $\tilde{\mathrm{P}}(\phi)$. Panels (a) and (b) show $\Delta_{opt}$ for estimation strategies maximizing the sharpness and the mutual information, respectively. Panels (c) and (d) show $|\beta_{opt}|$ for the same cases, respectively. We note that for small values of $\sigma_{pr}^{2}$, the optimal amplitude and phase values for the two different objective functions asymptote to the same value, which indicates that under the Gaussian approximation both estimation strategies are expected to show similar performances for large values of $k$.
\\
\\
\emph{\textbf{Asymptotic behavior}.---}
To further investigate the expected asymptotic behavior of the estimation strategies, we study the classical Fisher information (CFI) for measurements based on displacement and photon counting, which can provide insight into their performance compared to the fundamental lower bound, the Cramer-Rao lower bound. The CFI for displaced photon counting, with infinite photon number resolution (PNR), is given by:
\begin{align}
F(\phi) &= \sum_{n=0}^{\infty}\mathcal{L}(n|\phi) \Bigg( \frac{\partial}{\partial \phi}\mathrm{ln}\big(\mathcal{L}(n|\phi) \big) \Bigg) ^{2}
\nonumber
\\
\nonumber
\\
&= \frac{4|\alpha|^{2}|\beta|^{2}\mathrm{sin}^{2}(\phi)}{|\alpha|^{2} + |\beta|^{2} - 2|\alpha| |\beta| \mathrm{\cos}(\phi)}
\label{fieq}
\end{align}

Where $\mathcal{L}(n | \phi)$ is the likelihood function for detecting $n$ photons given the relative phase $\phi = \phi_{0} - \theta_{LO}$ between the input state and the LO. Note that when $|\beta|=|\alpha|$ the CFI of photon counting measurements is equal to the quantum Fisher Information (QFI) when $\phi=0$, which corresponds to performing displacement to the vacuum state. However, for estimation of an unknown phase, it is possible to find the optimal LO amplitude $|\beta|$ that maximizes $F(\phi)$:
\begin{equation}
|\beta_{FI}(\phi)|^{2} = \frac{|\alpha|^{2}}{\mathrm{cos}^{2}(\phi)}
\label{optcond}
\end{equation}

The evaluation of $F(\phi)$ in $|\beta_{FI}(\phi)|$ results in:
\begin{equation}
FI(\phi) \bigg \rvert_{|\beta|=|\beta_{FI}(\phi)|} = 4|\alpha|^{2} = QFI
\end{equation}

The CFI for displaced photon counting is guaranteed to be equal to the QFI for coherent states as long as Eq. (\ref{optcond}) is satisfied. Therefore, for any given relative phase $\phi$, there is a $|\beta_{FI}(\phi)|$ that saturates the QFI.

For the optimized estimation strategies based on adaptive measurements, after many adaptive steps the prior probability distributions approach Gaussian distributions with small variance $\sigma_{pr}^{2}$, and the estimate $\hat{\phi}$ approaches the true phase $\phi_{0}$. In this asymptotic limit, one may expect that the CFI of the estimation strategies approaches the QFI, and the amplitude $|\beta_{opt}|$ and phase $\Delta_{opt}$ approach values that satisfy the optimality condition in Eq. (\ref{optcond}), where $\Delta_{opt}$ describes the relative phase between the displacement field $\theta_{LO}$ and the actual phase $\phi_{0}$. That is:

\begin{equation}
\lim\limits_{\sigma_{pr} \to 0}|\beta_{opt}|^{2} = |\beta_{FI}(\Delta_{opt})|^{2} = \frac{|\alpha|^{2}}{\mathrm{cos}^{2}(\Delta_{opt})}
\label{optcond_delta}
\end{equation}

Figures \ref{varlut}(e) and \ref{varlut}(f) show the ratio $|\beta_{opt}|/|\beta_{FI}(\Delta_{opt})|$ of the optimal LO amplitude to the expected optimal amplitude from Eq. (\ref{optcond_delta}) given $\Delta_{opt}$ for strategies maximizing (e) the sharpness $\langle S(\beta, m) \rangle$ and (f) the mutual information $I(\beta, m)$. We observe that as $\sigma_{pr}$ decreases, the optimal amplitudes $|\beta_{opt}|$ obtained in the Gaussian approximation asymptote to $|\beta_{FI}(\Delta_{opt})|$ for $\mathrm{log_{10}}(\sigma_{pr}^{2}) < -1.5$ for all mean photon numbers. We note that while the optimal values $|\beta_{opt}|$ and $\Delta_{opt}$ in the estimation strategies without the Gaussian approximation during small number of steps $k$ do not satisfy the condition in Eq. (\ref{optcond_delta}), the value of $|\beta_{opt}|$ approaches $|\beta_{FI}(\Delta_{opt})|$ for large $k$.

\section{Preparation of initial phase}

In an ideal experiment for estimation of a completely unknown phase $\phi_{0}$, the initial phase of the  local oscillator (LO) displacement field $\theta_{LO}$ would be fixed, and the phase to be estimated would be chosen randomly $\phi_{0}\in[0,2\pi)$. This means that the initial relative phase $\theta_{LO}-\phi_{0}$ is random. However, since in our experimental setup we had arbitrary control of the phase $\theta_{LO}$ with high bandwidth for the adaptive measurements, we chose to fix $\phi_{0}$, and have the initial LO phase $\theta_{LO}$ be uniformly distributed, so that the initial relative phase is still random. In particular, in our experiment we set $\phi_{0} = \pi$ and randomly choose the initial LO phase $\theta_{LO}\in[0,2\pi)$, which is prepared with a 8-bit resolution DAC resulting in $2^8=256$ different possible initial phases.

To show that there is not any bias introduced by this procedure, we studied the phase estimator error $\hat{\phi} - \phi_{0}$ obtained from the estimation of $\phi_{0}$ with different initial LO phases $\theta_{LO}\in[0,2\pi)$. Figure \ref{error} shows the final phase error after $L=30$ adaptive steps for all experiments (total of $5\times 10^{4}$) for an estimation strategy maximizing the mutual information in each adaptive step. Each point represents one single estimation experiment, and the figure shows results for total mean photon numbers of $|\alpha|^{2}$ = 5 (orange), 10 (blue), 50 (red), 100 (green), and 500 (purple). We note the absence of any dependence of the final estimates $\hat{\phi}$ on the initial LO phase $\theta_{LO}$. Figure \ref{error} also shows an example of the distribution of final phase estimates $\hat{\phi}$ for $|\alpha|^{2} = 10$ showing normally distributed behavior. The variance Var[$\hat{\phi}$] of this distribution is approximately equal to the mean square error of the estimator for $|\alpha|^{2} = 10$.

\begin{figure}[!bp]
	\includegraphics[width = 8.5cm]{ErrorInitial_withhist_2}
	\caption{Estimator bias. Estimator error $\hat{\phi}-\phi_{0}$ for each initial LO phase $\theta_{LO}\in[0,2\pi)$ for every experimental run ($5\times 10^{4}$ samples) for different total mean photon numbers for a strategy maximizing the mutual information in each adaptive step. Note the independence of the error on the initial LO phase.}
	\label{error}
\end{figure}

\end{document}
