Non-Gaussian Photonic State Engineering

Abstract

Generation of a non-Gaussian quantum state uses a quantum optical frequency comb source comprising at least one nonlinear optical medium and providing spectral modes including a plurality of pairs of entangled spectral modes, each pair including: a first spectral mode centered at a first frequency, and a second spectral mode entangled with the first spectral mode and centered at a second frequency that is spaced from the first frequency at a multiple of the frequency-bin spacing. Electrically controllable optical transformation modules are connected in series with a first module receiving the spectral modes. Two or more of the modules each comprise: a spectral mode phase shifter applying respective phase shifts to different spectral modes based on at least a first electrical signal, and a spectral mode mixer coupled to the spectral mode phase shifter and coupling spectral modes centered at different frequencies based on at least a second electrical signal.

Description

TECHNICAL FIELD

This disclosure relates to non-gaussian photonic state engineering.

BACKGROUND

Gaussian states of bosonic modes can be described as quantum states of optical waves that can be prepared using quadrature squeezed optical waves and passive linear optical elements. Gaussian states and their feasibility of experimental production are useful for quantum information processing. However, Gaussian states and Gaussian measurements (e.g., homodyne and heterodyne detection) do not constitute a universal set, i.e., resources that would allow universal quantum computation. Thus, techniques for generating non-Gaussian states are useful.

SUMMARY

In one aspect, in general, an apparatus for generating a non-Gaussian quantum state associated with one or more spectral modes comprises: a quantum optical frequency comb source comprising at least one nonlinear optical medium and configured to provide a plurality of spectral modes spaced at a frequency-bin spacing, where the spectral modes include a plurality of pairs of entangled spectral modes, each pair of entangled spectral modes including: a first spectral mode centered at a first frequency, and a second spectral mode entangled with the first spectral mode and centered at a second frequency that is spaced from the first frequency at a multiple of the frequency-bin spacing; and a plurality of electrically controllable optical transformation modules connected in series with a first module in the series receiving spectral modes from the quantum optical frequency comb, where two or more of the modules each comprise: a spectral mode phase shifter configured to apply respective phase shifts to different spectral modes based on at least a first electrical signal, and a spectral mode mixer coupled to the spectral mode phase shifter and configured to couple spectral modes centered at different frequencies based on at least a second electrical signal.

Aspects can include one or more of the following features.

One or more of the electrically controllable optical transformation modules are configured to transform each pair of entangled spectral modes into single-mode squeezed vacuum states.

Each spectral mode phase shifter comprises a Fourier-transform pulse shaper.

Each spectral mode mixer comprises an electro-optic phase modulator.

The apparatus further comprises an output interface coupled to a last module in the series, the output interface including: a plurality of photon number resolving (PNR) detectors configured to detect respective spectral modes of more than one and fewer than all spectral modes output from the last module in the series, and at least one port providing at least one output spectral mode from the last module in the series not detected by any of the PNR detectors.

The output interface further comprises a wavelength-dependent element between the last module and the plurality of PNR detectors.

The wavelength-dependent element comprises at least one of: a grating or a prism.

Each PNR detector is configured to generate a detection signal that distinguishes a detected photon number equal to zero from a detected photon number equal to one, in each of a plurality of time slots; and distinguishes a detected photon number equal to one from at least one detected photon number greater than one, in each of the plurality of time slots.

The apparatus further comprises a trigger module configured to generate a trigger signal based on detection signals from the plurality of PNR detectors indicating one or more time slots in which the port providing at least one output spectral mode corresponds to a non-Gaussian quantum state.

The non-Gaussian quantum state comprises a quantum superposition of continuous variable (CV) optical wave quadrature states.

The quantum superposition comprises a superposition of optical waves with opposite phases.

The non-Gaussian quantum state is input to a module that includes one or more optical elements that provide at least one Gottesman-Kitaev-Preskill (GKP) qubit.

The one or more optical elements include at least one optical beamsplitter.

The trigger signal is based on detection signals that include a plurality of the detection signals each indicating a photon number of zero in at least one time slot.

The last module in the series consists essentially of a spectral mode mixer configured to couple spectral modes centered at different frequencies based on an electrical signal.

The first module in the series consists essentially of a spectral mode mixer configured to couple spectral modes centered at different frequencies based on an electrical signal.

A plurality of the modules in the series are integrated on a common photonic integrated circuit.

The frequency-bin spacing is a free spectral range of an optical parametric oscillator that includes the nonlinear optical medium.

In another aspect, in general, a method for generating a non-Gaussian quantum state associated with one or more spectral modes comprises: providing from a quantum optical frequency comb source, comprising at least one nonlinear optical medium, a plurality of spectral modes spaced at a frequency-bin spacing, where the spectral modes include a plurality of pairs of entangled spectral modes, each pair of entangled spectral modes including: a first spectral mode centered at a first frequency, and a second spectral mode entangled with the first spectral mode and centered at a second frequency that is spaced from the first frequency at a multiple of the frequency-bin spacing; and receiving spectral modes from the quantum optical frequency comb into a first module of a plurality of electrically controllable optical transformation modules connected in series, where two or more of the modules each comprise: a spectral mode phase shifter configured to apply respective phase shifts to different spectral modes based on at least a first electrical signal, and a spectral mode mixer coupled to the spectral mode phase shifter and configured to couple spectral modes centered at different frequencies based on at least a second electrical signal.

In another aspect, in general, an apparatus for generating a non-Gaussian quantum state associated with one or more spectral modes comprises: a quantum optical frequency comb source comprising at least one nonlinear optical medium and configured to provide a plurality of spectral modes spaced at a frequency-bin spacing, a plurality of electrically controllable optical transformation modules connected in series with a first module in the series receiving spectral modes from the quantum optical frequency comb, where two or more of the modules each comprise: a spectral mode phase shifter configured to apply respective phase shifts to different spectral modes based on at least a first electrical signal, and a spectral mode mixer coupled to the spectral mode phase shifter and configured to couple spectral modes centered at different frequencies based on at least a second electrical signal; and an output interface coupled to a last module in the series, the output interface including: a plurality of photon number resolving (PNR) detectors configured to detect respective spectral modes of more than one and fewer than all spectral modes output from the last module in the series, and at least one port providing at least one output spectral mode from the last module in the series not detected by any of the PNR detectors; wherein each PNR detector is configured to generate a detection signal that: distinguishes a detected photon number equal to zero from a detected photon number equal to one, in each of a plurality of time slots, and distinguishes a detected photon number equal to one from at least one detected photon number greater than one, in each of the plurality of time slots.

Aspects can include one or more of the following features.

The number of spectral modes spaced at a frequency-bin spacing is three or more.

The apparatus further comprises a trigger module configured to generate a trigger signal based on detection signals from the plurality of PNR detectors indicating one or more time slots in which the port providing at least one output spectral mode corresponds to a non-Gaussian quantum state.

The plurality of spectral modes provided by the quantum optical frequency comb source comprise a Gaussian state.

The plurality of spectral modes provided by the quantum optical frequency comb source comprise a state continuous-variable encoded state.

Aspects can have one or more of the following advantages.

The techniques described herein enable various implementations of non-Gaussian state production from input states populating discrete frequency bins. Some of the techniques include the use of controllable unitary operations with a quantum frequency processor, followed by photon-number-resolved (PNR) detection of ancilla modes. A quantum state of a single optical mode can be engineered at one reference frequency in a scalable way by leveraging Gaussian boson sampling (GBS) type state preparation (e.g., partial PNR detection on a multi-mode Gaussian state) across multiple frequencies using a comb-based continuous variables entanglement source and a quantum frequency processor to implement a linear mixing across the frequencies, and a heralding trigger obtained from detecting a subset of spectral modes.

Other features and advantages will become apparent from the following description, and from the figures and claims.

BRIEF DESCRIPTION OF THE DRAWINGS

The disclosure is best understood from the following detailed description when read in conjunction with the accompanying drawing. It is emphasized that, according to common practice, the various features of the drawing are not to-scale. On the contrary, the dimensions of the various features are arbitrarily expanded or reduced for clarity.

FIG. 1 is a schematic diagram of an example non-Gaussian quantum state generation system.

FIG. 2 is a schematic diagram of an example generation scheme for heralded non-Gaussian states.

FIG. 3A is a schematic diagram of an example non-Gaussian quantum state generation system.

FIG. 3B is a schematic diagram of an example non-Gaussian quantum state generation system.

FIG. 4A is a set of plots of wavefunctions in the quadrature basis (top panels) and photon number probabilities (bottom panels) for the generation of Schrödinger cat states.

FIG. 4B is a set of plots of wavefunctions in the quadrature basis (top panels) and photon number probabilities (bottom panels) for the generation of Schrödinger cat states.

