Fast Simulation of Bosonic Qubits via Gaussian Functions in Phase Space

Bosonic qubits are a promising route to building fault-tolerant quantum computers on a variety of physical platforms. Studying the performance of bosonic qubits under realistic gates and measurements is challenging with existing analytical and numerical tools. We present a novel formalism for simulating classes of states that can be represented as linear combinations of Gaussian functions in phase space. This formalism allows us to analyze and simulate a wide class of non-Gaussian states, transformations, and measurements. We demonstrate how useful classes of bosonic qubits—Gottesman-Kitaev-Preskill (GKP), cat, and Fock states—can be simulated using this formalism, opening the door to investigating the behavior of bosonic qubits under Gaussian channels and measurements, non-Gaussian transformations such as those achieved via gate teleportation, and important non-Gaussian measurements such as threshold and photon-number detection. Our formalism enables simulating these situations with levels of accuracy that are not feasible with existing methods. Finally, we use a method informed by our formalism to simulate circuits critical to the study of fault-tolerant quantum computing with bosonic qubits but beyond the reach of existing techniques. Specifically, we examine how finite-energy GKP states transform under realistic qubit phase gates; interface with a continuous-variable cluster state; and transform under non-Clifford t gate teleportation using magic states. We implement our simulation method as a part of the open-source Strawberry Fields python library.

Popular Summary

Bosonic qubits refer to encoding a qubit, i.e., a two-dimensional quantum system, into a continuous-variable (CV) quantum system like a pulse of light. Bosonic qubits are a promising path to building a quantum computer: first, they decrease the physical resources required to perform a computation, as qubit preparation, gates, and measurement may correspond to experimentally accessible operations; second, they are capable of correcting the inevitable physical errors that arise during a computation.

Motivated by the quest to build noise-resistant quantum computers that employ bosonic qubits, our paper provides mathematical and numerical tools to simulate bosonic qubits under a wide class of useful gates, measurements, and noise sources. The key to our technique is to represent the state of the CV quantum system as a linear combination of Gaussian functions, with the evolution of the state during a computation corresponding to transformations of the means, covariance matrices, and coefficients of the Gaussian functions. With our technique, we perform novel simulations of a leading type of bosonic qubit known as Gottesman-Kitaev-Preskill states undergoing realistic, noisy gates. We provide our simulator in the user-friendly Strawberry Fields software library.

Our results will be useful for optimizing the construction of device components used to prepare and manipulate bosonic qubits. Moreover, our results will inform more sophisticated models for predicting the behavior of scalable quantum computers employing millions of bosonic qubits.

Figure 1

A circuit representation of the Choi-Jamiolkowki method for applying the most general Gaussian CP map to a state. Here, $|0⟩$ denotes the vacuum state; the vertical line linking two modes with ellipses on the ends indicates the two-mode squeezing operation (in this case, we take the limit as squeezing goes to infinity); and $\hat{\mathrm{\Phi }}(\cdot )$ is a general Gaussian CP map applied to one half of the Choi state so that it can be teleported onto the target state $\hat{\rho }$. Since the two-mode squeezing operation is a Gaussian transformation, this circuit can be treated with conditional Gaussian dynamics when $\hat{\rho }$ is a Gaussian state.

Figure 2

Error syndrome measurement with normalizable GKP states following the approach from Ref. [11]. First, shifts in $q$ are corrected: an encoded data qubit ${\left|\psi \right.⟩}_{\mathrm{GKP}}$ and an ancilla ${\left|+\right.⟩}_{\mathrm{GKP}}$ are sent through a sum gate, and ${\left|\psi \right.⟩}_{\mathrm{GKP}}$ is displaced according to the result of a homodyne $q$ measurement on the ancilla, mod $\sqrt{\pi \hslash }$. A similar procedure follows for shifts in $p$.

Figure 3

Gaussian decomposition of the Wigner function of the (non-Gaussian) cat state. We use the labels $\mathcal{M}=\{+,-,z,\overline{z}\}$ for the four Gaussians needed. Note that the last two terms have complex coefficients and means and are also complex conjugates of each other; thus, the imaginary parts of the Gaussians cancel out and we only plot the real parts. The details of this particular decomposition are provided in Sec. 4b. For an introductory tutorial implemented by the authors on using the bosonic backend of Strawberry Fields to obtain this figure, see Ref. [34].

Figure 4

Heralding scheme to do (a) single- and (b) many-photon additions. The vertical line with two ellipses at the end is used to denote a two-mode squeezing operation. Note that we always postselect on a single click for all but the first mode.

Figure 5

Infidelity between the Gaussian approximation of a Fock state ${\hat{\rho }}_n$ with Wigner function described in Eqs. (61a) and an exact Fock state $|n⟩$ as a function of the squeezing parameter $r$. The fidelity is defined as $F=⟨n|{\hat{\rho }}_n|n⟩$.

Figure 6

