Fault-tolerant preparation of approximate GKP states

The phase estimation circuit for preparing an approximate $| \tilde{0}\rangle$ GKP state is given in figure 2. The H gate is the Hadamard gate and ${\rm{\Lambda }}({{\rm{e}}}^{{\rm{i}}\gamma })=\mathrm{diag}(1,{{\rm{e}}}^{{\rm{i}}\gamma })$. After applying several rounds of the circuit in figure 2, the input squeezed vacuum state (given in equation (2)) is projected onto an approximate eigenstate of S_p with some random eigenvalue⁷ ${{\rm{e}}}^{{\rm{i}}\theta }$. Additionally, an estimated value for the phase θ is obtained. After computing the phase, the state can be shifted back to an approximate +1 eigenstate of S_p.

Zoom In Zoom Out Reset image size

There are a variety of ways to choose the phases γ at each round in the gate ${\rm{\Lambda }}({{\rm{e}}}^{{\rm{i}}\gamma })$. In a non-adaptive phase estimation protocol, half of the values for γ can be chosen to be 0 and the other half can be chosen to be $\pi /2$. In an adaptive phase estimation protocol, if the phase to be estimated is θ and assuming no prior knowledge of θ, during the first round the phase can be chosen as ${\gamma }_{1}=0$. For later rounds, it was shown in [3] that the next phases can be chosen as

$$\begin{eqnarray}&&{\gamma }_{m}=\mathop{\arg \,\max }\limits_{\gamma }\displaystyle \sum _{{x}_{m}=0,1}| \int {\rm{d}}\theta {{\rm{e}}}^{{\rm{i}}\theta }{P}_{\gamma }(x[m]| \theta )| .\end{eqnarray} \tag{ 10 }$$

In equation (10), the probability ${P}_{\gamma }(x[m]| \theta )$ is the probability of obtaining measurement outcomes ${x}_{1},\,\cdots ,\,{x}_{m}$ when the produced output eigenstate $| {\psi }_{\theta }\rangle$ has eigenvalue ${{\rm{e}}}^{{\rm{i}}\theta }$. Since the measurement results for each round are independent, the estimated probability is given by

$$\begin{eqnarray}&&{P}_{\gamma ={\gamma }_{m}}(x[m]| \theta )=\displaystyle \prod _{i=1}^{m}{\cos }^{2}\left(\displaystyle \frac{\theta +{\gamma }_{i}}{2}+{x}_{i}\displaystyle \frac{\pi }{2}\right).\end{eqnarray} \tag{ 11 }$$

By analytically computing the evolution of the state, the exact probability of obtaining the outcomes $x[M]$ is derived in appendix D.

Given the final measurement record $x[M]$, the estimated phase $\tilde{\theta }$ is chosen as

$$\begin{eqnarray}&&\tilde{\theta }=\arg {\int }_{-\pi }^{\pi }{\rm{d}}\theta {{\rm{e}}}^{{\rm{i}}\theta }P(x[M]| \theta ).\end{eqnarray} \tag{ 12 }$$

The shift back correction based on the estimated phase $\tilde{\theta }$ is given by ${{\rm{e}}}^{{\rm{i}}\tfrac{\tilde{\theta }}{2\sqrt{\pi }}q}$. Note that for either the adaptive or non-adaptive protocol, the estimated phase and probabilities are computed using equations (11) and (12).

From the circuit of figure 2 and from equation (12), a measurement readout error can result in the wrong estimated phase which in turn results in an incorrect shift-back correction (application of ${{\rm{e}}}^{{\rm{i}}{v}{q}}$ where v is computed using equations (10) and (11)). Now suppose that the output state is afflicted by a shift error of the form ${{\rm{e}}}^{-{\rm{i}}{w}{q}}$. In this case, the final output state will have a total shift error of the form ${{\rm{e}}}^{-{\rm{i}}(v-w)q}$. If $| v-w| \mathrm{mod}\sqrt{\pi }\gt \sqrt{\pi }/6$, then the shift error will be uncorrectable and thus the phase estimation protocol will not be fault-tolerant⁸ . In what follows we will identify an output state of the phase estimation protocol as $x[m]$ if it arises from the measurement outcomes ${x}_{1},\,\cdots ,\,{x}_{m}$ in the fault-free case.

Consider the case where the phase estimation protocol is implemented in m rounds and let us assume that a single measurement readout error occurs and that all other operations are fault free. In this case the output state ${x}_{\mathrm{out}}[m]$ will have one bit which differs from the bit string ${x}_{\mathrm{corr}}[m]$ that would have been obtained in the fault free case (so that the Hamming distance ${d}_{{\rm{H}}}({x}_{\mathrm{out}}[m],{x}_{\mathrm{corr}}[m])=1$). The shift-back operation of equation (10) applied to the output state will be the shift correction of an output state that is one Hamming distance away from the actual state that was produced. Thus to ensure that the protocol is $(1,\tilde{\delta })$ fault-tolerant for some $\tilde{\delta }$, we should only accept output states ${x}_{\mathrm{out}}[m]$ which have the property that applying the shift back correction corresponding to any other state $x^{\prime} [m]$ with ${{d}}_{{\rm{H}}}(x^{\prime} [m],x[m])=1$ results in a correctable left-over shift error. These remarks motivate the following protocol to prepare an approximate $| \tilde{0}\rangle$ state which is fault-tolerant to a single measurement readout error.

• Calculation of acceptance set A_m and $\tilde{\delta }$. Consider all output states obtained from an m round fault-free phase estimation protocol of section 4.1. Let ${A}_{m}=\varnothing$ (which we call the acceptance set) and ${{\rm{\Gamma }}}_{m}=\varnothing$. For each output state $x[m]$, do the following

• $(1,\delta )$ fault-tolerant m round phase estimation protocol for measurement readout errors. Consider the acceptance set ${A}_{m}\ne \varnothing$ and $\tilde{\delta }$ computed from the procedure described above. During an implementation of the phase estimation protocol subject to measurement readout errors, if the obtained output state $x[m]\notin {A}_{m}$, abort the protocol and start anew.

Note that during the protocol, there can be more than one measurement readout error. However if more than one readout error occurs, it is possible that the output state is afflicted by a shift error of magnitude greater than $\sqrt{\pi }/6$. Our protocol only guarantees protection against a single readout error.

As an example, consider the four round phase estimation protocol. Applying the procedure described above for the adaptive phase estimation protocol, we find that A₄ is empty when choosing the initial phase to be ${\gamma }_{0}=0$ or ${\gamma }_{0}=\tfrac{\pi }{2}$. This indicates that the adaptive phase estimation protocol of four rounds described in section 4.1 is not fault-tolerant to single measurement readout errors⁹ . If the phases are computed using the non-adaptive phase estimation protocol and applying the protocol described above, we find that $\tilde{\delta }=\tfrac{\sqrt{\pi }}{6}-0.235\approx 0.0604$. The set A₄ of accepted measurement strings is given in table 1 and has $| {A}_{4}| =5$. The total probability of obtaining a state in A₄ is roughly $48.3 \%$. For the phase estimation protocol with more than four rounds, the acceptance set and the acceptance probability is very sensitive to the initial phase ${\gamma }_{0}$ and the domain of the optimized γ. For example, if we choose the initial phase ${\gamma }_{0}$ to be $\tfrac{\pi }{2}$ and the range of γ to be $[0,2\pi ]$, we get 18 accepted states (for the adaptive protocol). In this case the total probability of obtaining a state in A₈ is $6.25 \%$. If we choose the initial phase ${\gamma }_{0}$ to be 0, the acceptance set has four states and the probability of obtaining a state in A₈ is $1.3 \%$. For the non-adaptive phase estimation protocol, $| {A}_{8}| =66$ but the total probability of obtaining a state in A₈ is roughly $79.2 \%$. These results indicate that although the adaptive phase estimation protocol outperforms the non-adaptive protocol in the fault-free case [3, 23] the adaptive phase estimation protocol cannot be used in the presence of measurement readout errors for four rounds, and has significantly fewer states in A₈ for eight rounds. Therefore in a noisy implementation of phase estimation, the non-adaptive protocol is preferable. A list of the states belonging to the acceptance set A₈ for both the adaptive and non-adaptive protocols can be found at https://github.com/godott/GKP_phase_estimation.git.

Suppose now that for m rounds the set A_m is not empty. We then know there exists a $\tilde{\delta }$ such that the protocol is a $(1,\tilde{\delta })$ fault-tolerant m round phase estimation protocol for measurement readout errors. However there are other fault locations where a single fault resulting in a qubit error could potentially lead to an uncorrectable shift error on output states. Note that an X error prior to applying the controlled-$D(\sqrt{2\pi })$ gate will do nothing since the qubit is in the $| +\rangle$ state. A Z error will just change the sign of one of the peaks of the output state in p space. However a damping event occurring on the ancilla qubit before or during the application of the controlled-$D(\sqrt{2\pi })$ gate can result in incorrectly applying the $D(\sqrt{2\pi })$ displacement on the cavity. If several damping events occur during the protocol, the output state in p space could potentially be badly deformed. In [3] it was proposed to replace the ancilla qubit by a k-qubit cat state so that if a single qubit undergoes damping, a single shift error of size $D\left(\tfrac{\sqrt{2\pi }}{k}\right)$ would occur which can be made small for large k. However this approach would require the fault-tolerant preparation of a large k-qubit cat state which would significantly increase the required resources in addition to adding substantially more locations where errors can occur. Similar issues were observed in [24] for preparing cat codes. However large cavity displacements were mitigated by interacting the cavity with a three level system $| f\rangle$, $| e\rangle$ and $| g\rangle$ instead of a qubit. A transition from $| f\rangle$ to $| e\rangle$ would not cause the incorrect gate from being applied. Thus two damping events would be required to cause the incorrect gate from being applied.

Here we propose an alternative solution for mitigating damping errors during the application of the controlled-$D(\sqrt{2\pi })$ gate which uses a single extra flag qubit [25–30]. Consider the circuit in figure 3. If a bit flip error occurs before or during the application of the controlled-$D(\sqrt{2\pi })$ gate, the flag qubit will be measured as 1 instead of 0. In such a case the strategy will be simply to abort the protocol and start anew. The analysis for a general damping event using the flag qubit is given in appendix A. It is shown that if ever a damping event on the $| +\rangle$ ancilla qubit causes the wrong $D(\sqrt{2\pi })$ gate to be applied, the flag qubit will be measured as 1 instead of 0 in which case we abort the protocol and start anew. Thus with the flag qubit, the phase estimation protocol is robust to a damping error during the application of the controlled displacement gate. Note however that if a fault occurs in addition to a damping event, the shift error resulting from the damping event can potentially go undetected. For instance, if in addition to damping, the flag qubit is subject to a measurement error, the measurement outcome can be 0 instead 1 resulting in acceptance when the protocol should have been aborted.

Zoom In Zoom Out Reset image size