FIG. 4C is a set of plots of wavefunctions in the quadrature basis (top panels) and photon number probabilities (bottom panels) for the generation of Schrödinger cat states.

FIG. 5 is a set of plots of the fidelity (top panels) and success probability (bottom panels) for the generation of Schrödinger cat states.

DETAILED DESCRIPTION

Non-Gaussian quantum states of light are useful resources for optical quantum information processing, but methods to generate them efficiently can be challenging to implement. Here we describe a generic approach for non-Gaussian state production from input states populating discrete frequency bins by using controllable unitary operations with a quantum frequency processor, followed by photon-number-resolved detection of ancilla modes. Leveraging and refining the K-function representation of quantum states in the coherent basis, we develop a theoretical model amenable to numerical optimization and, as specific examples, design quantum frequency processor circuits for the production of Schrödinger cat states, exploring the performance tradeoffs for several combinations of ancilla modes and circuit depth. Our scheme provides a valuable general framework for producing complex quantum states in frequency bins, paving the way for single-spatial-mode, fiber-optic-compatible non-Gaussian resource states.

FIG. 1 shows an example non-Gaussian quantum state generation system 100 for generating a non-Gaussian quantum state associated with one or more spectral modes. The system 100 includes a quantum optical frequency comb source 102 comprising at least one nonlinear optical medium and configured to provide a plurality spectral modes spaced at a frequency-bin spacing (e.g., a free spectral range of an optical parametric oscillator that includes the nonlinear optical medium). In some implementations, the spectral modes include multiple pairs of entangled spectral modes. In such implementations, each pair of entangled spectral modes includes: a first spectral mode centered at a first frequency, and a second spectral mode entangled with the first spectral mode and centered at a second frequency that is spaced from the first frequency at a multiple of the frequency-bin spacing. The system 100 includes a module 104 called a “quantum frequency processor,” which includes multiple electrically controllable optical transformation modules connected in series. The module 104 includes a first module 106A in the series receiving spectral modes from the quantum optical frequency comb source 102, and a last module 106Z in the series coupled to an output interface 108. The modules in the series are not necessarily all identical. For example, the first and/or last module may include a subset of the sub-components of the other modules. In this example, two or more of the modules each include: a spectral mode phase shifter configured to apply respective phase shifts to different spectral modes based on at least a first electrical signal, and a spectral mode mixer coupled to the spectral mode phase shifter and configured to couple spectral modes centered at different frequencies based on at least a second electrical signal. The output interface 108 may include, or couple to, a wavelength-dependent element. The wavelength dependent element could be dispersive (e.g., a waveguide or a prism) or diffractive (e.g., a grating). The wavelength-dependent element can be used to spatially separate the spectral modes based upon their frequency, allowing for PNR detection of each spectral mode.

The distinction between discrete-variable (DV) and continuous-variable (CV) encodings offers a valuable lens through which to classify and understand photonic quantum information processing systems. Quantum information encoded on quantum states of optical waves can use DV encoding where a physical observable has one of multiple possible discrete values, or CV encoding where a physical observable has a value in a continuous range. Based on true (or approximate) finite-dimensional Hilbert spaces, DV optical designs are typically associated with qubits (or qudits) encoded in photons that are manipulated and subsequently measured with single-photon detectors. On the other hand, the infinite-dimensional Hilbert spaces of CV quantum information exploit collective photonic excitations (such as coherent or squeezed states) and homodyne/heterodyne detection with local oscillators (CV maintaining measurement techniques) and/or photon number resolving detectors (a measurement that can be used to convert CV into DV encoding), as fundamental resources. Gaussian states are an example of CV encoding and can be simpler to generate, compared to some realizations of DV encoding, which are non-Gaussian. The non-Gaussian quantum state generation techniques described herein can transform a Gaussian state input into a DV encoded output state. Other techniques have assumed that DV encoded states may be readily available as input states, which due to a variety of factors associated with practical implementations, is not always true. Thus, the ability to generate DV encoded states from CV encoded Gaussian states can be advantageous. From a technical side, the DV/CV divide can prove quite stark, and significant differences appear theoretically as well: for example, security proofs for CV quantum key distribution have generally proven much more challenging to establish due to the infinite dimensionality involved.

Yet this dichotomy is far from absolute, with features of particular quantum information processing approaches blurring the CV/DV distinction entirely. At the implementation level, many DV photonic systems utilize subspaces taken from a larger, intrinsically continuous Hilbert space-time and frequency bins forming representative examples of relevance. For encodings such as the Gottesman-Kitaev-Preskill (GKP) qubit, the logical quantum information is discrete, but the encoding occupies the full continuous Hilbert space. Here the CV aspects are not incidental features of the chosen Hilbert space; rather, they prove critical to the paradigm itself, providing the foundation for measuring and correcting continuous errors on the logical qubit state.

The potential of error-corrected photonic quantum information processing with GKP qubits makes them an appealing direction for research. But producing such states—and non-Gaussian CV states more generally—can be challenging. In some of the examples described herein, the frequency-bin degree of freedom (DoF) can exhibit several attractive features for scalable photonic quantum information processing, including wavelength parallelizability, compatibility with single-mode optical fiber, and CV state production with resonant parametric oscillators, both free-space and integrated.

A potential challenge of non-Gaussian state production with frequency-bin encoding, however, is the realization of arbitrary unitary operations. A quantum frequency processor (QFP) based on alternating applications of electro-optic phase modulators (EOMs) and pulse shapers can in principle synthesize any unitary frequency-bin operation in a scalable fashion.

Some of the prophetic examples described herein are based on simulation of a model for non-Gaussian frequency-bin state engineering on the QFP. We describe a resource-efficient method for computing the output of a QFP excited by Gaussian inputs and measured with photon-number-resolving (PNR) detectors applied to a subset of frequency modes. As examples of this general approach, we describe implementations of QFP circuits intended to produce Schrödinger cat states in one undetected bin and explore the impact of the number of components and ancilla modes on circuit performance, according to a cost function that balances both state fidelity and success probability. Some of the example implementations described herein can provide a general framework for non-Gaussian state production in frequency-bin quantum systems.

Without intending to be bound by theory, as an example of modeling such a system, we use a representation of Gaussian states in the coherent basis according to the K-function formalism. It has been shown that any N-mode pure Gaussian state |Ψ with covariance matrix (CM) V and displacement vector {right arrow over (x)}β can be written in the coherent basis |{right arrow over (α)} as

$\begin{array}{cc}❘\Psi \rangle =\int d^{2N}{\overrightarrow{x}}_{\alpha }K\left({\overrightarrow{x}}_{\alpha }\right)❘\overrightarrow{\alpha }\rangle ,& \left(1\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}K\left({\overrightarrow{x}}_{\alpha }\right)=\frac{e^{-\frac{1}{2}{\left({\overrightarrow{x}}_{\alpha }-{\overrightarrow{x}}_{\beta }\right)}^Tℬ\left({\overrightarrow{x}}_{\alpha }-{\overrightarrow{x}}_{\beta }\right)+\frac{1}{2}{\overrightarrow{x}}_{\alpha }^T\gamma {\overrightarrow{x}}_{\beta }}}{{\left(2\pi \right)}^N{\left(\det \Gamma \right)}^{1/4}},& \left(2\right)\end{array}$ $\mathrm{with}\Gamma =V+I/2,$ $\begin{array}{cc}ℬ=\frac{1}{2}\left(\begin{array}{cc}A+\frac{i}{2}\left(C+C^T\right)& C-\frac{i}{2}\left(A-B\right)\\ C^T-\frac{i}{2}\left(A-B\right)& B-\frac{i}{2}\left(C+C^T\right)\end{array}\right),& \left(3\right)\end{array}$ $\begin{array}{cc}\gamma =\left(\begin{array}{cc}0& iI\\ -\mathrm{iI}& 0\end{array}\right),& \left(4\right)\end{array}$ $\mathrm{where}A=A^T,B=B^T,\mathrm{and}C\mathrm{are}\mathrm{defined}\mathrm{as}\mathrm{the}\mathrm{blocks}\mathrm{of}{\Gamma }^{-1}\mathrm{as}\mathrm{follows}:$ $\begin{array}{cc}{\Gamma }^{-1}=\left(\begin{array}{cc}A& C\\ C^T& B\end{array}\right).& \left(5\right)\end{array}$

We note that since the CM V is symmetric, Γ and Γ⁻¹ are also symmetric. We work with the convention =1 (therefore the CM of vacuum is I/2) and consider the qqpp representation where vectors are defined as {right arrow over (x)}_α^T=({right arrow over (q)}_α^T,{right arrow over (p)}_α^T) with {right arrow over (q)}_α^T=(q_α1, . . . , q_αN) and {right arrow over (p)}_α^T=(p_α1, . . . , p_αN) the canonical position and momentum vectors. The volume element for integration is then defined as d^2N{right arrow over (x)}_α=dq_α1 . . . dq_αNdp_α1 . . . dp_αN, and α_i=(q_αi+ip_αi)/√{square root over (2)}.