Optical implementation of the GKP qubit t gate up to global phase, following the method from Ref. [11]. Here, in the ideal limit, ${|M⟩}_{\mathrm{GKP}}=(1/\sqrt{2})(e^{-i\pi /8}{|0⟩}_{\mathrm{GKP}}+e^{i\pi /8}{|1⟩}_{\mathrm{GKP}})$. The GKP Hadamard gate (with operator $\hat{H}$) is a $\pi /2$ phase shift, and the GKP $\mathrm{cz}$ gate is discussed in Appendix pp8. The output $p_0$ of the $p$-homodyne measurement is processed by the function $s(p_0)$, which rounds the value to the nearest $n\sqrt{\pi \hslash }$ and returns the parity of $n$. If $n$ is even (odd), no phase gate [a phase gate $\hat{P}(s)=e^{is\hslash {\hat{q}}^2/2}$ with $s=1$] is applied.

Figure 7

Simulation of two cat states passing through a lossy beam splitter, with the Wigner function of the first mode displayed for increasing values of $\alpha$. Note the difference in scale between $q$ and $p$. We compare the output of the simulation using the fock backend of Strawberry Fields with a cutoff of 50 photons (top row), and using a linear combination of Gaussians in the bosonic backend (bottom row). While the number of photons in the state increases with $\alpha$, saturating the cutoff for the Fock basis simulation and leading to an incorrect output, the number of Gaussian peaks in phase space that require tracking remains constant. Note additionally that the mixing introduced by loss necessitates the full density matrix being tracked in the Fock basis, while loss does not increase the complexity of simulation for the linear combination of Gaussians.

Figure 8

Wigner functions for a ${|{+}^{ϵ}⟩}_{\mathrm{GKP}}$ state with $ϵ=0.01$ (20 dB of per-peak squeezing) teleported onto a 15 dB squeezed state. We entangle the states via a CV $\mathrm{cz}$ gate and subject them to 1% loss in each mode; next, we measure the mode originally containing the GKP state in the $p$ basis, postselecting on $p=0$. In the top panel, we perform the simulation using the fock backend in Strawberry Fields with a cutoff of 50 photons. In the bottom row, we employ the bosonic backend. We see that the bosonic backend can maintain a correct description of the state even as the number of photons becomes very high due to the large per-peak squeezing.

Figure 9

Wigner functions for GKP, cat, and Fock states. All figures are produced using linear combinations of Gaussians in phase space, as presented in Sec. 4. For the GKP state, negative regions correspond directly to Gaussian peaks with negative coefficients, while for the cat state, negativity is produced by the sinusoidal oscillations of complex-valued Gaussian peaks. In the single-photon case, negativity is produced by the difference of two zero-mean Gaussians with slightly different covariance matrices.

Figure 10

Two thousand homodyne measurement samples in both quadratures for a cat (top row) and a GKP (bottom row) state obtained using Algorithm 19. For the probability density functions (PDFs), we compare histograms obtained from finite sampling to the asymptotic distributions, obtained easily from the linear combination of the marginals of each Gaussian function in the Wigner function, as in Eq. (29). For an intermediate tutorial implemented by the authors on how these results can be generated with the bosonic backend of Strawberry Fields, see Ref. [35].

Figure 11

For different levels of target squeezing on vacuum $r_{\mathrm{target}}$, a comparison of the states produced from ideal squeezing (red), the average map from measurement-based squeezing (green), and a single-shot occurrence of measurement-based squeezing (blue). The center and axis lengths are set by the mean and variances of the output state. We simulate using an ancillary squeezed state of $r=1.2$ (approximately $10.5$ dB), and a homodyne detection efficiency of 0.99. A key takeaway is that the single-shot instances have less antisqueezing and a distribution of locations for where they can be centered in phase space, indicating that the level of antisqueezing and the correct mean value observed in the average map are a result of the averaging procedure itself.

Figure 12

(a) Wigner function for ${|+i^{ϵ}⟩}_{\mathrm{GKP}}$ simulated directly with $ϵ=0.1$ (10 dB). We compare this to a method for ideally preparing the same state by generating ${|{+}^{ϵ}⟩}_{\mathrm{GKP}}$ with $ϵ=0.1$, then applying a GKP phase gate $\hat{P}=e^{i\hslash {\hat{q}}^2/2}$. Since the phase gate requires inline squeezing (see Appendix pp8), we simulate measurement-based squeezing with various levels of ancillary squeezing and ancilla detector efficiency, which we label in (b)–(d). We see that phase gates apply a shearing effect to the envelope and peaks of the finite-energy GKP state, with lower quality measurement-based squeezing adding additional broadening to the peaks. This simulation is performed with the bosonic backend of Strawberry Fields.

Figure 13

The marginal distribution for a homodyne measurement along (a) $(q-p)/\sqrt{2}$ and (b) $(q+p)/\sqrt{2}$, each followed by a rescaling of the results by $\sqrt{2}$. Binning the result to the nearest $n\sqrt{\pi }$ and taking the parity of (a) $n$ or (b) $n+1$ effects a GKP qubit $Y$ measurement. We present the marginals for Wigner functions (a), (b), and (d) from Fig. 12. While marginal distributions along either $q-p$ or $q+p$ could be sampled for a Pauli $Y$ measurement, the narrower peaks along $q-p$ yield higher likelihood of falling within the correct bin, emphasizing the value of tracking shearing effects to optimize readout fidelity.