Suppose that only a single damping event occurs during the four round phase estimation protocol, and that all other operations are fault-free. Using the analytic expressions derived in appendix A and for a damping rate of $p=0.5$, we performed simulations to numerically characterize the effects of qubit damping on the prepared GKP state. We found that the shift error resulting from a single damping event was negligible compared to the largest shift error resulting from a single measurement readout error.

Lastly, an X error prior to the ${\rm{\Lambda }}({{\rm{e}}}^{{\rm{i}}\gamma })$ gate will change the sign of γ. After propagating through the Hadamard gate, the X error will become a Z error which will not affect the measurement outcome of the ancilla qubit. For the four round non-adaptive phase estimation protocol, we performed a simulation which showed that the shift error resulting from a single X error prior to the ${\rm{\Lambda }}({{\rm{e}}}^{{\rm{i}}\gamma })$ gate is much smaller than the largest shift error arising from a measurement readout error on the ancilla qubit.

Let us assume that the noise model afflicting the oscillator can cause a shift error of size at most $\tilde{\delta }$ in p space (i.e. a shift of the form ${{\rm{e}}}^{-{\rm{i}}\tilde{\delta }q}$). From the above, the largest shift error that can arise from a single fault afflicting the qubit space during the four round fault-tolerant phase estimation protocol is $0.235\lt \sqrt{\pi }/6$. Following definition 5, we conclude that the protocol described in this section is a $(1,\tilde{\delta })$-fault-tolerant state preparation of an approximate logical $| \overline{0}\rangle$ GKP state with $\tilde{\delta }=\sqrt{\pi }/6-0.235\approx 0.06$.

In section 5 we will discuss a physical implementation of the phase estimation protocol using circuit QED. A much more detailed analysis of the noise afflicting the cavity and the resulting shift errors will be provided.

In section 5.2, the fault tolerant implementation of phase estimation described in this section is analyzed for a noise model which introduces gate noise, measurement readout errors and ancilla state preparation errors in addition to photon loss, amplitude damping and dephasing. Here we consider a noiseless implementation of the four round non-adaptive phase estimation protocol described in this section which will be used to benchmark the noisy implementation.

Since the shift correction obtained from non-adaptive phase estimation is in p space, for a range of values for δ, we computed using the integral in definition 1 after performing a shift back correction using equation (9). This allowed us to compute the probability of correcting shift errors in p of size $\tfrac{\sqrt{\pi }}{2}-\delta$. Plots for the states 0101 and 1000 which belong to A₄ (see table 1) are given in figure 4. The ${\left(\tfrac{\sqrt{\pi }}{2}-\delta ,\epsilon \right)}_{p}$ correctable properties of the other states in A₄ are similar to the ones shown in figure 4. It can be seen that for small values of size $\tfrac{\sqrt{\pi }}{2}-\delta$, both states can correct the shift error with probability close to 1.

Zoom In Zoom Out Reset image size

In this section we will describe a direct implementation of the controlled-$D(\sqrt{2\pi })$ gate. From [31, 32], the Hamiltonian describing the coupling between a qubit and a cavity in the dispersive regime is given by

$$\begin{eqnarray}&&{H}_{s}={\tilde{\omega }}_{r}{a}^{\dagger }a+{\tilde{\omega }}_{a}Z+\chi {a}^{\dagger }{aZ}-\phi {\left({a}^{\dagger }a\right)}^{2}Z,\end{eqnarray} \tag{ 13 }$$

where $\phi =\tfrac{{g}^{4}}{{{\rm{\Delta }}}^{3}}$, $\chi =\tfrac{{g}^{2}}{{\rm{\Delta }}}-\phi$, ${\tilde{\omega }}_{r}={\omega }_{r}+\phi$ and ${\tilde{\omega }}_{a}=\tfrac{{\omega }_{a}+\chi }{2}$. Here g is the coupling strength between the qubit and the cavity and ${\rm{\Delta }}=| {\omega }_{r}-{\omega }_{a}|$. The term $\phi {\left({a}^{\dagger }a\right)}^{2}Z$ corresponds to the nonlinear dispersive shift. The Hamiltonian in equation (13) can be derived by performing an exact diagonalization of the Jaynes–Cummings Hamiltonian and keeping only leading order terms in ϕ [32]. Note that a more systematic treatment of the qubit as an anharmonic oscillator leads to an additional term in equation (13) given by $-\tfrac{K}{2}{\left({a}^{\dagger }a\right)}^{2}$ which is referred to as the Kerr nonlinearity [33]. Hence, in our analysis, we choose the system Hamiltonian to be given by

$$\begin{eqnarray}&&{H}_{s}={\tilde{\omega }}_{r}{a}^{\dagger }a+{\tilde{\omega }}_{a}Z+\chi {a}^{\dagger }{aZ}-\phi {\left({a}^{\dagger }a\right)}^{2}Z-\displaystyle \frac{K}{2}{\left({a}^{\dagger }a\right)}^{2}.\end{eqnarray} \tag{ 14 }$$

The direct implementation of the controlled-$D(\sqrt{2\pi })$ gate can be achieved using the drive Hamiltonian

$$\begin{eqnarray}&&{H}_{d}(t)={ \mathcal E }(t){a}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\omega }_{d}t}+{{ \mathcal E }}^{* }(t){{a}{\rm{e}}}^{{\rm{i}}{\omega }_{d}t},\end{eqnarray} \tag{ 15 }$$

where ${ \mathcal E }(t)$ describes the pulse shape of the drive and ${\omega }_{d}$ is the drive frequency. Thus the total Hamiltonian describing the evolution of the qubit-cavity system during the implementation of the controlled-$D(\sqrt{2\pi })$ gate is given by

$$\begin{eqnarray}&&H(t)={H}_{s}+{H}_{d}(t).\end{eqnarray} \tag{ 16 }$$

Going into the rotating frame of the qubit and the cavity and defining ${\phi }_{\pm }=\phi \pm K$ and ${\omega }_{\pm }=\pm \chi$, in appendix B we show that

$$\begin{eqnarray}\begin{array}{rcl}{V}_{R}(0,T) & = & { \mathcal T }{{\rm{e}}}^{-{\rm{i}}{\displaystyle \int }_{0}^{T}H(t^{\prime} ){\rm{d}}{t}^{\prime} }\\ & = & {R}_{+}(T)D({A}_{+}){{\rm{e}}}^{{{\rm{i}}{B}}_{+}}| 0\rangle \langle 0| +{R}_{-}(T)D({A}_{-}){{\rm{e}}}^{{{\rm{i}}{B}}_{-}}| 1\rangle \langle 1| ,\end{array}\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}&&{R}_{\pm }(T)={{\rm{e}}}^{-{\rm{i}}{T}{\omega }_{\pm }{a}^{\dagger }a}\left(1\pm {\rm{i}}{T}{\phi }_{\pm }{\left({a}^{\dagger }a\right)}^{2}\right),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{A}_{\pm }=-{\rm{i}}{\int }_{0}^{T}{ \mathcal E }(t){{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }t}{\rm{d}}{t}\mp {\phi }_{\pm }(2{a}^{\dagger }a-1){\int }_{0}^{T}{ \mathcal E }(t){{t}{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }t}{\rm{d}}{t}+{ \mathcal O }({\phi }_{\pm }^{2}),\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}\begin{array}{rcl}{B}_{\pm } & = & -\displaystyle \frac{{\rm{i}}}{2}{\displaystyle \int }_{0}^{T}{\rm{d}}{t}{\displaystyle \int }_{t}^{T}{\rm{d}}{t}^{\prime} \left({{ \mathcal E }}^{* }(t^{\prime} ){ \mathcal E }(t){{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }(t^{\prime} -t)}-{ \mathcal E }(t^{\prime} ){{ \mathcal E }}^{* }(t){{\rm{e}}}^{-{\rm{i}}{\omega }_{\pm }(t^{\prime} -t)}\right)\\ & & \pm \displaystyle \frac{{\phi }_{\pm }(2{a}^{\dagger }a+1)}{2}{\displaystyle \int }_{0}^{T}{\rm{d}}{t}{\displaystyle \int }_{t}^{T}{\rm{d}}{t}^{\prime} \left({{ \mathcal E }}^{* }(t^{\prime} ){ \mathcal E }(t){{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }(t^{\prime} -t)}+{ \mathcal E }(t^{\prime} ){{ \mathcal E }}^{* }(t){{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }(t^{\prime} -t)}\right)(t^{\prime} -t)+{ \mathcal O }({\phi }_{\pm }^{2}).\end{array}\end{eqnarray} \tag{ 20 }$$

In deriving equations (18)–(20), we kept only leading order terms in ${\phi }_{\pm }$ since the Hamiltonian in equation (14) neglects higher order terms in ϕ and K. We keep only leading order terms since with current experimental parameter values, ${\phi }_{\pm }$ is roughly three orders of magnitude smaller than χ.

Firstly, notice that the terms ${R}_{\pm }(T)$ introduce a relative rotation between the $| 0\rangle$ and $| 1\rangle$ state of the qubit. Neglecting the terms proportional to ${\phi }_{\pm }$, it was shown in [3] that by choosing the total interaction time to be $T=\pi /\chi$, the relative rotation between $| 0\rangle$ and $| 1\rangle$ can be set to one. However due to the presence of the nonlinear dispersive shift and Kerr terms in equation (18), it is clear that the relative rotation cannot be completely removed. Further, ${R}_{\pm }(T)$ does not depend on the pulse shape ${ \mathcal E }(t)$ of the drive term. Therefore even with an optimized pulse shape (see below), the effects from the nonlinear dispersive shift and the Kerr will cause an undesired relative rotation between the qubit states. Regardless, we will still choose the interaction time T to be π / χ to ensure that we eliminate the relative rotation from the leading order terms in equation (18). In appendix B.4, we provide further details using the analytic expressions to show that to mitigate the effects due to the nonlinear dispersive shift and Kerr terms would require the protocol to take place on time scales of order $\tfrac{1}{{\phi }_{\pm }}$ which (for the parameter values in table 2) is four orders of magnitude larger than $\tfrac{\pi }{\chi }$. For such long time scales, noise terms such as photon loss, dephasing and damping would render the protocol impractical.

Table 2.  Parameter values for the two simulations of the non-adaptive noisy phase estimation protocol. The values for κ are chosen based on state of the art 3D cavities. The parameters in the second column results in a T₁ time given by ${T}_{1}=\tfrac{2\pi }{\kappa }=10\,\mathrm{ms}$. For a resonator frequency f_r = 7 GHz, the quality factor is given by $Q=\tfrac{{f}_{r}}{\kappa /(2\pi )}=7\times {10}^{7}$. For all simulations, the squeezing level of the input squeezed vacuum state is chosen to be 0.2.

Parameter values	Simulation 1	Simulation 2
p	5 × 10−3	10−3
$\tfrac{\kappa }{2\pi }$	1.59 × 10−5 MHz	1.59 × 10−6 MHz
$\tfrac{{\gamma }_{1}}{2\pi }$	1.06 × 10−3 MHz	1.06 × 10−4 MHz
$\tfrac{{\gamma }_{\phi }}{2\pi }$	$7.96\times {10}^{-4}$ MHz	7.96 × 10−5 MHz
g	8.92 MHz	5.09 MHz
Δ	0.32 GHz	0.32 GHz
K (Kerr)	10−4 MHz	10−5 MHz