The coherent basis representation is a valuable tool for working on photon-subtraction-based or, more generally, partial PNR detection schemes aimed at engineering Gaussian states into desired non-Gaussian states. Photon subtraction can be modelled either (i) as a beamsplitter whose two input ports are fed with the ith mode of |Ψ and vacuum |0, respectively, followed by PNR detection on the lower output port; or (ii) simply by acting the annihilation operator {circumflex over (α)}_i, where the index i refers to the mode, on |Ψ. Therefore, the photon subtraction operator will act only on the basis vectors of the state, i.e., coherent states in this instance. The action of beamsplitters or annihilation operators on coherent states is straightforward, making this basis particularly efficient for analytical or numerical evaluation. The situation is similar for partial PNR detection on a Gaussian state written as a coherent state expansion; the projection of a coherent state on a Fock state is the well known expression

$\langle n❘\alpha \rangle =\exp \left(-{❘\alpha ❘}^2/2\right){\alpha }^n/\sqrt{n!}.$

The probability of a length-N PNR pattern for an N-mode Gaussian state |Ψ with zero displacements, i.e., {right arrow over (x)}_β=0 in Eq. (1), is given by

$\begin{array}{cc}\begin{array}{c}{P}_{n_1\dots n_N}={❘\langle n_1\dots n_N❘\Psi \rangle ❘}^2\\ =\frac{1}{\det \mathscr{H}\sqrt{\det \Gamma }\stackrel{N}{\prod _{i=1}}{n}_i!2^{n_i}}{❘I_{n_1\dots n_N}❘}^2,\end{array}& \left(6\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{I}_{n_1\dots n_N}=\int d^{2N}{\overrightarrow{x}}_{\alpha }R\left({\overrightarrow{x}}_{\alpha }\right)\prod _{i=1}^N{\left(q_{{\alpha }_i}+i{p}_{{\alpha }_i}\right)}^{n_i},& \left(7\right)\end{array}$ $\begin{array}{cc}R\left({\overrightarrow{x}}_{\alpha }\right)=\frac{\sqrt{\det \mathscr{H}}}{{\left(2\pi \right)}^N}{e}^{-\frac{1}{2}{\overrightarrow{x}}_{\alpha }^T\mathscr{H}{\overrightarrow{x}}_{\alpha }},& \left(8\right)\end{array}$ $\mathrm{and}\mathscr{H}=ℬ+I/2.\mathrm{Equation}\left(7\right)\mathrm{can}\mathrm{be}\mathrm{rewritten}\mathrm{as}$ $\begin{array}{cc}{I}_{n_1\dots n_N}=\{\begin{array}{cc}0& \Sigma =\mathrm{odd},\\ \mathrm{Hf}\left(\sigma \right)& \Sigma =\mathrm{even},\end{array}& \left(9\right)\end{array}$

where Σ=Σ_i=1^Nn_i, Hf (σ) is the loop hafnian of the matrix σ with elements σ_ij=s_is_j, where 1≤i, j≤Σ and s_i=q_αi+ip_αi. The hafnian in Eq. (9) represents the mean value s₁ⁿ¹ . . . s_N^nN under the Gaussian distribution of Eq. (8).

We will derive the explicit relation of the matrix σ to the matrix and consequently to matrices Γ and V. We also give the expressions for the Fock expansion coefficients of the produced non-Gaussian states and simplify further the expressions.

The matrix Γ is defined as Γ=V+I/2, where V is the CM and I the identity matrix. Since V corresponds to a pure Gaussian state, it can be written as V=S_pV₀S_p^T where S_p is an orthogonal symplectic matrix for a general passive transformation (beamsplitters and phase rotations, but not squeezers) and V₀ is the CM for a product of N single mode squeezed vacuum states, i.e., the diagonal matrix

$\begin{array}{cc}{V}_0=\frac{1}{2}\mathrm{diag}\left(e^{2r_1},\dots ,e^{2r_N},e^{-2r_1},\dots ,e^{-2r_N}\right),& \left(10\right)\end{array}$

where r₁, . . . , r_N are the real and positive squeezing parameters for each of the N single-mode squeezed vacuum states (note that the phase of the squeezing has been absorbed into the orthogonal symplectic transformation S_p).

We have the following relation,

$\begin{array}{cc}\det \Gamma =\det \left[S_p\left(V_0+\frac{I}{2}\right)S_p^T\right]& \left(11\right)\end{array}$ $\begin{array}{cc}=\det S_p\det \left(V_0+\frac{I}{2}\right)\det S_p^T,& \left(12\right)\end{array}$ $\mathrm{from}\mathrm{which}\mathrm{we}\mathrm{write}$ $\begin{array}{cc}\det \Gamma =\det \left(V_0+\frac{I}{2}\right)& \left(13\right)\end{array}$

since detS_p=detS_p^T=1 as both S_p and S_p^T, are symplectic matrices. The right hand side of Eq. (13) is the determinant of a diagonal matrix from which we find

$\begin{array}{cc}\det \Gamma =\prod _{i=1}^N{\cosh }^2{r}_i.& \left(14\right)\end{array}$

Therefore, Eq. (6) is rewritten as

$\begin{array}{cc}{P}_{n_1\dots n_N}=\frac{{❘I_{n_1\dots n_N}❘}^2}{\det \mathscr{H}\stackrel{N}{\prod _{i=1}}{n}_i!2^{n_i}\cosh r_i}.& \left(15\right)\end{array}$

In the case where the input squeezing is the same among all single mode squeezed vacuum states, i.e. r₁= . . . =r_N=r, Eq. (14) reduces to detΓ=cosh^2Nr.

In order to simplify Eq. (5) we can write Γ=S_p (V₀+I/2) S_p^T, and since S_p^T−1=S_p is a symplectic orthogonal matrix we have

$\begin{array}{cc}{\Gamma }^{-1}={S_p\left(V_0+\frac{1}{2}\right)}^{-1}{S}_p^T.& \left(16\right)\end{array}$

The symplectic orthogonal matrix S_p has the following block matrix structure and properties:

$\begin{array}{cc}{S}_p=\left(\begin{array}{cc}{S}_A& S_B\\ -S_B& S_A\end{array}\right)& \left(17\right)\end{array}$ $\begin{array}{cc}{S}_A^T{S}_B=S_B^T{S}_A,& \left(18\right)\end{array}$ $\begin{array}{cc}{S}_A{S}_B^T=S_B{S}_A^T,& \left(19\right)\end{array}$ $\begin{array}{cc}{S}_A^T{S}_A+S_B^T{S}_B=I,& \left(20\right)\end{array}$ $\begin{array}{cc}{S}_A{S}_A^T+S_B{S}_B^T=I.& \left(21\right)\end{array}$

Moreover, since V₀ is diagonal we can write

$\begin{array}{cc}{\left(V_0+\frac{1}{2}\right)}^{-1}=I+\left(\begin{array}{cc}-T& 0\\ 0& T\end{array}\right),& \left(22\right)\end{array}$

where T=diag(tanh r₁, . . . , tanh r_N). In virtue of Eqs. (16), (17), and (19), we find that in Eq. (5)

$\begin{array}{cc}A=-S_{A}T{S}_A^T+S_{B}T{S}_B^T,& \left(23\right)\end{array}$ $\begin{array}{cc}C=C^T=S_{A}T{S}_B^T+S_{B}T{S}_A^T,& \left(24\right)\end{array}$ $\begin{array}{cc}A+B=2I.& \left(25\right)\end{array}$

Therefore, in the most general case possible, Eq. (5) is simplified to

$\begin{array}{cc}{\Gamma }^{-1}=\left(\begin{array}{cc}A& C\\ C& 2I-A,\end{array}\right)& \left(26\right)\end{array}$

where A and C are given in Eqs. (23) and (24), respectively, as functions of the passive symplectic transformation S_p and the input squeezing parameters.

Consequently, matrix of Eq. (3) simplifies to

$\begin{array}{cc}ℬ=\frac{1}{2}\left(\begin{array}{cc}A+\mathrm{iC}& C-i\left(A-I\right)\\ C-i\left(A-I\right)& 2I-A-\mathrm{iC}\end{array}\right).& \left(27\right)\end{array}$

The matrix appearing in Eq. (8) is defined as

$\begin{array}{cc}\mathscr{H}=ℬ+I/2.& \left(28\right)\end{array}$

We find it easier if we transform as =W^†W using the unitary matrix W defined as

$\begin{array}{cc}W=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}I& I\\ -\mathrm{iI}& iI\end{array}\right).& \left(29\right)\end{array}$