Figure 14

Teleportation of ${|0^{ϵ}⟩}_{\mathrm{GKP}}$ onto a $p$-squeezed state. Ideally, the teleported state should be ${|+⟩}_{\mathrm{GKP}}$, with Pauli $X$ ($Z$) measurement of the state yielding a qubit 0 outcome with 100% (50%) probability. However, with noisy teleportation circuits (loss and measurement-based squeezing in the cz gate), each instance of the circuit, conditioned on the probabilistic value of the teleportation feedforward, yields a slightly different state. We simulate the noisy teleportation circuit 500 times. In (a) and (b), we plot the mean and the (log of the) standard deviation for the probability of a Pauli $X$ measurement yielding 0, as a function of circuit loss and squeezing. Having found the mean probability of a Pauli $Z$ measurement yielding 0 to be roughly $50\mathrm{\%}$ across squeezing and loss levels, in (c) we provide the (log of the) standard deviation for the probability a Pauli $Z$ measurement yields 0.

Figure 15

(a) Wigner function for a GKP magic state ${|M^{ϵ}⟩}_{\mathrm{GKP}}=\hat{E}(ϵ)[(1/\sqrt{2})(e^{-i\pi /8}{|0⟩}_{\mathrm{GKP}}+e^{i\pi /8}{|1⟩}_{\mathrm{GKP}})]$ with $ϵ=0.1$ (10 dB). (b) Wigner function for t gate teleportation applied to ${|{+}^{ϵ}⟩}_{\mathrm{GKP}}$ with $ϵ=0.1$ (10 dB) using the circuit from Fig. 6 and postselecting on an outcome that does not require the feedforward phase gate. (c) Output from the same circuit, but postselecting on an outcome that does require the feedforward phase gate. The gate teleportation circuit for the t gate with realistic components is involved, but now simulatable using the bosonic backend of Strawberry Fields.

Figure 16

Mean probability of obtaining a qubit 0 outcome from different Pauli operator measurements after applying a t gate to ${|{+}^{ϵ}⟩}_{\mathrm{GKP}}$ with $ϵ=0.1$ (10 dB). We vary the per-peak squeezing of the GKP magic resource state and see how well the qubit readout matches the ideal readout values for the Pauli operators (dashed lines). We see that a per-peak squeezing of the order of 11 dB (close to that of the data state) is sufficient to retrieve correct readout to within 1%. The simulation provides another example of the usefulness of having two methods for performing Pauli $Y$ measurements, as seen in Fig. 13; homodyne detection and binning along $q-p$ yields a more faithful outcome.

Figure 17

Mean probability of obtaining a qubit 0 outcome from different Pauli operator measurements after applying a t gate to ${|{+}^{ϵ}⟩}_{\mathrm{GKP}}$ with $ϵ=0.1$ (10 dB). Here, the gate teleportation is effected using imperfect optical elements (see Figs. 6 and 18 for circuit breakdowns). We fix the ancillary magic state to have the same $ϵ$ as the data state; we employ 12 dB squeezed states for measurement-based squeezing, and set the homodyne efficiency to 99%. We see how readout changes as we vary $\eta$, which parameterizes the loss channel we apply after each of the five states are initialized and after each of the five beam splitters.

Figure 18

Review of CV gates and their realistic optical implementations. (a) The circuit for measurement-based squeezing [29]. An ancillary position eigenstate is combined at a beam splitter of angle $\theta$ with the target mode. The ancilla is then measured in $p$ quadrature with efficiency $\eta$, and informing a feedforward $q$ displacement of $p_0\tan \theta /\sqrt{\eta }$. The level of $q$-quadrature squeezing applied to the target mode is $r=\ln \cos \theta$. In practice, the quality of measurement-based squeezing is dependent on the level of squeezing of the ancillary state, and loss in the beam splitter. Squeezing in $p$ can be achieved with a momentum eigenstate, a $q$-quadrature measurement, and a $p$ displacement. Complex squeezing parameters can be implemented by sandwiching the circuit between two phase shifters. (b) A CV phase gate applied using rotations and inline squeezing [28]. Here, $r=\mathrm{arccosh}\sqrt{1+s^2/4}$, $\theta =\arctan (s/2)$, and $\varphi =-(\pi /2)\mathrm{sign}(s)-\theta$. (c) A CV $\mathrm{cz}$ gate [28], with $r=\mathrm{arcsinh}-s/2$, $\varphi =\pi /2$, $\sin 2\theta =(\cosh r{)}^{-1}$, and $\cos 2\theta =\tanh r$. For the GKP encoding, the qubit phase and $\mathrm{cz}$ gates correspond to $s=1$. In practice, the inline squeezing in the phase and $\mathrm{cz}$ gates can be implemented using measurement-based squeezing, and additional noise in the $\mathrm{cz}$ gate can be modeled by lossy beam splitters.

Figure 19

Simulating Gaussian measurement outcomes using rejection sampling.