The terms in equation (17) that perform the desired controlled displacement gate are given by A_±in equation (19). The goal is to choose a pulse shape that implements the desired gate while at the same time minimizes the contributions from the nonlinear dispersive shift and the Kerr term (terms proportional to ϕ_±in equation (19)). In order to be experimentally relevant, it is important to choose a pulse shape that is accessible to near term experiments using 2D and 3D cavities. We chose a Gaussian pulse with the following parameters

$$\begin{eqnarray}&&{ \mathcal E }(t)\approx -2.09562\left(\displaystyle \frac{\sqrt{2\pi }}{T}\right){{\rm{e}}}^{-\tfrac{{\left(t-\mu \right)}^{2}}{2{\sigma }^{2}}},\end{eqnarray} \tag{ 21 }$$

where $\mu =\tfrac{\pi }{2\chi }$ and $\sigma =\tfrac{\pi }{10\chi }$. With this Gaussian pulse we obtain

$$\begin{eqnarray}&&{A}_{\pm }\approx \mp \displaystyle \frac{\sqrt{2\pi }}{2}\mp {\phi }_{\pm }(2{a}^{\dagger }a-1){\int }_{0}^{T}{ \mathcal E }(t){{t}{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }t}{\rm{d}}{t}.\end{eqnarray} \tag{ 22 }$$

The amplitude of the Gaussian, given by $-2.09562\left(\tfrac{\sqrt{2\pi }}{T}\right)$, is chosen numerically to ensure that the appropriate gate is being performed.

Since both the A₊ and ${A}_{-}$ terms are symmetric with opposite signs, the gate $D\left(-\sqrt{\tfrac{\pi }{2}}\right)$ in figure 3 is not required. Furthermore, the first term in equation (22) is independent of χ. However the contribution from the second term in equation (22) will depend on the coupling strength and relative frequencies between the qubit and the cavity. Thus reducing the value of χ will minimize the undesired effects arising from the nonlinear dispersive shift and the Kerr term while still producing the desired gate.

The B_± terms of equation (20) result in a phase difference between the $| 0\rangle$ and $| 1\rangle$ qubit states. However, computing the first integrals in equation (20) (i.e. terms before ϕ_±), we find that for our chosen pulse, the phase difference remains fixed during every round of the phase estimation protocol. Therefore after applying the controlled displacement gate, we apply an additional phase gate which removes the phases introduced by the left most integrals of B_±. Note however that the ϕ_±terms in B_±will not be canceled and will thus introduce additional shift errors.

In this section we perform a numerical analysis for the noisy implementation of the phase estimation protocol described in section 4.2. The simulation is performed in several steps that we now describe. First, the controlled displacement gate is modeled using the following master equation

$$\begin{eqnarray}&&\dot{\rho }=-{\rm{i}}[H(t),\rho ]+\kappa { \mathcal D }[a]\rho +{\gamma }_{1}{ \mathcal D }[{\sigma }_{-}]\rho +{\gamma }_{\phi }{ \mathcal D }[{\sigma }_{z}]\rho ,\end{eqnarray} \tag{ 23 }$$

where H(t) is given by equation (16) after going into the rotating frame and ${ \mathcal D }[L]\rho =(2L\rho {L}^{\dagger }-{L}^{\dagger }L\rho -\rho {L}^{\dagger }L)/2$. The density matrix corresponds to the joint state of the cavity, ancilla qubit and flag qubit. The parameters κ, γ₁ and γ_ϕ correspond to the photon loss rate, the qubit decay rate and qubit dephasing. The pulse shape of the drive is given by equation (21). The total evolution time of the controlled displacement is given by $\tfrac{\pi }{\chi }$ with $\chi =\tfrac{{g}^{2}}{{\rm{\Delta }}}-\phi$.

Both before and after the controlled displacement gate, the state of the cavity is subject to photon loss and its evolution is computed by solving a master equation. Based on current gate, state preparation and measurement times, we chose the evolution time from the preparation of the ancilla to the first CNOT gate (after which the controlled-displacement gate is performed) to be 0.14 μs. Similarly, the evolution time after the controlled-displacement gate to the time the ancilla is measured is also chosen to be 0.14 μs. Thus during one round, the cavity freely evolves with photon loss for a total of 0.28 μs. In addition, before and after the controlled displacement gate, we allow all qubit locations to fail with the following depolarizing noise model

• 1.  
  With probability p, each two-qubit gate is followed by a two-qubit Pauli error drawn uniformly and independently from $\{I,X,Y,Z\}{}^{\otimes 2}\setminus \{I\otimes I\}$.
• 2.  
  With probability $\tfrac{2p}{3}$, the preparation of the $| 0\rangle$ state is replaced by $| 1\rangle =X| 0\rangle$.
• 3.  
  With probability $\tfrac{2p}{3}$, any single qubit measurement has its outcome flipped.
• 4.  
  With probability p/10, each Hadamard gate is followed by a Pauli error drawn uniformly and independently from { X, Y, Z }.

We chose p/10 for Hadamard gate failures since for current superconducting architectures, single qubit gate fidelities are about an order of magnitude higher than two-qubit gate fidelities. In practice, the gate errors applied to the qubit Hilbert space will depend strongly on the circuit QED architecture and should also be modeled using a master equation as in equation (23). However, we chose a depolarizing model in order to reduce the computation time and simplify the analysis. We also mention that in our analysis we assumed that g is tunable. Thus both before and after the controlled-displacement gate, the cavity and qubit system is decoupled and can be treated separately.

The master equation was numerically solved using Qutip [34], our code can be accessed at https://github.com/godott/GKP_phase_estimation.git. Due to the long computation time during the controlled displacement gate, performing a full Monte Carlo simulation to take into account gate errors with the depolarizing model was unfeasible (i.e. gate locations both before and after the controlled displacement gate). Instead, we analytically computed the error probabilities for each Pauli operator (using the depolarizing noise model described above) immediately after the first CNOT gate, before both measurement locations and before the phase gate¹⁰ . At each location, all possible Pauli operators based on their associated probabilities were added. For a given Pauli error, the probabilities (which are expressed as functions of p) were then used (in addition to the state evolution obtained from the master equation) to compute the final probability of obtaining a given output state. Instances with two or more faults occurring on the qubit Hilbert space before and after the controlled-displacement gate were neglected. However our analysis is still more complete than previous implementations which considered only measurement errors and errors during the controlled-displacement gate. More details of the Pauli simulation can be found in appendix D.

We performed two different simulations where for each simulation, the chosen parameter values are given in table 2. The parameters chosen for the first simulation (middle column in table 2) are based on current experimental values for 2D and 3D cavities [24, 35–39]. The parameters chosen for the second simulation are based on values that might be obtained with improved future technologies. Plots showing the probability of correcting shift errors $\tfrac{\sqrt{\pi }}{2}-\delta$ after performing a shift correction of the output states of the noisy phase estimation protocol using equation (9) are given in figure 5. In what follows we will refer to these plots as $\epsilon -\delta$ plots. Note that noise during the phase estimation protocol introduces more shift errors in p space. Therefore in figure 5, only  − δ plots in p space are shown.

Zoom In Zoom Out Reset image size

For the four round phase estimation protocol, we find that the state 1010 has a similar  − δ plot to the one for the state 0101 whereas the state 0010 has a similar  − δ plot to the one for the state 1000. The probability of obtaining each state for the two noisy simulations (parameters in table 2) are given in table 3. It can be seen that the probability of the state 0101 to correct shift errors is significantly more affected by the noise than the state 1000 or 1111. To understand why this is the case, it is useful to look at the wave function densities for these three states in the q basis when no noise is present (see figure 6). Comparing the wave function densities, it can be observed that the state 1000 has three peaks with similar amplitudes, whereas the states 0101 and 1111 have two smaller peaks compared to the peak at the center. Performing a numerical analysis, it is found that when only damping is present (with all other noise terms including Kerr and nonlinear dispersive shift set to zero), damping has a negligible effect on the height of the peaks for all considered states. However performing a numerical simulation with only the Kerr and nonlinear dispersive shift terms present, these terms significantly reduce the height of the peaks for the state 0101 and 1010. Therefore, these states are left with one large peak in q with all other peaks close to zero which results in a wave function density with very low resolution in p space (which thus affects the ability for these states to correct shift errors in p). Interestingly, the state 1111 is more robust to the Kerr and nonlinear dispersive shift contributions since there is a much smaller reduction in the amplitudes of the smaller peaks compared to those for the states 0101 and 1010 (see figure 7). Furthermore, the states 1000 and 0010 have three large peaks in q and thus even with a reduced amplitude, these states have a higher wave function density resolution in p space.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we will derive the unitary operator describing the evolution of the qubit-cavity state during the application of the control-displacement gate. In the dispersive regime, the qubit-cavity interaction can be described as:

$$\begin{eqnarray}&&H(t)={H}_{s}+{H}_{d}(t),\end{eqnarray} \tag{ B1 }$$

where

$$\begin{eqnarray}&&{H}_{s}={\tilde{\omega }}_{r}{a}^{\dagger }a+{\tilde{\omega }}_{a}Z+\chi {a}^{\dagger }{aZ}-\phi {\left({a}^{\dagger }a\right)}^{2}Z.\end{eqnarray} \tag{ B2 }$$

and

$$\begin{eqnarray}&&{H}_{d}(t)={ \mathcal E }(t)(a+{a}^{\dagger }).\end{eqnarray} \tag{ B3 }$$

In equation (B2), $\phi =\tfrac{{g}^{4}}{{{\rm{\Delta }}}^{3}}$, ${\tilde{\omega }}_{r}={\omega }_{r}+\phi$, $\chi =\tfrac{{g}^{2}}{{\rm{\Delta }}}-\phi$ and ${\tilde{\omega }}_{a}=\tfrac{{\omega }_{a}+\chi }{2}$. The term $\phi {\left({a}^{\dagger }a\right)}^{2}Z$ corresponds to the nonlinear dispersive shift. Note that we have not included the Kerr term $-\tfrac{K}{2}{\left({a}^{\dagger }a\right)}^{2}$. At the end of this section we will modify our results to take into account its effect. We also note that in writing the drive Hamiltonian in equation (B3), we neglected a term of the form λX (where λ = (g/Δ)) and all higher powers in λ. For parameter values considered in this paper, we found the effect of this term to be negligible. Additionally, we chose a pulse shape which is real valued.

In equation (B3), we represent the drive pulse ${ \mathcal E }(t)$ as

$$\begin{eqnarray}&&{ \mathcal E }(t)={{\rm{\Omega }}}_{x}(t)\cos ({\omega }_{d}t)+{{\rm{\Omega }}}_{y}(t)\sin ({\omega }_{d}t),\end{eqnarray} \tag{ B4 }$$