Utilizing Eqs. (27), (28), and (29) we find

$\begin{array}{cc}\tilde{\mathscr{H}}=\left(\begin{array}{cc}I& A-I+\mathrm{iC}\\ 0& I\end{array}\right),& \left(30\right)\end{array}$

from which we see that det =det I=1. Since |det W|²=1, we have det =det and conclude that

$\begin{array}{cc}\det \mathscr{H}=1.& \left(31\right)\end{array}$

Therefore, Eqs. (8) and (15) are further simplified to

$\begin{array}{cc}{P}_{n_1\dots n_N}=\frac{{❘I_{n_1\dots n_N}❘}^2}{\stackrel{N}{\prod _{i=1}}{n}_i!2^{n_i}\cosh r_i},& \left(32\right)\end{array}$ $\begin{array}{cc}R\left({\overrightarrow{x}}_{\alpha }\right)=\frac{1}{{\left(2\pi \right)}^N}{e}^{-\frac{1}{2}{\overrightarrow{x}}_{\alpha }^T\mathscr{H}{\overrightarrow{x}}_{\alpha }}.& \left(33\right)\end{array}$

In order to derive a convenient expression for , we work with and observe that

$\begin{array}{cc}{\tilde{\mathscr{H}}}^{-1}=\left(\begin{array}{cc}I& -\left(A-I+\mathrm{iC}\right)\\ 0& I\end{array}\right)& \left(34\right)\end{array}$

is indeed the inverse of , i.e., it satisfies =1. Since =W^†W we find that =WW^\ and finally

$\begin{array}{cc}{\mathscr{H}}^{-1}=\frac{1}{2}\left(\begin{array}{cc}3I-A-iC& i\left(A-I+iC\right)\\ i\left(A-I+iC\right)& I+A+iC\end{array}\right).& \left(35\right)\end{array}$

Therefore, using Eqs. (23), (24), and (35), any given passive symplectic transformation S_p, and input squeezing parameters, one can readily write .

Making use of Eq. (33), we can express the matrix elements of σ as

$\begin{array}{cc}{\sigma }_{\mathrm{ij}}=\langle \left(q_{{\alpha }_i}+{\mathrm{ip}}_{{\alpha }_i}\right)\left(q_{{\alpha }_j}+i{p}_{{\alpha }_j}\right)\rangle =\frac{1}{{\left(2\pi \right)}^N}\int d^{2N}{\overrightarrow{x}}_{\alpha }\exp \left(-\frac{1}{2}{\overrightarrow{x}}_{\alpha }^T\mathscr{H}{\overrightarrow{x}}_{\alpha }\right)\times \left(q_{{\alpha }_i}+i{p}_{{\alpha }_i}\right)\left(q_{{\alpha }_j}+i{p}_{{\alpha }_j}\right)=\frac{d}{d{\lambda }_i}\frac{d}{d{\lambda }_j}\exp \left(\frac{1}{2}{\overrightarrow{\Lambda }}^T{\mathscr{H}}^{-1}\overrightarrow{\Lambda }\right){❘}_{\overrightarrow{\Lambda }=\overrightarrow{0}},& \left(36\right)\end{array}$

where {right arrow over (Λ)}^T=({right arrow over (λ)}^T, i{right arrow over (λ)}^T) is a 2N-dimensional vector with {right arrow over (λ)}^T=(λ₁, . . . , λ_N) a real N-dimensional vector. Viewing

$\frac{1}{2}{\overrightarrow{\Lambda }}^T{\mathscr{H}}^{-1}\overrightarrow{\Lambda }$

in the exponential of the right hand side of Eq. (36) as a polynomial in λ_i, Eq. (36) is equal to the coefficient of λ_iλ_j. This way, we can write

$\begin{array}{cc}{\sigma }_{ij}=2\left({\mathscr{H}}_{\mathrm{ij}}^{-1}-{\mathscr{H}}_{i+Nj+N}^{-1}\right).& \left(37\right)\end{array}$

From the covariance matrix V, one can find matrix Γ⁻¹ and therefore matrix σ using Eqs. (35) and (37), which is used in the calculation in Eq. (9).

The Gaussian moment problem of Eq. (7) represents a hafnian calculation and is related to the Gaussian boson sampling paradigm. When the indices i, j are equal this corresponds to a loop, i.e., matching an object with itself. Therefore, it is typically referred to as a loop hafnian. Equation (32) is the probability of finding n_i photons in each one of the i=1, . . . , N modes. If we wish to engineer the N-mode Gaussian state into an M-mode (M<N) non-Gaussian, one as shown in FIG. 2, we leave M modes undetected; without loss of generality we assume the undetected modes are the M upper modes. The probability of the PNR pattern (n_M+1, . . . , n_N) on the lower detected modes is precisely the probability P_{nM+1, . . . , nN} of producing the corresponding non-Gaussian state. This probability is

$\begin{array}{cc}P\equiv P_{n_{M+1},\dots ,n_N}=\stackrel{\infty }{\sum _{n_1,\dots ,n_M=0}}{P}_{n_1,\dots ,n_N}.& \left(38\right)\end{array}$

For numerical simulations, the above sum is truncated to a finite upper limit, which should be chosen with care to ensure that it encompasses all Fock coefficients of nonnegligible probability. This condition can be verified in practice by successively increasing the limits and observing no change to P.

FIG. 2 shows an example system 200 configured for heralding an M-mode state |Φ from N single-mode, zero-displacement squeezed resource states and N×N unitary operation U. Partial PNR detection on the N−M lower modes by photodetectors 202, each configured as a PNR detector to distinguish among different numbers of detected photons, produces a non-Gaussian state on the undetected M modes. For example, a PNR detector can be implemented using detectors incorporating cryogenic sensors, and/or superconducting materials. One example of a PNR detector is a transition-edge sensor (TES). Another example of a PNR detector is a segmented detector based on single-photon avalanche-photodiodes.

In some implementations of a PNR detector, the PNR detector can be configured to be able to distinguish a detected photon number equal to zero from a detected photon number equal to one, in each of a plurality of time slots, and to also be able to distinguish a detected photon number equal to N from a detected photon number of N+1, in each of the plurality of time slots, for each value of N up to at least some relatively large number (e.g., N at least 10 or greater). In some cases, it may not be necessary for the PNR detector to be able to distinguish any possible number of detected photons from any other possible number of detected photons. For example, for a particular system configuration, it may only be necessary to distinguish among photon numbers that are more likely to be detected.

The non-Gaussian state |Φ on the M undetected modes can be written as a partial projection on Fock states of the detected modes:

$\begin{array}{cc}❘\Phi \rangle =\frac{1}{\sqrt{P}}\langle n_{M+1}\dots n_N❘\Psi \rangle ,& \left(39\right)\end{array}$

where P is given in Eq. (38) and |Ψ is the input N-mode Gaussian state.

The Fock expansion coefficients of heralded state |Φ are c_{n1. . . nM}=n₁ . . . n_M|Φ. Using Eqs. (7) and (39) we find

$\begin{array}{cc}{c}_{n_1\dots n_M}=\frac{I_{n_1\dots n_M{n}_{M+1}\dots n_N}}{\sqrt{P}\prod _{i=1}^N\sqrt{n_i!2^{n_i}\cos h{r}_i}},& \left(40\right)\end{array}$

where the numerator is given by Eq. (7). Therefore, for any given partial PNR pattern (n_M+1, . . . , n_N) one can compute the Fock expansion coefficients of the produced state |Φ, which can be benchmarked against a target non-Gaussian state |Φ_t through direct comparison of Fock coefficients or collectively through fidelity =|Φ_t|Φ|².

Regarding the example model framework thus far described, we note the following. First, we note that our formalism provides an approach to computing Gaussian states in the Fock basis. By incorporating the reduced dimensionality of a pure state directly, our approach may call for calculation of fewer expansion coefficients to fully characterize the output, in the case of pure state evolution.

Second, extra care may be needed when dealing with loop hafnians. Let us give an example. Say that one wants to calculate (s₁²s₂s₃). To apply Wick's formula for implementing a practical application of the modeled system, one may rewrite the mean value as containing four different objects, i.e., s₁²s₂s₃=g₁g₂g₃g₄. Wick's formula gives the computed matchings as g₁g₂g₃g₄=g₁g₂g₃g₄+g₁g₃g₂g₄+g₁g₄g₂g₃, and then we substitute back g₁=g₂=s₁, g₃=s₂, and g₄=s₃, which gives

$\langle s_1^2{s}_2{s}_3\rangle =\langle s_1^2\rangle \langle s_2{s}_3\rangle +\langle s_1{s}_2\rangle \langle s_1{s}_3\rangle +\langle s_1{s}_3\rangle \langle s_1{s}_2\rangle .$

For the calculation of loop hafnians, we give a formula which is the nonzero result of Eq. (9):

$\begin{array}{cc}\mathrm{Hf}\left(\sigma \right)=\frac{1}{\left(\frac{\sum }{2}\right)!}\sum _{v_1=0}^{n_1}\dots \sum _{v_N=0}^{n_N}{\left(-1\right)}^{v_1+\dots +v_N}\times \left(\begin{array}{c}{n}_1\\ v_1\end{array}\right)\dots \left(\begin{array}{c}{n}_N\\ v_N\end{array}\right){\left(\frac{1}{2}{\overrightarrow{h}}^T\sigma \overrightarrow{h}\right)}^{\frac{\Sigma }{2}}& \left(41\right)\end{array}$

where {right arrow over (h)}^T=(n₁/2−v₁, . . . , n_N/2−v_N). Equations (7) and (41) can be used directly in Eq. (38) for the probability of finding any non-Gaussian state in the undetected modes and in Eq. (40) for the Fock expansion coefficients of such a state. This tailored expression for the loop hafnian may provide significant computational speed up for non-Gaussian state engineering work. Essentially, the improvement can be obtained when the dominant bottleneck in Wick's formula stems from repeated factors (e.g., s₁ in the example above) rather than many non-repeated factors, (e.g., s₂ and s₃ in the example above).