where ${\omega }_{d}={\tilde{\omega }}_{r}-\chi$ is the drive frequency.

Since the Hamiltonian does not commute at different times, the unitary operator describing the time evolution under the Hamiltonian of equation (B1) is given by

$$\begin{eqnarray}&&V(0,T)={ \mathcal T }{{\rm{e}}}^{-{\rm{i}}{\int }_{0}^{T}H(t^{\prime} ){\rm{d}}{t}^{\prime} }.\end{eqnarray} \tag{ B5 }$$

The right-hand side of equation (B5) can be computed using the Suzuki–Trotter decomposition, so that

$$\begin{eqnarray}\begin{array}{rcl}V(0,T) & = & \left(\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \prod _{j=1}^{n}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}({\tilde{\omega }}_{r}+(\chi -\phi {a}^{\dagger }a)Z){a}^{\dagger }a}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{ \mathcal E }({t}_{j})(a+{a}^{\dagger })}\right){{\rm{e}}}^{-{\rm{i}}{T}{\tilde{\omega }}_{a}Z}\\ & = & \left(\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \prod _{j=1}^{n}{{\rm{e}}}^{-{\rm{i}}\tilde{t}({\omega }_{+}-\phi {a}^{\dagger }a){a}^{\dagger }a}{{\rm{e}}}^{-{\rm{i}}\tilde{t}{ \mathcal E }({t}_{j})(a+{a}^{\dagger })}\right){{\rm{e}}}^{-{\rm{i}}{\tilde{\omega }}_{a}T}| 0\rangle \langle 0| \\ & & +\left(\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \prod _{j=1}^{n}{{\rm{e}}}^{-{\rm{i}}\tilde{t}({\omega }_{-}+\phi {a}^{\dagger }a){a}^{\dagger }a}{{\rm{e}}}^{-{\rm{i}}\tilde{t}{ \mathcal E }({t}_{j})(a+{a}^{\dagger })}\right){{\rm{e}}}^{{\rm{i}}{\tilde{\omega }}_{a}T}| 1\rangle \langle 1| ,\end{array}\end{eqnarray} \tag{ B6 }$$

where $\tilde{t}\equiv T/n$, ${\omega }_{\pm }\equiv \tilde{{\omega }_{r}}\pm \chi$.

Now, the two terms on the right-hand side of equation (B6) can be decomposed as

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \prod _{j=1}^{n}{{\rm{e}}}^{-{\rm{i}}\tilde{t}({\omega }_{+}-\phi {a}^{\dagger }a){a}^{\dagger }a}{{\rm{e}}}^{-{\rm{i}}\tilde{t}{ \mathcal E }({t}_{j})(a+{a}^{\dagger })}=\mathop{\mathrm{lim}}\limits_{n\to \infty }{R}_{n}^{n}\left({R}_{n}^{-n}{D}_{{t}_{n}}{R}_{n}^{n}\right)\cdots \left({R}_{n}^{-2}{D}_{{t}_{2}}{R}_{n}^{2}\right)\left({R}_{n}^{-1}{D}_{{t}_{1}}{R}_{n}^{1}\right),\end{eqnarray} \tag{ B7 }$$

where

$$\begin{eqnarray}&&{R}_{n}={{\rm{e}}}^{-{\rm{i}}\tilde{t}({\omega }_{+}-\phi {a}^{\dagger }a){a}^{\dagger }a},\end{eqnarray} \tag{ B8 }$$

and

$$\begin{eqnarray}&&{D}_{{t}_{j}}={{\rm{e}}}^{-{\rm{i}}\tilde{t}{ \mathcal E }({t}_{j})(a+{a}^{\dagger })}.\end{eqnarray} \tag{ B9 }$$

Hence from the above, we see that we need to compute terms of the form

$$\begin{eqnarray}&&{R}_{n}^{-k}{D}_{{t}_{k}}{R}_{n}^{k}={{\rm{e}}}^{{\rm{i}}\tilde{t}k\left({\omega }_{+}-\phi {a}^{\dagger }a\right){a}^{\dagger }a}{{\rm{e}}}^{-{\rm{i}}\tilde{t}{ \mathcal E }({t}_{k})\left(a+{a}^{\dagger }\right)}{{\rm{e}}}^{-{\rm{i}}\tilde{t}k\left({\omega }_{+}-\phi {a}^{\dagger }a\right){a}^{\dagger }a}.\end{eqnarray} \tag{ B10 }$$

In order to compute the products appearing in equation (B10), we will use the identity ${{A}{e}}^{B}{A}^{-1}={{e}}^{{{ABA}}^{-1}}$. Defining

$$\begin{eqnarray}&&H^{\prime} =k({\omega }_{+}-\phi {a}^{\dagger }a){a}^{\dagger }a,\end{eqnarray} \tag{ B11 }$$

and

$$\begin{eqnarray}&&{A}_{\tilde{t}}={{\rm{e}}}^{-{\rm{i}}\tilde{t}{ \mathcal E }({t}_{k})(a+{a}^{\dagger })},\end{eqnarray} \tag{ B12 }$$

with

$$\begin{eqnarray}&&A(\tilde{t})={{\rm{e}}}^{{\rm{i}}{H}^{\prime} \tilde{t}}{A}_{\tilde{t}}{{\rm{e}}}^{-{\rm{i}}{H}^{\prime} \tilde{t}}\end{eqnarray} \tag{ B13 }$$

we have that

$$\begin{eqnarray}&&A(\tilde{t})=\exp \left(-{\rm{i}}\tilde{t}{ \mathcal E }({t}_{k}){{\rm{e}}}^{{\rm{i}}{H}^{\prime} \tilde{t}}(a+{a}^{\dagger }){{\rm{e}}}^{-{\rm{i}}{H}^{\prime} \tilde{t}}\right),\end{eqnarray} \tag{ B14 }$$

The term ${{\rm{e}}}^{{\rm{i}}{H}^{\prime} \tilde{t}}(a+{a}^{\dagger }){{\rm{e}}}^{-{\rm{i}}{H}^{\prime} \tilde{t}}$ in equation (B14) can be computed using Heisenberg's equation of motion. We obtain

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{a}(\tilde{t})}{{\rm{d}}\tilde{t}}={\rm{i}}[H^{\prime} ,a(\tilde{t})],\end{eqnarray} \tag{ B15 }$$

with

$$\begin{eqnarray}&&[H^{\prime} ,a]=k(\phi ^{\prime} +2\phi {a}^{\dagger }a)a,\end{eqnarray} \tag{ B16 }$$

where we defined

$$\begin{eqnarray}&&\phi ^{\prime} \equiv \phi -{\omega }_{+}.\end{eqnarray} \tag{ B17 }$$

We thus have the following set of coupled differential equations:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{a}(\tilde{t})}{{\rm{d}}\tilde{t}}={\rm{i}}{k}\left(\phi ^{\prime} +2\phi {a}^{\dagger }(\tilde{t})a(\tilde{t})\right)a(\tilde{t}),\end{eqnarray} \tag{ B18 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{a}}^{\dagger }(\tilde{t})}{{\rm{d}}\tilde{t}}=-{{\rm{i}}{k}{a}}^{\dagger }(\tilde{t})(\phi ^{\prime} +2\phi {a}^{\dagger }(\tilde{t})a(\tilde{t})),\end{eqnarray} \tag{ B19 }$$

The solutions to equations (B18) and (B19) are given by

$$\begin{eqnarray}&&a(\tilde{t})={{a}{\rm{e}}}^{{\rm{i}}{k}\tilde{t}\left(\phi ^{\prime} +2\phi ({a}^{\dagger }a-1)\right)},\end{eqnarray} \tag{ B20 }$$

and

$$\begin{eqnarray}&&{a}^{\dagger }(\tilde{t})={{\rm{e}}}^{-{\rm{i}}{k}\tilde{t}\left(\phi ^{\prime} +2\phi ({a}^{\dagger }a-1)\right)}{a}^{\dagger }.\end{eqnarray} \tag{ B21 }$$

To see this, we take a derivative of equation (B20) to obtain

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{a}(\tilde{t})}{{\rm{d}}\tilde{t}}={\rm{i}}{k}\phi ^{\prime} a(\tilde{t})+2{\rm{i}}{k}\phi a(\tilde{t})({a}^{\dagger }a-1).\end{eqnarray} \tag{ B22 }$$

Comparing equations (B18) and (B22), we must have that

$$\begin{eqnarray}&&a(\tilde{t})({a}^{\dagger }a-1)={a}^{\dagger }(\tilde{t}){a}^{2}(\tilde{t}).\end{eqnarray} \tag{ B23 }$$

Using equations (B20) and (B21) and defining ${H}_{{te}}\equiv -{kt}\left(\phi ^{\prime} +2\phi ({a}^{\dagger }a-1)\right)$, we have that

$$\begin{eqnarray}\begin{array}{rcl}{a}^{\dagger }(\tilde{t}){a}^{2}(\tilde{t}) & = & {{\rm{e}}}^{{{\rm{i}}{H}}_{{te}}\tilde{t}}{a}^{\dagger }{{\rm{e}}}^{-{{\rm{i}}{H}}_{{te}}\tilde{t}}{{\rm{e}}}^{{{\rm{i}}{H}}_{{te}}\tilde{t}}{{a}{\rm{e}}}^{-{{\rm{i}}{H}}_{{te}}\tilde{t}}{{a}{\rm{e}}}^{-{{\rm{i}}{H}}_{{te}}\tilde{t}}\\ & = & {a}^{\dagger }{{\rm{e}}}^{-2{\rm{i}}\phi k\tilde{t}}{{a}{\rm{e}}}^{2{\rm{i}}\phi k\tilde{t}}{{a}{\rm{e}}}^{-{{\rm{i}}{H}}_{{te}}\tilde{t}}\\ & = & {a}^{\dagger }{{a}{a}{\rm{e}}}^{-{{\rm{i}}{H}}_{{te}}\tilde{t}}\\ & = & {a}^{\dagger }{aa}(\tilde{t}),\end{array}\end{eqnarray} \tag{ B24 }$$

as desired. A similar calculation can be done for ${a}^{\dagger }(\tilde{t})$. Hence we have that

$$\begin{eqnarray}&&A(\tilde{t})=\exp \left(-{\rm{i}}\tilde{t}{ \mathcal E }({t}_{k})\left({{a}{\rm{e}}}^{{\rm{i}}{{\rm{\Phi }}}_{k}\tilde{t}}+{{\rm{e}}}^{-{\rm{i}}{{\rm{\Phi }}}_{k}\tilde{t}}{a}^{\dagger }\right)\right),\end{eqnarray} \tag{ B25 }$$

for

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{k}=k\left(\phi ^{\prime} +2\phi ({a}^{\dagger }a-1)\right).\end{eqnarray} \tag{ B26 }$$

Writing V(0, T) (equation (B6)) as

$$\begin{eqnarray}&&V(0,T)={V}_{+}(0,T)| 0\rangle \langle 0| +{V}_{-}(0,T)| 1\rangle \langle 1| ,\end{eqnarray} \tag{ B27 }$$

using equation (B25) and the fact that ${R}_{n}^{n}={{\rm{e}}}^{-{\rm{i}}{T}\left({\omega }_{+}-\phi {a}^{\dagger }a\right){a}^{\dagger }a}\equiv {R}_{+}(T)$, we can write

$$\begin{eqnarray}&&{V}_{+}(0,T)={R}_{+}(T)\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \prod _{k=1}^{n}\exp \left\{{A}_{k}{a}^{\dagger }-{{aA}}_{k}^{\dagger }\right\},\end{eqnarray} \tag{ B28 }$$

where we defined

$$\begin{eqnarray}&&{A}_{k}\equiv -\displaystyle \frac{{\rm{i}}{T}}{n}{ \mathcal E }({t}_{k}){{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}.\end{eqnarray} \tag{ B29 }$$

Notice that we can write V₊(0, T) as products of displacements D(A_k). However, now A_k is an operator instead of a complex number. In order to compute the products in (B28), we will need to obtain a relation for terms of the form D(A_k)D(A_j) (where A_k is given in equation (B29)). Since A_k is proportional to T/n and remembering that we will take the limit where $n\to \infty$, using Baker–Campbell–Hausdorff, we have the exact relation

$$\begin{eqnarray}&&\exp \left\{{A}_{k}{a}^{\dagger }-{{aA}}_{k}^{\dagger }+{A}_{j}{a}^{\dagger }-{{aA}}_{j}\right\}=D({A}_{k})D({A}_{j})\exp \left\{-\displaystyle \frac{1}{2}[{A}_{k}{a}^{\dagger }-{{aA}}_{k}^{\dagger },{A}_{j}{a}^{\dagger }-{{aA}}_{j}]\right\}.\end{eqnarray} \tag{ B30 }$$

The commutator on the right-hand side of equation (B30) has four terms

$$\begin{eqnarray}&&\left[{A}_{k}{a}^{\dagger }-{{aA}}_{k}^{\dagger },{A}_{j}{a}^{\dagger }-{{aA}}_{j}\right]=\left[{A}_{k}{a}^{\dagger },{A}_{j}{a}^{\dagger }\right]-\left[{A}_{k}{a}^{\dagger },{{aA}}_{j}^{\dagger }\right]-\left[{{aA}}_{k}^{\dagger },{A}_{j}{a}^{\dagger }\right]+\left[{{aA}}_{k}^{\dagger },{{aA}}_{j}^{\dagger }\right].\end{eqnarray} \tag{ B31 }$$

As we will show, two of these terms vanish when taking the limit $n\to \infty$.

First, we compute

$$\begin{eqnarray}&&\left[{A}_{k}{a}^{\dagger },{A}_{j}{a}^{\dagger }\right]={\left(\displaystyle \frac{-{\rm{i}}{T}}{n}\right)}^{2}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j})\left[{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger },{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }\right],\end{eqnarray} \tag{ B32 }$$

with

$$\begin{eqnarray}&&\begin{array}{l}\left[{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger },{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }\right]={{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger }{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }-{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger }\\ \ \ \ =\ {{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger }{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}{a}^{\dagger }{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k+j}}-{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}{a}^{\dagger }{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}.\end{array}\end{eqnarray} \tag{ B33 }$$

Now, terms such as ${{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}$ can be computed by defining

$$\begin{eqnarray}&&{H}_{\mathrm{temp}}\equiv -j\phi ^{\prime} \left(1+\tilde{\delta }({a}^{\dagger }a-1)\right),\end{eqnarray} \tag{ B34 }$$

and invoking Heisenberg's equation of motion. Using $[{H}_{\mathrm{temp}},{a}^{\dagger }]=-j\phi ^{\prime} \tilde{\delta }{a}^{\dagger }$ and solving the differential equation, we obtain

$$\begin{eqnarray}&&{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}={{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}j\phi ^{\prime} \tilde{\delta }}{a}^{\dagger }.\end{eqnarray} \tag{ B35 }$$

Using the result of equation (B35) into equation (B33), we have

$$\begin{eqnarray}\begin{array}{rcl}\left[{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger },{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }\right] & = & {{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}k\phi ^{\prime} \tilde{\delta }}{a}^{\dagger }{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}(k+j)\phi ^{\prime} \tilde{\delta }}{a}^{\dagger }{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}-{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}j\phi ^{\prime} \tilde{\delta }}{a}^{\dagger }{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}(j+k)\phi ^{\prime} \tilde{\delta }}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}\\ & = & {{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}(j+k)\phi ^{\prime} \tilde{\delta }}\left({{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}k\phi ^{\prime} \tilde{\delta }}-{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}j\phi ^{\prime} \tilde{\delta }}\right){\left({a}^{\dagger }\right)}^{2}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}.\end{array}\end{eqnarray} \tag{ B36 }$$

Using $\phi ^{\prime} \tilde{\delta }=2\phi$ and expanding equation (B36) to leading order in ϕ, we obtain

$$\begin{eqnarray}\begin{array}{rcl}\left[{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger },{{\rm{e}}}^{-\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }\right] & \approx & \left(1-\displaystyle \frac{2{\rm{i}}{T}}{n}(j+k)\phi \right)\left(\displaystyle \frac{-2{\rm{i}}{T}}{n}\phi \right){\left({a}^{\dagger }\right)}^{2}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}\\ & \approx & \displaystyle \frac{-2{\rm{i}}{T}}{n}(k-j)\phi {\left({a}^{\dagger }\right)}^{2}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}.\end{array}\end{eqnarray} \tag{ B37 }$$

Inserting equation (B37) into equation (B32), we obtain

$$\begin{eqnarray}&&\left[{A}_{k}{a}^{\dagger },{A}_{j}{a}^{\dagger }\right]=2{\left(\displaystyle \frac{-{\rm{i}}{T}}{n}\right)}^{3}(k-j)\phi { \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){\left({a}^{\dagger }\right)}^{2}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}+{ \mathcal O }({\phi }^{2}).\end{eqnarray} \tag{ B38 }$$

A similar calculation shows that

$$\begin{eqnarray}&&\left[{{aA}}_{k}^{\dagger },{{aA}}_{j}^{\dagger }\right]=2{\left(\displaystyle \frac{{\rm{i}}{T}}{n}\right)}^{3}(k-j)\phi { \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}{a}^{2}+{ \mathcal O }({\phi }^{2}).\end{eqnarray} \tag{ B39 }$$

Next we compute the cross terms. We first have