Third and finally, the formulas above enable calculation of the coefficients (n₁ . . . n_N|Ψ for any diagonal input covariance matrix V₀ and passive symplectic mode transformation S_p—i.e., any covariance matrix for a pure Gaussian state—without numerical evaluation of a single matrix inverse or determinant: these expressions have all been reduced to matrix or scalar operations in the above examples. This simplification has an impact on the efficiency of the numerical procedure, eliminating potentially time-consuming inverse calculations from the optimization loop.

Up to this point, the mathematical formulation has been general with respect to the underlying optical modes, applicable equally well to any photonic DoF. We now refine our focus to frequency bins specifically. Fundamentally, the QFP is designed to realize arbitrary unitary operations on a discrete set of equispaced, clearly separated frequency modes, or bins.

In order to understand the basic principles of operation, consider a discrete set of frequency modes, each centered at ω_n=ω₀+nΔω(n∈) and associated with an annihilation operator {circumflex over (α)}_n. The corresponding output operators {circumflex over (b)}_n relate to the inputs {circumflex over (α)}_n as

$\begin{array}{cc}{\hat{b}}_n=e^{i{\varphi }_n}{\hat{a}}_n& \left(42\right)\end{array}$

for a line-by-line pulse shaper and

$\begin{array}{cc}{\hat{b}}_n=\sum _{k=-\infty }^{\infty }{f}_{n-k}{\hat{a}}_k& \left(43\right)\end{array}$

for an EOM driven with phase function φ(t) periodic at the inverse mode spacing

$T=\frac{2\pi }{\mathrm{\Delta \omega }},$

so that

$e^{i\phi \left(t\right)}={\sum }_n{f}_n{e}^{-\mathrm{in}\Delta \omega t}\mathrm{and}{f}_n=\frac{1}{T}{\int }_0^T\mathrm{dt}{e}^{i\phi \left(t\right)}{e}^{\mathrm{in}\Delta \omega t}.$

As written, this formulation contains an infinite number of frequency bins; in the interests of numerical tractability, though, we can limit the total number of considered modes to N and discretize the temporal period as

$t_n=\frac{nT}{N}$

(n∈{0, 1, . . . , N−1}). Under this approximation, the total N×N unitary for a sequence of Q components becomes

$\begin{array}{cc}U=\left(F{D}_Q{F}^{†}\right)D_{Q-1}\dots \left({\mathrm{FD}}_3{F}^{†}\right)D_2\left(F{D}_1{F}^{†}\right),& \left(44\right)\end{array}$

where F is the discrete Fourier transform with elements

$F_{mn}=\frac{1}{\sqrt{N}}{e}^{2\pi \mathrm{imn}/N}$

(m,n∈{0, 1, . . . , N−1}). Each D_q is a diagonal unitary matrix; the odd-numbered q signify an EOM with elements

${\left(D_q\right)}_{nn}=e^{i{\phi }^{\left(q\right)}\left(t_n\right)},$

and the even-numbered q a pulse shaper with

${\left(D_q\right)}_{nn}=e^{i{\varphi }_n^{\left(q\right)}}.$

We bookend the QFP with EOMs in our example, rather than pulse shapers, based on previous experience where we have observed no increase in circuit performance with the addition of a front- or back-end pulse shaper. The form in Eq. (44) accurately reflects the physical situation as long as N is sufficiently large so that photon probability amplitudes do not reach the edge of the truncated simulated domain and artificially “wrap around” to the other side; in practice, this situation can be avoided by limiting the maximum EOM modulation index or applying bandpass filters to the pulse shaper matrices.

Diagonal unitary decompositions in the form of Eq. (44) can be used in a variety of photonic DoFs, including position/momentum, parallel waveguides, and time bins-whenever the physical system can be modeled as the application of phase shifts in alternating Fourier-transform pairs. One can analytically design such systems by starting with the beamsplitter/phase-shifter decomposition of path encoding, and then expressing each beamsplitter layer as six alternating phase masks, for example. However, this may introduce significant resource overhead, which may present a challenge to computing the D_q matrix elements used to synthesize a desired target matrix U optimally—i.e., without an intermediate conversion step to an equivalent path circuit. Accordingly, numerical optimization can be used in some QFP implementations for basic gates such as the Hadamard, controlled-NOT, cyclic hop, and arbitrary single-qubit unitaries. In some examples described herein, the non-Gaussian CV cases considered deal with many-photon states inherently.

FIG. 3A and FIG. 3B provide an overview of an example non-Gaussian state engineering system. The case of N_s=3 input squeezed modes is shown for concreteness. FIG. 3A shows a schematic view of an example system 300, where squeezed states in distinct frequency modes from an input module 302 traverse the sequence of EOMs and pulse shapers in a QFP module 304, which is connected to an output module 306. In this example, the condition for successful heralding is the detection of n_s photons each in all but one of the central N_s bins and zero photons in all adjacent bins. The undetected mode is left in state |Φ. FIG. 3B shows a schematic view of the logic of system 300, represented by a transformation U, where each rail denotes an individual frequency bin, with the QFP functioning as a complex interferometer. As previously mentioned, our mathematical formalism applies to a system like FIG. 3A where the M undetected modes can be any of the total N modes without loss of generality. In the following example we choose M=1 and select this single undetected mode as the Kth mode, which is at the center of a set of N_s modes that are populated with single-mode squeezed vacuum states at the input; the remaining N−N_s modes are initially vacuum. For our simulations, we take the phase of the squeezing to be zero for all cases, although incorporating non-zero phases is possible in alternative examples. After application of U, the Kth mode is left undetected and n_s photons are detected in each of the other central N_s−1 modes. Production of the desired state in the Kth mode is heralded by simultaneously detecting vacuum in the remaining N−N_s modes: in essence, a bucket detector for all remaining modes. The use of such vacuum postselection is a consequence of the presence of an infinite set of ancilla modes in the frequency-bin DoF, which can be detected to ensure that the output state is pure.

Given the massive design space available for non-Gaussian QFP circuits—in terms of input states, unitaries, and output patterns—we have attempted in the specific configuration of FIG. 3A to provide a relatively simple example that nevertheless retains key features anticipated for successful circuits. By placing the output mode of interest in the center of the squeezed inputs, we increase opportunities for multiphoton interference with relatively weak modulation amplitudes, and selecting from the modes initially populated with photons for PNR detection with n_s>0 should permit reasonable success probabilities.

The use of PNR detection in the output module 306 after the QFP module 304 allows us to attain non-Gaussian states that are not available within other Gaussian cluster state models. In some alternative implementations, a GKP qubit encoding can be used that has a single photon in a discrete grid of spectro-temporal modes. While similar in that this also leverages the frequency DoF, many of the examples described herein provide for construction of non-Gaussian states in which quantum information resides in the field quadratures of optical modes, making our analysis multi-rather than single-photon in nature.

Having detailed the mathematical formalism and highlighted the specific features of the QFP, we now apply an example framework toward the design of quantum circuits that produce desired non-Gaussian output states, according to the configuration presented in FIG. 3A.

As noted previously, quantum system design with the QFP lacks an optimal analytical unitary decomposition procedure, so that numerical optimization may be used to obtain a QFP configuration realizing a desired unitary. In the context of non-Gaussian state design, the need for numerical optimization in itself is not unique, but has proven a fixture in path encoding as well.

In one approach to non-Gaussian state design, rather than designing a quantum circuit to implement some target state Φ_t directly, a Fock-truncated core state |Φ_core is sought instead—related to |Φ_t via a squeezing and displacement operation, |Φ_t=S ()D(β)|Φ_core. Suppose that the mode unitary found to produce |Φ_core is U; then, by absorbing the displacement and squeezing operation into a new set of inputs and mode unitary U′, an interferometer for the desired full state |Φ_t can be produced immediately via an analytical decomposition scheme. In the QFP case, however, if a set of EOM and pulse shaper solutions are found that can implement the core state preparation circuit U, there is no straightforward connection from it [i.e., the D_q matrices in Eq. (44)] to the modified configuration that would realize U′; instead, numerical optimization may be employed on U′, effectively doubling the rounds of numerical design compared to path encoding. In what follows we concentrate on synthesizing circuits that produce the full target state |Φ_t immediately, avoiding this intermediate core state step.

To begin the optimization process we first define the target state |Φ_t in the Fock basis, i.e., the coefficients τ_n=n|Φ_t. We employ MATLAB's particle swarm optimization (PSO) tool to find the N_s nonzero input squeezing values of the total length-N vector of inputs

$\begin{array}{cc}\overrightarrow{r}=\left(0,\dots ,0,r_{K-\text{?}\frac{N_s}{2}]},\dots ,r_{K-1},r_K,r_{K+1},\dots ,r_{K+\lfloor \frac{N_s}{2}\rfloor },0,\dots ,0\right),& \left(45\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

and a QFP unitary, U, that when applied to the N-mode input followed by detection of n_s photons in each of the remaining N_s−1 squeezed modes, produces a state |Φ. Letting n_c denote the photon number at which we truncate the state for numerical simulations, we therefore compute a total of n_c+1 coefficients (including vacuum) to fully describe the heralded output. To find the optimal squeezing values and U, PSO varies the phase shifts applied to the N QFP modes by the pulse shapers, each EOM's phase modulation function φ(t), and the N_s nonzero elements of {right arrow over (r)} in order to minimize the cost function

$\begin{array}{cc}C=P{\log }_{10}\left(1-ℱ\right),& \left(46\right)\end{array}$

where and P are the fidelity of |Φ with respect to |Φ_t and the probability of producing |Φ, respectively. We have found a logarithmic cost function of this form useful for penalizing <1 more strongly than P<1, emulating the effect of a constraint on without the computational cost associated with a strict constraint function. Other cost functions may also be employed without materially altering our approach. The Fock coefficients of |Φ in the Kth mode can be expressed as Eq. (40) which we write in the form

$\begin{array}{cc}{c}_{n_K}=\frac{I_{\overrightarrow{n}}}{\sqrt{P}\prod _{i=1}^N\sqrt{n_i!2^{n_i}cosh{r}_i}},& \left(47\right)\end{array}$

where {right arrow over (n)}=(0, . . . , 0, n_s, . . . , n_s, n_K, n_s, . . . , n_s, 0, . . . , 0) is the vector of photon numbers over all output modes, so that P and can be written as

$\begin{array}{cc}P=\sum _{n_K=0}^{n_c}\left(\frac{{❘I_{\stackrel{\rightharpoonup }{n}}❘}^2}{\prod _{i=1}^N{n}_i!2^{n_i}\cosh r_i}\right)\mathrm{and}& \left(48\right)\end{array}$ $\begin{array}{cc}ℱ={❘\langle {\Phi }_t|\Phi \rangle ❘}^2={❘\sum _{n_K=0}^{n_c}{\tau }_{n_K}^{*}{c}_{n_K}❘}^2& \left(49\right)\end{array}$

With the cost function defined we now lay out an example procedure for evaluating and P at each PSO iteration. First, we calculate

${\left(V_0+\frac{I}{2}\right)}^{-1}$

using {right arrow over (r)} in Eq. (22). U is calculated by plugging the N phase shifts for each pulse shaper and each EOM's φ(t) into Eq. (44), which we convert to symplectic form, S_p, via

$\begin{array}{cc}{S}_p=W\left(\begin{array}{cc}U& 0\\ 0& U^{*}\end{array}\right)W^{†},& \left(50\right)\end{array}$

where W is defined by Eq. (29) and U* corresponds to element-by-element conjugation (no transpose). Γ⁻¹ is then calculated by Eq. (16), and the blocks A and extracted per Eq. (26). A and are used to find with Eq. (35). The matrix elements of σ are found by using in Eq. (37). Because we detect vacuum in all the QFP modes except for the center N_s,

${\overrightarrow{h}}^T=\left(0,\dots ,0,\frac{n_s}{2}-v_{K-\lfloor \frac{N_s}{2}\rfloor },\dots ,\frac{n_K}{2}-v_K,\dots ,\frac{n_s}{2}+v_{K-\lfloor \frac{N_s}{2}\rfloor },0,\dots ,0\right)$

in Eq. (41) renders unimportant all the elements of σ other than the center N_s×N_s block. Therefore we proceed to evaluate Eq. (47) using only the center block of σ, for n_K∈{0, 1, . . . , n_c} and are left with the Fock coefficients c_nK of |Φ.

We choose to compute Γ⁻¹ in this manner for computational reasons. As a large matrix—2N×2N in general and 128×128 in our case—Γ is time-consuming to invert. We bypass this time sink by calculating Γ⁻¹ directly with Eq. (22), rather than performing Γ=S_pV₀S_p^T+I/2 and inverting Γ.

We take further action to streamline the n_c+1 calculations of I_{{right arrow over (n)}} needed to find the Fock coefficients of |Φ, dominated by Hf(σ) in Eqs. (9) and (41). σ is the only quantity in Eq. (41) that will change in the successive iterations of PSO; therefore, we can precompute a number of the elements of Eq. (41) outside the optimization loop and use them for every PSO iteration. Calculating these elements upfront proves imperative to expediting the optimization process when n_c becomes large.

Consider a given value n_K. First, we define the length-N_s vector

$\begin{array}{cc}{\overrightarrow{s}}^T=\left(n_s,\dots ,n_s,n_K,n_s,\dots ,n_s\right)& \left(51\right)\end{array}$

the sum of photons in the output modes

$\begin{array}{cc}\Sigma =\sum _{i=1}^{N_s}{s}_i=\left(N_s-1\right)n_s+n_K,& \left(52\right)\end{array}$

and index vectors for each mode {right arrow over (v)}_i^T=(0, 1, . . . , s_i), where the maximum s_i for each mode is taken from Eq. (51).

Then we find all combinations of the entries of the v vectors and store them in a matrix D, where each row corresponds to a unique length-N_s listing of elements, one drawn from each {right arrow over (v)}_i. Because we choose to detect the same number of photons, n_s, in the N_s−1 modes, D will be of dimension (n_s+1)^Ns−1(n_K+1)×N_s and will take on the role of the nested summations that appear in Eq. (41). We calculate the exponent of the (−1) factor in Eq. (41) for all terms of the nested summation and store them in whose elements are defined as

$\begin{array}{cc}{𝒲}_i=\sum _{j=1}^{N_s}{D}_{\mathrm{ij}}.& \left(53\right)\end{array}$

Similarly, the product of binomials in Eq. (41) is calculated for all terms in the nested summation and stored in {right arrow over (X)},

$\begin{array}{cc}{X}_i=\prod _{j=1}^{N_s}\left(\begin{array}{c}{s}_j\\ D_{\mathrm{ij}}\end{array}\right).& \left(54\right)\end{array}$

{right arrow over (h)}^T for all the terms in the nested summation are stored in vectors

$\begin{array}{cc}{\overrightarrow{Z}}_i^T=\left(\frac{s_1}{2}-D_{i1},\dots ,\frac{s_{N_s}}{2}=D_{{\mathrm{iN}}_s}\right).& \left(55\right)\end{array}$

The Hf(σ) calculation is then reduced to a single summation over these precomputed elements,

$\begin{array}{cc}\mathrm{Hf}\left(\sigma \right)=\sum _{i=1}^{\kappa }\frac{1}{\left(\frac{\Sigma }{2}\right)!}{\left(-1\right)}^{{𝒲}_i}{X_i\left(\frac{1}{2}{\overrightarrow{Z}}_i^T\sigma {\overrightarrow{Z}}_i\right)}^{\frac{\Sigma }{2}},& \left(56\right)\end{array}$

where κ=(n_s+1)^Ns−1 (n_K+1). This process can be repeated for all values n_K∈{0, 1, . . . , n_c} giving us all necessary precomputed {right arrow over (s)}, Σ, {right arrow over (W)}, {right arrow over (X)}, and {right arrow over (Z)}_i, elements.

As examples of our method, we seek to generate even Schrödinger cat states with coherent amplitudes α ranging from 0.5 to 3 in steps of 0.25. For each α value |Φ_t is therefore set to

$\begin{array}{cc}❘{\Phi }_t\rangle =\frac{❘\alpha \rangle +❘-\alpha \rangle }{\sqrt{2\left(1+e^{-2|\alpha {|}^2}\right)}},& \left(57\right)\end{array}$

where

$❘\pm \alpha \rangle \approx e^{-\frac{1}{2}{❘\pm \alpha ❘}^2}{\sum }_{n=0}^{n_c}\frac{{\left(\pm \alpha \right)}^n}{\sqrt{n!}}❘n\rangle .$

We truncate |±α at n_c=40 for all α values, which encompasses all the Fock support to high precision at α=3, and therefore for any |Φ_t with α<3 as well. Indeed, the truncation error defined as ε_nc=1−Σ_n=0^nc|n|Φ_t |² is less than 10⁻¹⁴ for α≤3 and n_c=40. This choice is highly conservative, as one could likely consider smaller n_c values such as n_c=20 or 30 for added computational speed up, for which the errors remain small: ε₂₀<10⁻³ and ε₃₀<10⁻⁸ at α=3.

As an example, to make these results as tractable as possible for implementation, we limit φ(t) to a single sinewave and constrain

$\left(r_{K-\lfloor \frac{N_s}{2}\rfloor },\dots ,r_K,\dots ,r_{K+\lfloor \frac{N_s}{2}\rfloor }\right)$

to a maximum of 1.5 (corresponding to a squeezing value of approximately 13 dB). We proceed to optimize with Q∈{3,5, 7} total QFP elements, N=64 QFP modes, and N_s∈{3,5} input squeezed states, along with a 32-mode bandpass filter on each pulse shaper to prevent unphysical solutions that reach the edge of the N=64-mode truncation. The nonzero PNR detectors are set to herald on n_s=1, which ensures that Σ in Eq. (9) will be even when computing even Fock coefficients in the undetected mode K (n_K∈{0, 2 . . . , 40}). The target cat state coefficients are real numbers; however, the coefficients found by optimization are in general complex. Therefore if the state found by optimization has fidelity equal to one, it should have a constant phase for all Fock coefficients.

FIG. 4A, FIG. 4B and FIG. 4C show, for prophetic examples, wavefunctions in the quadrature basis q|Φ (top panels) and photon number probabilities |n|Φ |² (bottom panels) for example target and QFP output states. FIG. 4A has parameters: α∈{1, 1.5, 2}, Q=3, and N_s=3. FIG. 4B has parameters: α=2.25, Q∈{3, 5, 7}, and N_s=5. FIG. 4C has parameters: α=1.75, Q=3, and N_s∈{3,5}.

To elucidate how the size of the cat state changes with α we plot, in FIG. 4A, |Φ_t (target) and |Φ (circuit), with N_s=3 and Q=3, for α∈{1, 1.5, 2}. The plots in FIG. 4B and FIG. 4C illustrate how the quality of |Φ changes with N_s and Q for single α values, whereas FIG. 5 shows the overall trends. While running PSO it became apparent that our chosen cost function [Eq. (46)] was not ideal in that the |Φ with the lowest was not always the “best” state as viewed by our intuition. For example, in FIG. 5 for α=2.5, Q=3, and N_s=5, the output |Φ with the lowest is a state with <<1. Consequently, we include in FIG. 5 not only the states |Φ with the lowest cost , but also higher cost solutions, found with different initial conditions, that “pass the eye test” (corrected). We emphasize that this distinction does not reflect any issues in the optimization procedure itself, but rather in our selection of the cost function; to encourage PSO to find even higher fidelity states, future tests could increase the base of the logarithm beyond 10 in Eq. (46).

FIG. 5 shows, for prophetic examples, plots of fidelity (top panels) and success probability (bottom panels) plotted as functions of α (in steps of 0.25) for various combinations of N_s and Q. Corrected results are shown for the N_s=5, Q=3 casc. It is evident from FIG. 5 that for a set amount of resources, constant N, and Q, the output state quality found by PSO finds decreases as α increases. This can be attributed to the fact that |Φ_t becomes noticeably more non-Gaussian as α is increased, see FIG. 4A. FIG. 5 also shows that while increasing the complexity of the QFP through the number of elements Q can moderately improve the success probability, see FIG. 4B, it does not lead to markedly higher fidelities in these examples. In contrast, for any Q, the addition of more ancilla resource states (larger N_s) can substantially improve fidelity, particularly for larger values of a, albeit with about an order of magnitude reduction in success probability, see FIG. 5. This behavior makes sense; the extra photons available provide a greater variety of interference possibilities in design, yet also reduce the success probability through additional PNR detector conditions that are to be satisfied.

Single-mode squeezed vacuum states are especially convenient because of their zero displacement in phase space ({right arrow over (x)}_β=0) and diagonal covariance matrix V₀. The latter facilitated closed-form expressions for det Γ [Eq. (14)] and Γ⁻¹ [Eq. (26)] can simplify calculations for the numerical optimizer. However, the covariance matrix just before the partial PNR detection (i.e., just after the passive Gaussian unitary operator U of FIG. 2) nevertheless remains completely general for a pure Gaussian state. Indeed, the covariance matrix of a pure Gaussian state is of the form SV_vacS^T, where V_vac is the covariance matrix of vacuum and S is any symplectic matrix which includes squeezing. i.e., S=S_pS_sS_p′, where S_p and S_p′ are passive transformations and S_s is the symplectic transformation for squeezing. Since passive transformations have no effect on vacuum, the most general covariance matrix for a pure Gaussian state can thus be written as SV_vacS^T=S_pV₀S_p^T, which is precisely the covariance matrix assumed in our analysis. For example, any two-mode squeezed state (like those produced in quantum frequency combs described herein) can be expressed as the interference of two single-mode squeezed states on a beamsplitter, whose unitary can be carefully incorporated in the front-end of the QFP module 304 in FIG. 3A.

Although the diagonal input covariance matrix V₀ does not reduce the generality of our formulation, the absence of displacement is potentially significant. Incorporating nonzero displacements will not affect the covariance matrix we have used; it will, however, introduce additional variables into the optimization procedure for generating desired non-Gaussian states.

It is possible to move beyond unitary operations as well. For example, by coupling each frequency bin to additional environmental modes, then tracing these out, photon loss can be added into the formulation. Such an extension can be useful for some implementations of a QFP.

A long-standing challenge in CV encoding is the realization of GKP qubit states for error correction. An example of the value of GKP qubits, call them |0 and |1 in the logical basis, lies in their infinite series of equispaced delta functions, |1 being displaced from |0 by √{square root over (π)} when plotted in the q-quadrature basis. Since these ideal states are unphysical, approximate states |{tilde over (0)} and |{tilde over (1)} were presented in the original GKP proposal, which consist of a sum of Gaussian peaks with standard deviation Δ, all under another Gaussian envelope with standard deviation

$\frac{1}{k}.$

Δ=k=0.15 is used for |{tilde over (0)} and |{tilde over (1)} to maintain a 99% error correction rate. A determination of how n_s and ancilla mode placement affect the quality of the output state can aid in finding effective QFP circuits for direct GKP state production.

An alternative path to quality approximations of GKP qubits, for which our system is already well suited, is the so-called “cat breeding” protocol. In some examples of the protocol, two cat states are squeezed by some amount r, where r=−lnΔ. The squeezed cat states are combined on a balanced beamsplitter, and a homodyne measurement is made on one of the output modes. When the result of the homodyne measurement for a single output mode's p-quadrature is zero, the other output mode is left in a state with three equispaced peaks. The height of the peaks follows a binomial distribution, and the width of the peaks is determined by the amount of squeezing applied to the initial cat states. Successive iterations of the protocol, where the beamsplitter inputs are the states produced by the previous iteration, yield higher-order binomial states. To ensure that the final state has the correct spacing associated with GKP states, the initial cat states have a coherent amplitude

$\alpha ={\sqrt{2}}^{m-1}\sqrt{\pi }{e}^r,$

where m is the number of iterations of the protocol to be executed. The larger r and m are, the more closely the resulting state will resemble the approximate GKP state, making access to large cat states vital to the protocol. The cat breeding protocol can be modified to replace the homodyne measurement by PNR detection. Because PNR detection neglects the phase of the output, fine control over the relative phase of the input states may be used to achieve the same comb-like output as in the homodyne approach. In this example, by detecting four photons at one output mode after a single iteration of the protocol, numerically generated states with a fidelity of 0.996 with respect to an approximate GKP state (Δ=k=0.545) at a success probability of 0.09.

In this example, our system can generate cat states up to a size α=2 with 99.87% fidelity when N_s=5 and Q=7. These capabilities make our non-Gaussian state engineering system a viable candidate to meet the resource state demands set by cat breeding protocols.

In summary, we have described a system for the production of non-Gaussian quantum states using the QFP, a device designed to implement arbitrary linear-optic transformations on discrete spectral modes. Our formulation using the K function expansion enables efficient calculation of multimode Gaussian states in the photon-number basis, providing a valuable framework for analysis in any photonic DoF. Applying this to the QFP specifically, we have designed quantum circuits that produce non-Gaussian cat states with a variety of amplitudes, revealing a clear fidelity/success-probability tradeoff with the number of squeezed ancillas. The techniques described herein provide tools for designing CV quantum systems in frequency bins, and facilitate the ultimate realization of fiber-compatible, single-spatial-mode, and massively parallelizable quantum information processors based on non-Gaussian photonic states.

The techniques described herein enable implementation of a compact photonic source for arbitrary non-Gaussian states of light at a selected frequency (e.g., a selected frequency of a frequency comb source). For example, a multi-frequency-mode (spectral comb type) light source (in some cases a squeezed light source) can be used at the input. A QFP receives the input and provides a frequency-domain linear-optic unitary transformation enacted by iterating between a dispersive element and a phase modulator. At the output of the QFP a frequency-dependent element and PNR detectors can enable spectrally-resolved photon number resolving detection on frequencies other than the selected frequency (e.g., all but the selected frequency). The undetected frequency is where a particular generated quantum state is produced, in the event that a particular “click pattern” occurs on the PNR detectors.

Thus, in a compact form-factor (e.g., in a photonic integrated device), the system can provide a non-Gaussian quantum state that can be used for any of a variety of applications. For example, the source can provide high-fidelity high-rate GKP qubits, or high-quality photonic cat states, which can be used for quantum computing or all-photonic quantum repeaters. Quantum optical receiver systems for optical communications or sensing can be implemented by realizing non-Gaussian operations on a collection of optical modes, by generating off-line non-Gaussian states, and using them in conjunction with linear optics and photon detectors, to act on the received modulated signal (e.g., for communications) or on the target-return signal (e.g., sensing).

While the disclosure has been described in connection with certain embodiments, it is to be understood that the disclosure is not to be limited to the disclosed embodiments but, on the contrary, is intended to cover various modifications and equivalent arrangements included within the scope of the appended claims, which scope is to be accorded the broadest interpretation so as to encompass all such modifications and equivalent structures as is permitted under the law.

Claims

1. An apparatus for generating a non-Gaussian quantum state associated with one or more spectral modes, the apparatus comprising: a quantum optical frequency comb source comprising at least one nonlinear optical medium and configured to provide a plurality of spectral modes spaced at a frequency-bin spacing, where the spectral modes include a plurality of pairs of entangled spectral modes, each pair of entangled spectral modes including: a first spectral mode centered at a first frequency, anda second spectral mode entangled with the first spectral mode and centered at a second frequency that is spaced from the first frequency at a multiple of the frequency-bin spacing; anda plurality of electrically controllable optical transformation modules connected in series with a first module in the series receiving spectral modes from the quantum optical frequency comb, where two or more of the modules each comprise: a spectral mode phase shifter configured to apply respective phase shifts to different spectral modes based on at least a first electrical signal, anda spectral mode mixer coupled to the spectral mode phase shifter and configured to couple spectral modes centered at different frequencies based on at least a second electrical signal.

2. The apparatus of claim 1, wherein one or more of the electrically controllable optical transformation modules are configured to transform each pair of entangled spectral modes into single-mode squeezed vacuum states.

3. The apparatus of claim 1, wherein each spectral mode phase shifter comprises a Fourier-transform pulse shaper.

4. The apparatus of claim 1, wherein each spectral mode mixer comprises an electro-optic phase modulator.

5. The apparatus of claim 1, further comprising an output interface coupled to a last module in the series, the output interface including: a plurality of photon number resolving (PNR) detectors configured to detect respective spectral modes of more than one and fewer than all spectral modes output from the last module in the series, andat least one port providing at least one output spectral mode from the last module in the series not detected by any of the PNR detectors.

6. The apparatus of claim 5, wherein the output interface further comprises a wavelength-dependent element between the last module and the plurality of PNR detectors.

7. The apparatus of claim 6, wherein the wavelength-dependent element comprises at least one of: a grating or a prism.

8. The apparatus of claim 5, wherein each PNR detector is configured to generate a detection signal that distinguishes a detected photon number equal to zero from a detected photon number equal to one, in each of a plurality of time slots; anddistinguishes a detected photon number equal to one from at least one detected photon number greater than one, in each of the plurality of time slots.

9. The apparatus of claim 8, further comprising a trigger module configured to generate a trigger signal based on detection signals from the plurality of PNR detectors indicating one or more time slots in which the port providing at least one output spectral mode corresponds to a non-Gaussian quantum state.

10. The apparatus of claim 9, wherein the non-Gaussian quantum state comprises a quantum superposition of continuous variable (CV) optical wave quadrature states.

11. The apparatus of claim 10, wherein the quantum superposition comprises a superposition of optical waves with opposite phases.

12. The apparatus of claim 9, wherein the non-Gaussian quantum state is input to a module that includes one or more optical elements that provide at least one Gottesman-Kitaev-Preskill (GKP) qubit.

13. The apparatus of claim 12, wherein the one or more optical elements include at least one optical beamsplitter.

14. The apparatus of claim 9, wherein the trigger signal is based on detection signals that include a plurality of the detection signals each indicating a photon number of zero in at least one time slot.

15. The apparatus of claim 5, wherein the last module in the series consists essentially of a spectral mode mixer configured to couple spectral modes centered at different frequencies based on an electrical signal.

16. The apparatus of claim 15, wherein the first module in the series consists essentially of a spectral mode mixer configured to couple spectral modes centered at different frequencies based on an electrical signal.

17. The apparatus of claim 1, wherein a plurality of the modules in the series are integrated on a common photonic integrated circuit.

18. The apparatus of claim 1, wherein the frequency-bin spacing is a free spectral range of an optical parametric oscillator that includes the nonlinear optical medium.

19. A method for generating a non-Gaussian quantum state associated with one or more spectral modes, the method comprising: providing from a quantum optical frequency comb source, comprising at least one nonlinear optical medium, a plurality of spectral modes spaced at a frequency-bin spacing, where the spectral modes include a plurality of pairs of entangled spectral modes, each pair of entangled spectral modes including: a first spectral mode centered at a first frequency, anda second spectral mode entangled with the first spectral mode and centered at a second frequency that is spaced from the first frequency at a multiple of the frequency-bin spacing; andreceiving spectral modes from the quantum optical frequency comb into a first module of a plurality of electrically controllable optical transformation modules connected in series, where two or more of the modules each comprise: a spectral mode phase shifter configured to apply respective phase shifts to different spectral modes based on at least a first electrical signal, anda spectral mode mixer coupled to the spectral mode phase shifter and configured to couple spectral modes centered at different frequencies based on at least a second electrical signal.

20. An apparatus for generating a non-Gaussian quantum state associated with one or more spectral modes, the apparatus comprising: a quantum optical frequency comb source comprising at least one nonlinear optical medium and configured to provide a plurality of spectral modes spaced at a frequency-bin spacing,a plurality of electrically controllable optical transformation modules connected in series with a first module in the series receiving spectral modes from the quantum optical frequency comb, where two or more of the modules each comprise: a spectral mode phase shifter configured to apply respective phase shifts to different spectral modes based on at least a first electrical signal, anda spectral mode mixer coupled to the spectral mode phase shifter and configured to couple spectral modes centered at different frequencies based on at least a second electrical signal; andan output interface coupled to a last module in the series, the output interface including: a plurality of photon number resolving (PNR) detectors configured to detect respective spectral modes of more than one and fewer than all spectral modes output from the last module in the series, andat least one port providing at least one output spectral mode from the last module in the series not detected by any of the PNR detectors;wherein each PNR detector is configured to generate a detection signal that: distinguishes a detected photon number equal to zero from a detected photon number equal to one, in each of a plurality of time slots, anddistinguishes a detected photon number equal to one from at least one detected photon number greater than one, in each of the plurality of time slots.

21. The apparatus of claim 20, wherein the number of spectral modes spaced at a frequency-bin spacing is three or more.

22. The apparatus of claim 20, further comprising a trigger module configured to generate a trigger signal based on detection signals from the plurality of PNR detectors indicating one or more time slots in which the port providing at least one output spectral mode corresponds to a non-Gaussian quantum state.

23. The apparatus of claim 20, wherein the plurality of spectral modes provided by the quantum optical frequency comb source comprise a Gaussian state.

24. The apparatus of claim 20, wherein the plurality of spectral modes provided by the quantum optical frequency comb source comprise a state continuous-variable encoded state.