$$\begin{eqnarray}&&\left[{A}_{k}{a}^{\dagger },{{aA}}_{j}^{\dagger }\right]={\left(\displaystyle \frac{{\rm{i}}{T}}{n}\right)}^{2}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j})\left[{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger },{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}\right],\end{eqnarray} \tag{ B40 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}\left[{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger },{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}\right] & = & {{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger }{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}-{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger }\\ & = & {{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger }{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}-{{a}{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{a}^{\dagger }\\ & = & {a}^{\dagger }{{a}{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}-{{a}{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{a}^{\dagger }{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}\\ & = & \left({a}^{\dagger }a-{{aa}}^{\dagger }{{\rm{e}}}^{\tfrac{-2{\rm{i}}{T}}{n}\phi (k-j)}\right){{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}\\ & = & \left({a}^{\dagger }a\left(1-{{\rm{e}}}^{\tfrac{-2{\rm{i}}{T}}{n}\phi (k-j)}\right)-{{\rm{e}}}^{\tfrac{-2{\rm{i}}{T}}{n}\phi (k-j)}\right){{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}.\end{array}\end{eqnarray} \tag{ B41 }$$

Expanding the term proportional to ${a}^{\dagger }a$ to leading order in ϕ, equation (B41) becomes

$$\begin{eqnarray}&&\left[{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger },{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}\right]\approx \left(\displaystyle \frac{2{\rm{i}}{T}}{n}\phi {a}^{\dagger }a(k-j)-{{\rm{e}}}^{\tfrac{-2{\rm{i}}{T}}{n}\phi (k-j)}\right){{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}.\end{eqnarray} \tag{ B42 }$$

Inserting equation (B42) into equation (B40), we obtain

$$\begin{eqnarray}&&\left[{A}_{k}{a}^{\dagger },{{aA}}_{j}^{\dagger }\right]\approx 2{\left(\displaystyle \frac{{\rm{i}}{T}}{n}\right)}^{3}\phi (k-j){ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){a}^{\dagger }{{a}{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}-{\left(\displaystyle \frac{{\rm{i}}{T}}{n}\right)}^{2}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){{\rm{e}}}^{\tfrac{-2{\rm{i}}{T}}{n}\phi (k-j)}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}.\end{eqnarray} \tag{ B43 }$$

The last commutator to compute is

$$\begin{eqnarray}&&\left[{{aA}}_{k}^{\dagger },{A}_{j}{a}^{\dagger }\right]={\left(\displaystyle \frac{-{\rm{i}}{T}}{n}\right)}^{2}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j})\left[{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}},{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}{a}^{\dagger }}\right],\end{eqnarray} \tag{ B44 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}\left[{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}},{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }\right] & = & {{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{a}^{\dagger }-{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{j}}{a}^{\dagger }{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}\\ & = & {{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{a}^{\dagger }-{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{a}^{\dagger }{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}{{a}{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}\\ & = & {{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}\left({a}^{\dagger }a\left({{\rm{e}}}^{\tfrac{2{\rm{i}}{T}}{n}\phi (k-j)}-1\right)+{{\rm{e}}}^{\tfrac{2{\rm{i}}{T}}{n}\phi (k-j)}\right)\\ & \approx & {{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}\left(\displaystyle \frac{2{\rm{i}}{T}}{n}\phi (k-j){a}^{\dagger }a+{{\rm{e}}}^{\tfrac{2{\rm{i}}{T}}{n}\phi (k-j)}\right)\end{array}\end{eqnarray} \tag{ B45 }$$

so that

$$\begin{eqnarray}&&\left[{{aA}}_{k}^{\dagger },{A}_{j}{a}^{\dagger }\right]\approx 2{\left(\displaystyle \frac{{\rm{i}}{T}}{n}\right)}^{3}\phi (k-j){ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{a}^{\dagger }a+{\left(\displaystyle \frac{-{\rm{i}}{T}}{n}\right)}^{2}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{{\rm{e}}}^{\tfrac{2{\rm{i}}{T}}{n}\phi (k-j)}\end{eqnarray} \tag{ B46 }$$

Inserting equations (B38), (B39), (B43) and (B46) into equation (B30), we have

$$\begin{eqnarray}\begin{array}{rcl}D({A}_{k}+{A}_{j}) & = & D({A}_{k})D({A}_{j})\exp \left[-\displaystyle \frac{1}{2}\left\{2{\left(\displaystyle \frac{-{\rm{i}}{T}}{n}\right)}^{3}(k-j)\phi { \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){\left({a}^{\dagger }\right)}^{2}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}\right.\right.\\ & & +2{\left(\displaystyle \frac{{\rm{i}}{T}}{n}\right)}^{3}(k-j)\phi { \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k+j)}}{a}^{2}-2{\left(\displaystyle \frac{{\rm{i}}{T}}{n}\right)}^{3}\phi (k-j){ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){a}^{\dagger }{{a}{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}\\ & & +{\left(\displaystyle \frac{{\rm{i}}{T}}{n}\right)}^{2}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){{\rm{e}}}^{\tfrac{-2{\rm{i}}{T}}{n}\phi (k-j)}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}-2{\left(\displaystyle \frac{{\rm{i}}{T}}{n}\right)}^{3}\phi (k-j){ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{a}^{\dagger }a\\ & & \left.\left.-{\left(\displaystyle \frac{-{\rm{i}}{T}}{n}\right)}^{2}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}{{\rm{e}}}^{\tfrac{2{\rm{i}}{T}}{n}\phi (k-j)}\right\}\right].\end{array}\end{eqnarray} \tag{ B47 }$$

Since the product D(A_n) D(A_n − 1) ⋯ D(A₂)D(A₁) will produce sums in the exponents proportional to n², all terms in equation (B47) proportional to ${\left(\tfrac{T}{n}\right)}^{3}$ will vanish in the limit where $n\to \infty$. Hence equation (B47) simplifies to

$$\begin{eqnarray}\begin{array}{rcl}D({A}_{k}+{A}_{j}) & = & D({A}_{k})D({A}_{j})\exp \left[{\rm{i}}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){\left(\displaystyle \frac{T}{n}\right)}^{2}\left\{\displaystyle \frac{{{\rm{e}}}^{\tfrac{-2{\rm{i}}{T}}{n}\phi (k-j)}{{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}-{{\rm{e}}}^{\tfrac{2{\rm{i}}{T}}{n}\phi (k-j)}{{\rm{e}}}^{\tfrac{{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{(k-j)}}}{2{\rm{i}}}\right\}\right]\\ & = & D({A}_{k})D({A}_{j})\exp \left[-{\rm{i}}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){\left(\displaystyle \frac{T}{n}\right)}^{2}\sin \left(\displaystyle \frac{T}{n}{\overline{{\rm{\Phi }}}}_{(k-j)}\right)\right],\end{array}\end{eqnarray} \tag{ B48 }$$

where ${\overline{{\rm{\Phi }}}}_{(k-j)}$ is defined as

$$\begin{eqnarray}&&{\overline{{\rm{\Phi }}}}_{(k-j)}\equiv (k-j)\left(\phi -{\omega }_{+}+2\phi {a}^{\dagger }a\right).\end{eqnarray} \tag{ B49 }$$

We thus have that the product of the modified displacement operators is given by

$$\begin{eqnarray}&&D({A}_{k})D({A}_{j})=D({A}_{k}+{A}_{j})\exp \left[{\rm{i}}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){\left(\displaystyle \frac{T}{n}\right)}^{2}\sin \left(\displaystyle \frac{T}{n}{\overline{{\rm{\Phi }}}}_{(k-j)}\right)\right].\end{eqnarray} \tag{ B50 }$$

We now wish to compute

$$\begin{eqnarray}&&{P}_{n}\equiv D({A}_{n})D({A}_{n-1})D({A}_{n-2})\cdots D({A}_{2})D({A}_{1}).\end{eqnarray} \tag{ B51 }$$

First notice that for any operators A and B, to leading order e^Ae^B = e^Be^Ae^{[A, B]}. Now for terms of the form $\exp \left[{\rm{i}}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){\left(\tfrac{T}{n}\right)}^{2}\sin \left(\tfrac{T}{n}{\overline{{\rm{\Phi }}}}_{(k-j)}\right)\right]D({\sum }_{m}{A}_{m})$, the commutator of the two exponents will be proportional to ${\left(\tfrac{T}{n}\right)}^{3}$ which will vanish in the limit where $n\to \infty$. Hence, we can commute all terms $\exp \left[{\rm{i}}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j}){\left(\tfrac{T}{n}\right)}^{2}\sin \left(\tfrac{T}{n}{\overline{{\rm{\Phi }}}}_{(k-j)}\right)\right]$ to the right-hand side of P_n. Hence we have

$$\begin{eqnarray}&&{P}_{n}=D\left(\displaystyle \sum _{k=1}^{n}{A}_{k}\right)\exp \left[{\rm{i}}\displaystyle \sum _{k\lt j}^{n}\displaystyle \frac{{T}^{2}}{{n}^{2}}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j})\sin \left(\displaystyle \frac{T}{n}{\overline{{\rm{\Phi }}}}_{(k-j)}\right)\right].\end{eqnarray} \tag{ B52 }$$

Now, taking the limit where $n\to \infty$ of P_n, we obtain

$$\begin{eqnarray}\begin{array}{rcl}\mathop{\mathrm{lim}}\limits_{n\to \infty }{P}_{n} & = & D\left(\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \sum _{k=1}^{n}\left(\displaystyle \frac{-{\rm{i}}{T}}{n}\right){ \mathcal E }({t}_{k}){{\rm{e}}}^{\tfrac{-{\rm{i}}{T}}{n}{{\rm{\Phi }}}_{k}}\right)\exp \left[{\rm{i}}\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \sum _{k\lt j}^{n}\displaystyle \frac{{T}^{2}}{{n}^{2}}{ \mathcal E }({t}_{k}){ \mathcal E }({t}_{j})\sin \left(\displaystyle \frac{T}{n}{\overline{{\rm{\Phi }}}}_{(k-j)}\right)\right]\\ & = & D\left(-{\rm{i}}{\displaystyle \int }_{0}^{T}{ \mathcal E }(t){{\rm{e}}}^{-{\rm{i}}{\rm{\Phi }}t}{\rm{d}}{t}\right)\exp \left[-{\rm{i}}{\displaystyle \int }_{0}^{T}{\rm{d}}{t}{\displaystyle \int }_{t}^{T}{\rm{d}}{t}^{\prime} { \mathcal E }(t){ \mathcal E }(t^{\prime} )\sin \left(\overline{{\rm{\Phi }}}(t^{\prime} -t)\right)\right]\end{array}\end{eqnarray} \tag{ B53 }$$

We conclude that

$$\begin{eqnarray}&&{V}_{+}(0,T)={R}_{+}(T)D({A}_{+}){{\rm{e}}}^{{{\rm{i}}{B}}_{+}},\end{eqnarray} \tag{ B54 }$$

where

$$\begin{eqnarray}&&{R}_{+}(T)={{\rm{e}}}^{-{\rm{i}}{T}({\omega }_{+}-\phi {a}^{\dagger }a){a}^{\dagger }a},\end{eqnarray} \tag{ B55 }$$

$$\begin{eqnarray}&&{A}_{+}=-{\rm{i}}{\int }_{0}^{T}{ \mathcal E }(t){{\rm{e}}}^{-{\rm{i}}{\rm{\Phi }}t}{\rm{d}}{t},\end{eqnarray} \tag{ B56 }$$

and

$$\begin{eqnarray}&&{B}_{+}=-{\int }_{0}^{T}{\rm{d}}{t}{\int }_{t}^{T}{\rm{d}}{t}^{\prime} { \mathcal E }(t){ \mathcal E }(t^{\prime} )\sin \left(\overline{{\rm{\Phi }}}(t^{\prime} -t)\right).\end{eqnarray} \tag{ B57 }$$

Hence to conclude, the unitary evolution of the Hamiltonian described in equation (B1) is given by

$$\begin{eqnarray}&&V(0,T)={R}_{+}(T)D({A}_{+}){{\rm{e}}}^{{{\rm{i}}{B}}_{+}}{{\rm{e}}}^{-{\rm{i}}{\tilde{\omega }}_{a}T}| 0\rangle \langle 0| +{R}_{-}(T)D({A}_{-}){{\rm{e}}}^{{{\rm{i}}{B}}_{-}}{{\rm{e}}}^{{\rm{i}}{\tilde{\omega }}_{a}T}| 1\rangle \langle 1| .\end{eqnarray} \tag{ B58 }$$

Recall that in deriving the expression for V₊(0, T), in several steps of the calculation we expanded to leading order in ϕ. Hence the expressions obtained in equations (B55)–(B57) are only valid to leading order in ϕ. Hence to leading order in ϕ, we have

$$\begin{eqnarray}&&{R}_{\pm }(T)={{\rm{e}}}^{-{\rm{i}}{T}{\omega }_{\pm }{a}^{\dagger }a}\left(1\pm {\rm{i}}{T}\phi {\left({a}^{\dagger }a\right)}^{2}\right)+{ \mathcal O }({\phi }^{2}),\end{eqnarray} \tag{ B59 }$$

$$\begin{eqnarray}&&{A}_{\pm }=-{\rm{i}}{\int }_{0}^{T}{ \mathcal E }(t){{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }t}{\rm{d}}{t}\mp \phi (2{a}^{\dagger }a-1){\int }_{0}^{T}{ \mathcal E }(t){{t}{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }t}{\rm{d}}{t}+{ \mathcal O }({\phi }^{2}),\end{eqnarray} \tag{ B60 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{B}_{\pm } & = & {\int }_{0}^{T}{\rm{d}}{t}{\int }_{t}^{T}{\rm{d}}{t}^{\prime} { \mathcal E }(t){ \mathcal E }(t^{\prime} )\sin \left({\omega }_{\pm }(t^{\prime} -t)\right)\pm \phi (2{a}^{\dagger }a+1){\int }_{0}^{T}{\rm{d}}{t}\\ & & \times {\int }_{t}^{T}{\rm{d}}{t}^{\prime} { \mathcal E }(t){ \mathcal E }(t^{\prime} )\cos \left({\omega }_{\pm }(t^{\prime} -t)\right)(t^{\prime} -t)+{ \mathcal O }({\phi }^{2}).\end{array}\end{eqnarray} \tag{ B61 }$$

When including the Kerr nonlinearity, the Hamiltonian during the control-displacement gate is now given by

$$\begin{eqnarray}&&H(t)={\tilde{\omega }}_{r}{a}^{\dagger }a+{\tilde{\omega }}_{a}Z+\chi {a}^{\dagger }{aZ}-\phi {\left({a}^{\dagger }a\right)}^{2}Z-\displaystyle \frac{K}{2}{\left({a}^{\dagger }a\right)}^{2}.\end{eqnarray} \tag{ B62 }$$

In this case, we have that

$$\begin{eqnarray}\begin{array}{rcl}V(0,T) & = & \left(\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \prod _{j=1}^{n}{{\rm{e}}}^{-{\rm{i}}\tilde{t}\left({\omega }_{+}-(\phi +\tfrac{K}{2}){a}^{\dagger }a\right){a}^{\dagger }a}{{\rm{e}}}^{-{\rm{i}}\tilde{t}{ \mathcal E }({t}_{j})(a+{a}^{\dagger })}\right){{\rm{e}}}^{-{\rm{i}}{\tilde{\omega }}_{a}T}| 0\rangle \langle 0| \\ & & +\left(\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \prod _{j=1}^{n}{{\rm{e}}}^{-{\rm{i}}\tilde{t}\left({\omega }_{-}+(\phi -\tfrac{K}{2}){a}^{\dagger }a\right){a}^{\dagger }a}{{\rm{e}}}^{-{\rm{i}}\tilde{t}{ \mathcal E }({t}_{j})(a+{a}^{\dagger })}\right){{\rm{e}}}^{{\rm{i}}{\tilde{\omega }}_{a}T}| 1\rangle \langle 1| ,\end{array}\end{eqnarray} \tag{ B63 }$$

Hence, we see that the analysis of section B.1 is exactly the same, with $\phi \to \phi \pm \tfrac{K}{2}$ in R_± (T), A_±and B_±.

In this section, we consider the same Hamiltonian as in equation (B1) but with a drive term of the form

$$\begin{eqnarray}&&{H}_{d}(t)={ \mathcal E }(t){a}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\omega }_{d}t}+{{ \mathcal E }}^{* }(t){{a}{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}.\end{eqnarray} \tag{ B64 }$$

We will go into the rotating frame of both the qubit and the cavity. We define

$$\begin{eqnarray}&&{H}_{s}^{{\prime} }={\tilde{\omega }}_{r}{a}^{\dagger }a+{\tilde{\omega }}_{a}Z.\end{eqnarray} \tag{ B65 }$$

With $U(t)={{\rm{e}}}^{{{\rm{i}}{H}}_{s}^{{\prime} }t}$ and applying the transformation $U(t)H(t){U}^{-1}(t)-{\rm{i}}{U}(t)\tfrac{\partial }{\partial t}{U}^{-1}(t)\equiv {H}_{R}(t)$ to the Hamiltonian in equation (B1) with H_d(t) given in equation (B64), we obtain

$$\begin{eqnarray}&&{H}_{R}(t)=\chi {{Za}}^{\dagger }a-\phi {\left({a}^{\dagger }a\right)}^{2}Z+{ \mathcal E }(t){a}^{\dagger }{{\rm{e}}}^{{\rm{i}}\left({\tilde{\omega }}_{r}-{\omega }_{d}\right)t}+{{ \mathcal E }}^{* }(t){{a}{\rm{e}}}^{-{\rm{i}}\left({\tilde{\omega }}_{r}-{\omega }_{d}\right)t}\end{eqnarray} \tag{ B66 }$$

The frequency dependence in the drive term of equation (B66) can be eliminated by choosing ${\omega }_{d}={\tilde{\omega }}_{r}$. Again, we wish to compute the unitary evolution

$$\begin{eqnarray}&&{V}_{R}(0,T)={ \mathcal T }{{\rm{e}}}^{-{\rm{i}}{\int }_{0}^{T}{H}_{R}(t^{\prime} ){\rm{d}}{t}^{\prime} }.\end{eqnarray} \tag{ B67 }$$

Using the Suzuki–Trotter decomposition

$$\begin{eqnarray}\begin{array}{rcl}{V}_{R}(0,T) & = & \left(\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \prod _{j=1}^{n}{{\rm{e}}}^{-{\rm{i}}\tilde{t}({\omega }_{+}-\phi {a}^{\dagger }a){a}^{\dagger }a}{{\rm{e}}}^{-{\rm{i}}\tilde{t}({ \mathcal E }({t}_{j}){a}^{\dagger }+{{ \mathcal E }}^{* }({t}_{j})a)}\right)| 0\rangle \langle 0| \\ & & +\left(\mathop{\mathrm{lim}}\limits_{n\to \infty }\displaystyle \prod _{j=1}^{n}{{\rm{e}}}^{-{\rm{i}}\tilde{t}({\omega }_{-}+\phi {a}^{\dagger }a){a}^{\dagger }a}{{\rm{e}}}^{-{\rm{i}}\tilde{t}\left({ \mathcal E }({t}_{j}){a}^{\dagger }+{{ \mathcal E }}^{* }({t}_{j})a\right)}\right)| 1\rangle \langle 1| ,\end{array}\end{eqnarray} \tag{ B68 }$$

where now we have that ω_± = ±χ.

Comparing equation (B68) to equation (B6), we see that we can follow the same steps as in appendix B.1 by using the new values for ω_±and replacing ${ \mathcal E }$ by ${{ \mathcal E }}^{* }$ in the conjugate expressions. Doing so, we obtain

$$\begin{eqnarray}&&{V}_{R}(0,T)={R}_{+}(T)D({A}_{+}){{\rm{e}}}^{{{\rm{i}}{B}}_{+}}| 0\rangle \langle 0| +{R}_{-}(T)D({A}_{-}){{\rm{e}}}^{{{\rm{i}}{B}}_{-}}| 1\rangle \langle 1| ,\end{eqnarray} \tag{ B69 }$$

where

$$\begin{eqnarray}&&{R}_{\pm }(T)={{\rm{e}}}^{-{\rm{i}}{T}{\omega }_{\pm }{a}^{\dagger }a}\left(1\pm {\rm{i}}{T}\phi {\left({a}^{\dagger }a\right)}^{2}\right),\end{eqnarray} \tag{ B70 }$$

$$\begin{eqnarray}&&{A}_{\pm }=-{\rm{i}}{\int }_{0}^{T}{ \mathcal E }(t){{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }t}{\rm{d}}{t}\mp \phi (2{a}^{\dagger }a-1){\int }_{0}^{T}{ \mathcal E }(t){{t}{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }t}{\rm{d}}{t}+{ \mathcal O }({\phi }^{2}),\end{eqnarray} \tag{ B71 }$$

$$\begin{eqnarray}\begin{array}{rcl}{B}_{\pm } & = & -\displaystyle \frac{{\rm{i}}}{2}{\displaystyle \int }_{0}^{T}{\rm{d}}{t}{\displaystyle \int }_{t}^{T}{\rm{d}}{t}^{\prime} \left({{ \mathcal E }}^{* }(t^{\prime} ){ \mathcal E }(t){{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }(t^{\prime} -t)}-{ \mathcal E }(t^{\prime} ){{ \mathcal E }}^{* }(t){{\rm{e}}}^{-{\rm{i}}{\omega }_{\pm }(t^{\prime} -t)}\right)\\ & & \pm \displaystyle \frac{\phi (2{a}^{\dagger }a+1)}{2}{\displaystyle \int }_{0}^{T}{\rm{d}}{t}{\displaystyle \int }_{t}^{T}{\rm{d}}{t}^{\prime} \left({{ \mathcal E }}^{* }(t^{\prime} ){ \mathcal E }(t){{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }(t^{\prime} -t)}+{ \mathcal E }(t^{\prime} ){{ \mathcal E }}^{* }(t){{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }(t^{\prime} -t)}\right)(t^{\prime} -t)+{ \mathcal O }({\phi }^{2}).\end{array}\end{eqnarray} \tag{ B72 }$$

In this section we will show that regardless of the chosen pulse shape (even with numerical tools such as optimal control), to leading order in ϕ_±, the changes to the unitary evolution of the qubit-cavity system due to the nonlinear dispersive shift and Kerr terms cannot be completely removed for time scales $T\ll \tfrac{1}{{\phi }_{\pm }}$. Using equations (B55), (B60) and (B61), we can write (for instance choosing the terms affecting the $| 0\rangle \langle 0|$)

$$\begin{eqnarray}&&{R}_{+}(T)D({A}_{+})={{\rm{e}}}^{-{\rm{i}}{\omega }_{+}{{Ta}}^{\dagger }a}{{\rm{e}}}^{{\rm{i}}{\omega }_{+}T{\phi }_{+}{\left({a}^{\dagger }a\right)}^{2}}{{\rm{e}}}^{{I}_{1}{a}^{\dagger }-{I}_{1}^{* }a-{I}_{2}{\phi }_{+}\left(2{a}^{\dagger }a-1\right){a}^{\dagger }+{I}_{2}^{* }{\phi }_{+}a\left(2{a}^{\dagger }a-1\right)},\end{eqnarray} \tag{ B73 }$$

where we define

$$\begin{eqnarray}&&{I}_{1}\equiv -{\rm{i}}{\int }_{0}^{T}{ \mathcal E }(t){{\rm{e}}}^{{\rm{i}}{\omega }_{+}t}{\rm{d}}{t},\end{eqnarray} \tag{ B74 }$$

and

$$\begin{eqnarray}&&{I}_{2}\equiv {\int }_{0}^{T}{ \mathcal E }(t){{t}{\rm{e}}}^{{\rm{i}}{\omega }_{+}t}{\rm{d}}{t}.\end{eqnarray} \tag{ B75 }$$

Let $A={\rm{i}}{\omega }_{+}T{\phi }_{+}{\left({a}^{\dagger }a\right)}^{2}$ and $B={I}_{1}{a}^{\dagger }-{I}_{1}^{* }a-{I}_{2}{\phi }_{+}(2{a}^{\dagger }a-1){a}^{\dagger }+{I}_{2}^{* }{\phi }_{+}a(2{a}^{\dagger }a-1)$. Using the Baker–Campbell–Hausdorff lemma and keeping only leading order terms in ϕ₊, we have

$$\begin{eqnarray}&&{e}^{A}{e}^{B}={e}^{A+B+\tfrac{1}{2}[A,B]-\tfrac{1}{12}[B[A,B]]-\tfrac{1}{720}\left[B,[B,[B,[B,A]]]\right]}.\end{eqnarray} \tag{ B76 }$$

Computing the commutators, we find

$$\begin{eqnarray}&&\begin{array}{l}\,{e}^{A}{e}^{B}\,=\,{e}^{\displaystyle \frac{1}{6}i{\omega }_{+}T{\phi }_{+}\left(-\displaystyle \frac{1}{5}| {I}_{1}{| }^{4}+| {I}_{1}{| }^{2}\right)}{e}^{\left[\left({I}_{2}{\phi }_{+}+\left(1+\displaystyle \frac{1}{2}i{\omega }_{+}T{\phi }_{+}\right){I}_{1}\right){a}^{\dagger }-\left({I}_{2}^{\dagger }{\phi }_{+}+\left(1-\displaystyle \frac{1}{2}i{\omega }_{+}T{\phi }_{+}\right){I}_{1}^{\dagger }\right)a+\displaystyle \frac{1}{6}i{\omega }_{+}T{\phi }_{+}({I}_{1}^{2}{\left({a}^{\dagger }\right)}^{2}+4| {I}_{1}{| }^{2}{a}^{\dagger }a+{\left({I}^{\dagger }\right)}^{2}{a}^{2})\right.}\\ \,{}^{\left.+i{\omega }_{+}T{\phi }_{+}({I}_{1}{a}^{\dagger }{a}^{\dagger }a+{I}_{1}^{\dagger }{a}^{\dagger }{aa})-2{I}_{2}{\phi }_{+}{a}^{\dagger }{{aa}}^{\dagger }+2{I}_{2}^{\dagger }{\phi }_{+}{{aa}}^{\dagger }a+{\omega }_{+}T{\phi }_{+}{\left({a}^{\dagger }a\right)}^{2}\right]}\end{array}\end{eqnarray} \tag{ B77 }$$

Notice that the last term in equation (B77) is proportional to ${\left({a}^{\dagger }a\right)}^{2}$ and independent of the pulse shape (performing the same calculation as above including the B₊ term will not change the conclusion). The terms dependent on the pulse shape are expressed as lower powers of a and ${a}^{\dagger }$. One cannot eliminate all terms proportional to ϕ₊ unless one chooses a time scale comparable to $\tfrac{1}{{\phi }_{+}}$. For time scales on the order of $\tfrac{1}{{\phi }_{+}}$, higher order terms in ϕ₊ will be relevant and therefore it might possible to choose pulse (with numerical techniques such as optimal control) which can eliminate effects from the nonlinear dispersive shift and Kerr terms However using the parameters of table 2, $\tfrac{1}{{\phi }_{+}}$ is roughly four orders of magnitude longer than the chosen time scale ($T=\tfrac{\pi }{\chi }$) of our protocol. For such long time scales, effects due to photon loss, damping and dephasing would render the protocol impractical.

From figure 2, the operator describing the evolution of a single round of phase estimation in terms of the measurement operator set $\left\{{M}_{x}=\tfrac{D\left(-\sqrt{2\pi }\right)+{\left(-1\right)}^{x}{{\rm{e}}}^{{\rm{i}}\gamma }D\left(\sqrt{2\pi }\right)}{2}| x=0,1\right\}$ can be written as

$$\begin{eqnarray}&&{U}_{\mathrm{measure}}^{\gamma }=\left\{| 0\rangle \langle 0| \otimes {M}_{0}+| 1\rangle \langle 1| \otimes {M}_{1}\right\}\end{eqnarray} \tag{ C1 }$$

After the measurement, the state become

$$\begin{eqnarray}\begin{array}{rcl}| {\psi }_{\mathrm{output}}\rangle & = & \displaystyle \frac{{U}_{\mathrm{measure}}^{\gamma }| {\psi }_{\mathrm{input}}\rangle }{\sqrt{{{\Pr }}_{x}}}\\ & = & \displaystyle \frac{1}{2\sqrt{{{\Pr }}_{x}}}\left(D\left(-\sqrt{2\pi }\right)| {\psi }_{\mathrm{input}}\rangle +{\left(-1\right)}^{x}{{\rm{e}}}^{{\rm{i}}\gamma }D\left(\sqrt{2\pi }\right)| {\psi }_{\mathrm{input}}\rangle \right)\\ & = & \displaystyle \frac{1}{2\sqrt{{{\Pr }}_{x}}}\left(D\left(-\sqrt{2\pi }\right)| {\psi }_{\mathrm{input}}\rangle +{{\rm{e}}}^{{\rm{i}}(\gamma +x\pi )}D\left(\sqrt{2\pi }\right)| {\psi }_{\mathrm{input}}\rangle \right)\end{array}\end{eqnarray} \tag{ C2 }$$

where Pr_x is the probability of measuring x ∈ {0, 1}.

After M arounds of the phase estimation circuit, the following state will be generated:

$$\begin{eqnarray}&&| {\rm{\Phi }}(x[M])\rangle =\displaystyle \frac{1}{{2}^{M}\sqrt{N}}\displaystyle \prod _{m=1}^{M}\left[D\left(-\sqrt{\displaystyle \frac{\pi }{2}}\right)+{{\rm{e}}}^{{\rm{i}}({\gamma }_{m}+{x}_{m}\pi )}D\left(\sqrt{\displaystyle \frac{\pi }{2}}\right)\right]| \mathrm{sq}.\mathrm{vac}\rangle ,\end{eqnarray} \tag{ C3 }$$

where $N={\prod }_{i}{{\Pr }}_{{x}_{i}}$ is the normalization factor. From the product of M sums, if we choose j of them to be positive shifts $\left(D\left(\sqrt{\tfrac{\pi }{2}}\right)\right)$, M − j to be negative shifts $\left(D\left(-\sqrt{\tfrac{\pi }{2}}\right)\right)$, we can produce a peak at $\sqrt{\pi /2}(j-(M-j))=\sqrt{2\pi }(j-M/2)$ in q space, and the produced phase is ${\prod }_{l\in {S}_{j}}{{\rm{e}}}^{{\rm{i}}({\gamma }_{l}+{x}_{l}\pi )}$ (S_j is the set of rounds we choose to be positive shifts). There are $\left(\begin{array}{c}M\\ m\end{array}\right)$ ways to choose such an S_j, by choosing any of such S_j we can produce a peak at $\sqrt{2\pi }(j-M/2)$. Thus, $| {\rm{\Phi }}(x[M])\rangle$ can be written as:

$$\begin{eqnarray}&&| {\rm{\Phi }}(x[M])\rangle =\displaystyle \frac{1}{{2}^{M}\sqrt{N}}\displaystyle \sum _{j=0}^{M}{c}_{j}(x[M])D\left(\sqrt{2\pi }(j-M/2)\right)| \mathrm{sq}.\mathrm{vac}\rangle ,\end{eqnarray} \tag{ C4 }$$

where $N\approx {\sum }_{j=0}^{M}| {c}_{j}(x[M]){| }^{2}$ with ${c}_{j}(x[M])={\sum }_{\{{S}_{j}\}}{\prod }_{l\in {S}_{j}}{{\rm{e}}}^{{\rm{i}}({\gamma }_{l}+{x}_{l}\pi )}$.

Assume we already performed M − 1 rounds of phase estimation and the measurement results ${x}_{1},...,{x}_{M-1}$ and angle parameters we chose (denoted as ${\gamma }_{1},\,\ldots {\gamma }_{M-1}$) are known. Then from (C4), the input state to the Mth round is (including the qubit):

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{{\rm{input}}}{\rangle }_{M}=\displaystyle \frac{1}{{2}^{M-1}\sqrt{{{\Pr }}_{{x}_{1}}...{{\Pr }}_{{x}_{M-1}}}}\displaystyle \sum _{j=0}^{M-1}{c}_{j}(x[M-1])D\left(\sqrt{2\pi }(j-(M-1)/2)\right)| 0\rangle \otimes | \mathrm{sq}.\mathrm{vac}\rangle .\end{eqnarray} \tag{ C5 }$$

The probability of measuring x_M is thus given by

$$\begin{eqnarray}\begin{array}{rcl}{{\Pr }}_{{x}_{M}} & = & \langle {{\rm{\Phi }}}_{\mathrm{input}}| {\left({U}_{\mathrm{measure}}^{\gamma }\right)}^{\dagger }{U}_{\mathrm{measure}}^{\gamma }| {{\rm{\Phi }}}_{\mathrm{input}}{\rangle }_{M}\\ & = & {\left|\displaystyle \frac{1}{{2}^{M}\sqrt{{{\Pr }}_{{x}_{1}}...{{\Pr }}_{{x}_{M-1}}}}\displaystyle \sum _{j=0}^{M}{c}_{j}(x[M])D\left(\sqrt{2\pi }(j-(M)/2)\right)| \mathrm{sq}.\mathrm{vac}\rangle \right|}^{2}.\end{array}\end{eqnarray} \tag{ C6 }$$

Notice that the transfer probability from one peak to another is almost 0 (< 10⁻³⁵) assuming an initial squeezing level of 0.2. Therefore, we can make the following approximation

$$\begin{eqnarray}&&\approx \displaystyle \frac{1}{{2}^{2M}{{\Pr }}_{{x}_{1}}...{{\Pr }}_{{x}_{M-1}}}\displaystyle \sum _{j=0}^{M}{\left|{c}_{j}(x[M])D(\sqrt{2\pi }(j-(M)/2))| \mathrm{sq}.\mathrm{vac}\rangle \right|}^{2}\end{eqnarray} \tag{ C7 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{{2}^{2M}{{\Pr }}_{{x}_{1}}...{{\Pr }}_{{x}_{M-1}}}\displaystyle \sum _{j=0}^{M}| {c}_{j}(x[M]){| }^{2}\mathop{\underbrace{\langle \mathrm{sq}.\mathrm{vac}| {D}^{\dagger }(\sqrt{2\pi }(j-(M)/2))| D(\sqrt{2\pi }(j-(M)/2))| \mathrm{sq}.\mathrm{vac}\rangle }}\limits_{1}\end{eqnarray} \tag{ C8 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{{2}^{2M}{{\Pr }}_{{x}_{1}}...{{\Pr }}_{{x}_{M-1}}}\displaystyle \sum _{j=0}^{M}| \displaystyle \sum _{\{{S}_{j}\}}\displaystyle \prod _{l\in {S}_{j}}{{\rm{e}}}^{{\rm{i}}({\gamma }_{l}+{x}_{l}\pi )}{| }^{2}\end{eqnarray} \tag{ C9 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{{2}^{2M}{{\Pr }}_{{x}_{1}}...{{\Pr }}_{{x}_{M-1}}}\displaystyle \sum _{j=0}^{M}\left[\left(\displaystyle \sum _{\{{S}_{j}^{m}\}}\displaystyle \prod _{l\in {S}_{j}^{m}}{{\rm{e}}}^{-{\rm{i}}({\gamma }_{l}+{x}_{l}\pi )}\right)\left(\displaystyle \sum _{\{{S}_{j}^{n}\}}\displaystyle \prod _{k\in {S}_{j}^{n}}{{\rm{e}}}^{{\rm{i}}({\gamma }_{k}+{x}_{k}\pi )}\right)\right].\end{eqnarray} \tag{ C10 }$$

The above expression is recursive, hence substituting ${{\Pr }}_{{x}_{M-1}}$ into it, we find

$$\begin{eqnarray}&&{{\Pr }}_{{x}_{M}}\approx \displaystyle \frac{{\displaystyle \sum }_{j=0}^{M}\left[\left({\displaystyle \sum }_{\{{S}_{j}^{m}\}}{\displaystyle \prod }_{l\in {S}_{j}^{m}}{{\rm{e}}}^{-{\rm{i}}({\gamma }_{l}+{x}_{l}\pi )}\right)\left({\displaystyle \sum }_{\{{S}_{j}^{n}\}}{\displaystyle \prod }_{k\in {S}_{j}^{n}}{{\rm{e}}}^{{\rm{i}}({\gamma }_{k}+{x}_{k}\pi )}\right)\right]}{4{\displaystyle \sum }_{j=0}^{M-1}\left[\left({\displaystyle \sum }_{\{{S}_{j}^{m}\}}{\displaystyle \prod }_{l\in {S}_{j}^{m}}{{\rm{e}}}^{-{\rm{i}}({\gamma }_{l}+{x}_{l}\pi )}\right)\left({\displaystyle \sum }_{\{{S}_{j}^{n}\}}{\displaystyle \prod }_{k\in {S}_{j}^{n}}{{\rm{e}}}^{{\rm{i}}({\gamma }_{k}+{x}_{k}\pi )}\right)\right]}.\end{eqnarray} \tag{ C11 }$$