Taming numerical errors in simulations of continuous variable non-Gaussian state preparation.

Numerical simulation of continuous variable quantum state preparation is a necessary tool for optimization of existing quantum information processing protocols. A powerful instrument for such simulation is the numerical computation in the Fock state representation. It unavoidably uses an approximation of the infinite-dimensional Fock space by finite complex vector spaces implementable with classical digital computers. In this approximation we analyze the accuracy of several currently available methods for computation of the truncated coherent displacement operator. To overcome their limitations we propose an alternative with improved accuracy based on the standard matrix exponential. We then employ the method in analysis of non-Gaussian state preparation scheme based on coherent displacement of a two mode squeezed vacuum with subsequent photon counting measurement. We compare different detection mechanisms, including avalanche photodiodes, their cascades, and photon number resolving detectors in the context of engineering non-linearly squeezed cubic states and construction of qubit-like superpositions between vacuum and single photon states.

Introduction

Quantum information theory exploits fundamental features of quantum physics to design protocols and algorithms that offer significant improvements over their classical counterparts[1-4]. There are several candidate physical systems suitable for these applications, each with distinct advantages. Continuous variable quantum information processing with light offers feasible and fast generation and manipulation of entangled Gaussian quantum states that are at the core of the information protocols[5-12]. However, truly universal quantum information processing also requires elements of quantum non-Gaussianity[13-17]. Protocols based on Gaussian states and Gaussian operations are not universal[13] and can be efficiently simulated on a classical device[18].

For continuous variables of light, the non-Gaussianity is commonly introduced by photon number counting detectors, either the most basic on-off detectors capable of discerning presence of light[19], or the more advanced detectors truly distinguishing the photon numbers[20-29]. Such detectors can be employed for direct conditional implementation of non-Gaussian operations[30-37], or for conditional preparation of non-Gaussian quantum states[38-50]. The latter can be then used as a resource in deterministic implementation of non-Gaussian gates[14,37,51]. One thing these approaches have in common is the inherent probabilistic nature of measurement that results in several trade-offs between quality of the implemented operation or the prepared quantum state, the rate with which the desired operation succeeds, and the experimental challenges of the photon number resolving detector[26-29,52]. For any given set of realistic detectors and any desired task we then need the ability to faithfully simulate the optical circuit to find out the required parameters leading to the optimal performance, or to find out whether the task is even feasible.

However, numerical simulation of simple quantum optical circuits, even though it is often employed in continuous variable quantum information processing[53-57], is not a straightforward task. It is burdened by various difficulties, including discretization errors in numerical models relying on continuous representation, truncation errors in discrete models[53], the omnipresent rounding errors due to finite precision of arithmetics[58-62] and numerical truncation errors occurring in finite approximations of infinite processes[60,62]. If not prevented by rigorous analysis, these numerical artifacts can dominate the computed values and lead to rapid divergence from correct results.

In this paper we evaluate the numerical errors arising when an optical circuit for probabilistic preparation of non-Gaussian quantum states of light[14,38] is simulated on a classical digital computer. We then propose an alternative method for construction of truncated unitary operators aiming to curtail these errors. Finally we take advantage of these tools to fully simulate the circuit for preparation of resource states for the cubic phase gate[51], and single mode qubit-like superpositions of zero and one photon. The goal is to find the optimal trade-offs between the quality of the states and the probability of success for a range of available photon counting detectors[20,26-29,52].

This paper is structured as follows. In the State preparation circuit section we review the state preparation circuit. In the Perils of numerical simulation of CV systems section we describe the errors naturally occurring in numerical simulations. In the The curious case of coherent displacement section we focus on coherent displacement and identify the numerical errors appearing in different methods of its calculation. In the Truncated approximate matrix exponential (TAME) section we propose an alternative method for its calculation, followed by an overview of verification process in the Verification of approximated matrices section. We then proceed with the Numerical simulation of the preparation circuit section, where we describe the methodology of the actual simulation and present the results of its applications in sections eight and nine.

State preparation circuit

The most common method of conditional state preparation is based on suitable manipulation of ENPR state with coherent displacement and subsequent photon counting measurement[14,39,63,64]. In Fig. 1 we present a variant of the circuit which can be used for preparation of simple non-Gaussian quantum states, including the qubit-like and superpositions. Our circuit accounts for basic imperfections limited to detection inefficiencies and propagation losses. A physical EPR resource, generating the two mode squeezed vacuum (TMSV), lies at its very heart and serves as a source of perfectly correlated photons. One of the entangled modes is then displaced with controllable amplitude and phase and consequently measured. The detection can use either an avalanche photodiode (APD), a photon number resolving detector (PNRD) or its approximation employing an APD cascade[52]. The resulting marginal stateconditioned on the detection outcome , characterized by the POVM element , is obtained with the probability of successIn both the expressions (1) and (2) we use the lower right indices to emphasize which modes the operators and channels act on. Starting from the inner-most component, the initial TMSV state is denoted with with coefficients , where the parameter  sets the experimentally controllable squeezing strength. We model the overall losses and inefficiencies in the preparation scheme as attenuation of the measured mode prior to its displacement. This can represented by a Gaussian quantum channel with its action on the mode given in terms of Kraus operators[65] as with , where defines the photon number operator and denotes the annihilation operator respective to the measured mode. The parameter  describes the efficiency of the preparation circuit. Subsequently the converse  characterizes the overall losses and inefficiencies in the preparation scheme. The displacement of the second mode is given by the unitary operator , where is the displacement amplitude[66]. In a more realistic analysis of the preparation circuit it would be straightforward to include the propagation losses affecting the mode carrying the resulting state. This form of decoherence can be accounted for by modifying the squeezing strength of the non-linearly squeezed state[67]. Consequently we do not consider this additional attenuation since it does not influence the fundamental properties of these non-Gaussian states.

Figure 1

Variation of the conditional preparation scheme. We start with a two mode squeezed vacuum state . One of its modes is then displaced with and measured, using either APD, PNRD or an APD cascade approximating PNRD. The detection outcome is characterized by the POVM element . We model overall losses and inefficiencies within the scheme using a beam splitter with intensity transmittance  to represent attenuation of the signal state in the setup.

From the experimental perspective the parameters and can be fine tuned to engineer a desired state with optimal performance given particular experimental configuration characterized by the efficiency and conditioning on the detection outcome with respective POVM element .

We can utilize this scheme to prepare a variety of quantum states. Consider now a lossless configuration employing an ideal PNRD. Its POVM elements correspond to projectors onto individual Fock states . The output state, conditioned on the detection of a particular Fock state , is then proportional to where the coefficients follow from the definition of the TMSV state and are matrix elements of the displacement operator. By tuning the parameters and we can construct a set of states parametrized by the possible combinations of the and coefficients.

Possible applications of this scheme include construction of generally non-classical superpositions of Fock states[38,39] and, in particular, non-linearly squeezed non-Gaussian states[14,63]. Every application can be translated into constrained optimization of the tunable parameters with the constraint and objective functions embodying the nature of the particular application.

For example, if one were to construct a specific state , the optimization objective could be to maximize some metric of similarity with the target state, e.g., fidelity. It would also be practical to construct the state with non-negligible probability of success. This requirement could be expressed either as a constraint allowing only solutions with the probability greater than some threshold, or as an additional optimization objective in multi-objective optimization.

The constraint and objective functions generally involve the success probability (2) and the resulting density operator (1). Both of which can be obtained by simulating the preparation procedure numerically on a classical digital computer. But alas, numerical simulations come with their own hurdles which will be identified and subsequently addressed in the following sections.

Perils of numerical simulation of CV systems

Classical digital computers[68] encode information into finite sequences of bits and it is therefore impossible to represent arbitrary real numbers. The standard approach[59-62,69] is to approximate real numbers with floating point (FP) numbers. Real numbers are then rounded to their closest representable FP neighbors. This generally introduces rounding errors. To make matters worse, FP arithmetic with FP numbers does not necessarily produce exactly representable floating point numbers. Results of FP arithmetic must be rounded, possibly introducing additional rounding errors[59-62,69]. Consequently complex sequences of arithmetic operations possess the potential to accumulate and even amplify rounding errors. Even the most straightforward tasks such as adding up a sequence of FP numbers can produce widely different results with varying degrees of accuracy based on the algorithm of choice[59-61]. Rounding error analysis is therefore a crucial part of algorithm design[58-61] and commonly used numerical algorithms are frequently accompanied by rigorous rounding error analysis. Nevertheless numerical simulation cannot be considered completely accurate as the error analysis only establishes upper bounds on the numerical errors[58,60-62].

The practical concerns, when dealing with numerical simulation, are therefore always related to size of the errors, rather than to their presence. This is a familiar concept in physics, a discipline which is well acquainted with limited precision of measured quantities[70,71]. Numerical simulation of CV systems suffers from further issues related to the fundamental representation of quantum states and quantum operations. CV states reside in infinite-dimensional Hilbert spaces and can be, in principle, described in two distinct ways. The first description employs continuous functions, either wave functions given in position or momentum representation, or quasi-probability distributions[5-7] which combine the two quadratures. The practical issue with this approach is the continuous nature and generally infinite support of these functions as their support must be limited to finite intervals and both their domains and ranges discretized during numerical integration[61,62], introducing additional numerical errors.

Alternatively we can take the advantage of the discrete Fock basis spanned by eigenstates of the number operator. This basis is still infinite but, unlike in the case of basis spanned by eigenstates of continuous operators, the number of its elements is countable. While exact representation of CV states in Fock basis remains impossible, we can truncate the basis to a finite number of elements and approximate the original Hilbert space with this truncated, finite-dimensional, restriction. We can thusly avoid discretization errors and deal with truncation errors instead. Consequently numerical simulations utilizing truncated Hilbert spaces spanned by truncated Fock basis are often employed in detailed analysis of CV quantum circuits.

Formal definition of truncated Fock spaces

Let denote the original Hilbert space and let be the original Fock basis (FB). In this basis the vector components of individual Fock states satisfy , that is, Fock states form an orthonormal basis. We take the first elements of FB, and truncate their vector forms to the first components, forming the truncated Fock basis (TFB) where we use the upper right indices in to denote dimensions of said vectors. Vector components of TFB elements satisfy . The basis therefore remains orthonormal. The linear hull of forms the dimensional truncated Fock space (TFS) .

So far we have only defined TFS itself and the transition from FB to TFB. In the following we define the transition of vectors from into and linear operators from to . Let be an arbitrary state expressed as (where ) with coefficients . The expression (where ) then defines its truncated variant from . Let be a linear operator on expressed as (where ) with matrix elements . Then (where ) defines its truncated analogue on . A natural extension of this approach allows for transitions from higher-dimensional spaces to lower-dimensional spaces.

Navigating truncated Fock spaces

In this description, pure quantum states become complex dimensional vectors of numbers, linear operators turn into complex matrices and the operations we would otherwise perform, reduce to linear algebraic expressions such as matrix multiplication, Kronecker products and matrix traces. There is, however, a hefty price to be paid for this simplification, manifesting in the form of truncation errors with several distinct effects on the simulation.

Firstly, it is impossible to represent general quantum states exactly. Take an arbitrary quantum state and its truncated variant . The quality of the truncated state can be determined from its normalization, or rather the lack of it, using the cutoff errorwhere are the vector components of the state in Fock representation. In essence the quality of the representation is loosely given by the support of the state relative to the dimension of the TFS. This is not the only conceivable metric, but it is a convenient one as it is straightforward to calculate.

Secondly, the algebraic structure of the space changes with the transition to finite dimension. As a result the usual commutation rules no longer apply since for any pair of operators and the relation does not necessarily hold anymore. We can illustrate the change in algebraic structure on bosonic creation and annihilation operators. In the regular infinite-dimensional case we have , that is, the two operators commute to identity. With the truncated commutator the result remains the same , which is an identity matrix of the corresponding dimension . Conversely the commutator of the truncated annihilation and creation operators differs from identity in the final element on the diagonalwhich can be understood as a truncation error due to the product of two truncated matrices.

Thirdly and finally, replacing infinite-dimensional operators in arguments of operator functions with their truncated versions may not be without consequences. Consider an operator function . In principle for general operator arguments. This has grave consequences for numerical simulation of unitary evolution. It is customary to approximate the exponential operator, , with the matrix exponential of the truncated operator argument[53]. However, this method can not be relied upon as . We must therefore seek alternative approaches: there are three primary techniques available for numerical simulation. The first one relies on the knowledge of a closed form formula for elements of the unitary operator. It has to be derived analytically and is not always attainable. The second method, proposed in the recent paper[53], is numerical and derives individual elements of unitaries by recurrent formulae. In the third approach the matrix exponential is simply computed with the truncated matrix argument as and the dimension of the computation space is chosen large enough so that the errors are irrelevant in the particular simulation.

Neither approach is perfect. Each suffers from specific numerical errors. This is a valid concern even for the first method which uses analytical forms: it is because mathematical expressions, especially those involving factorials, large powers of non-negligible numbers or relying on special functions, which are often defined using similar expressions or recurrent formulae, still need to be evaluated numerically with finite precision in floating point arithmetic, leading to introduction and eventual accumulation of rounding errors. The numerical errors cannot be straightforwardly calculated without a priori knowledge of the ideal operator or without thorough numerical analysis of rounding errors, an area of expertise that is mostly out of the scope of theoretical physics and therefore scarcely present in research reports.

In the following section we apply these methods of construction to the simplest experimentally testable example, coherent displacement, and use this particular case study to demonstrate the fundamental shortcomings of each approach.

The curious case of coherent displacement

Coherent displacement is a fundamental Gaussian operation in quantum optics used in a broad range of quantum protocols for quantum state preparation, manipulation, and measurement[5-7,13,14,66]. Coherent displacement is represented by the unitary operatorwhere gives the displacement amplitude and represent the annihilation and creation operators. It is one of the operations for which a closed form formula exists[66], given aswhere denotes the associated Laguerre polynomial function[72]. This relation only covers the lower triangular matrix; the rest of the matrix can be easily recovered from (6) usingThe formula (6) can be computed in multiple different ways with varying numerical accuracy impacted by the simplifications made in the expression and the order of their evaluation. When implemented exactly as it stands in (6), it is plagued by the limitations of FP arithmetic. Its first term underflows for comparatively large m, while the second term overflows for and large enough difference . When both the numerical underflow and the overflow coincide, the ill-defined expression is evaluated, resulting in error. We discuss the circumstances in detail in Sect. S1 of the Supplementary material and establish a set of acceptable combinations of the m, n and parameters such that formula (6) is always well defined.

We can utilize the recurrent formulae[53] or the plain matrix exponential[73,74] with a truncated argument instead of the closed form formula (6). While we can not ascertain their accuracy without a priori knowledge of the ideal operator, we can easily determine whether the generated matrices G are outright incorrect by checking the normalisationof displaced truncated Fock states . It corresponds to the sum of squared absolute values of elements in the jth column of the truncated displacement matrix or its approximation employing the matrix exponential with truncated argument where we set .

In Fig. 2 we show the normalisation (8) for  constructed using the closed form formula (6), represented with a blue solid line, the recurrent formula[53] shown with a red dash-dotted line, and approximated with the matrix exponential (black dashed line). We utilize double precision[60,69] in the computation and try to avoid numerical issues plaguing the direct method (6) by keeping the working dimension sufficiently low. There are two regions of qualitatively distinct behaviour in the plot. The first region, spanning the first 40 Fock states, shows correct normalization for all three methods of construction. In the following region the normalisation dwindles for both the closed form and the recurrent formulae whilst the matrix exponential remains incorrectly normalized. It remains normalized only because the matrix exponential function, by definition, produces unitary matrices from anti-Hermitian arguments. Unitarity is not necessarily the desired outcome here since the goal is to obtain the correct matrix rather than the computed approximation .

Figure 2

Normalisation (8) of individual displaced truncated Fock states with . The displacement operator is constructed on 101 dimensional TFS using the closed form formula (blue solid), the recurrent formula (red dash-dotted), and approximated with the matrix exponential (black dashed line).

Let us explicitly discuss the issue at hand. The displacement operator (5) is unitary by definition. Columns of its matrix representation can be understood as coefficient vectors of displaced Fock states. In the infinite-dimensional case these states should be normalized, that is the vector 2–norm[75] of each column should satisfy . However, this will not generally hold in finite dimension where we can find a threshold state that, when displaced, will not be properly represented on the TFS. The states will suffer from non-negligible errors (3), making their normalization .

The plot in Fig. 2 reveals that when the matrix is constructed via (6), the higher states are correctly denormalized. Conversely the matrix exponential produces incorrectly normalized states. In this context such behavior can be considered a manifestation of truncation errors.

The normalisation of the recurrently computed matrix starts to rise exponentially somewhere around  due to accumulation of rounding errors. This behavior depends on the chosen and the breakdown is more prominent when is large. Here the displacement  was chosen to emphasize this effect. For instance, when , a similar exponential breakdown appears for  instead.

Truncated approximate matrix exponential (TAME)

So far we have seen that, when it comes to numerically generating truncated representations of unitary operators, both direct calculation and the recurrent formulae have fundamental issues leading to significant rounding errors or numerically invalid expressions. The matrix exponential function avoids these issues mostly at the cost of truncation errors and their subsequent amplification. However, the observations in Fig. 2 also suggest that these errors tend to be significant only in higher regions of said matrices.

This opens up a new possibility of approximating the exponential operators. We can use the matrix exponential on a sufficiently higher dimension and only then truncate the result to the required , thus avoiding the erroneous areas, while, at the same time, keeping the computational dimension low enough to avoid needlessly increasing the time of computation. We call this approach truncated approximate matrix exponential (TAME). Consider the approximation of the truncated displacement operator, , constructed in such a way,Here represents the initial working dimension and the final dimension of the target TFS. Following (5) we set .

In Fig. 3 we compare constructed using the closed form formula (6) and approximated with TAME. We chose the dimension and the displacement magnitude to accommodate the limits established in Sect. S1 of the Supplementary material. The secondary dimension was chosen high enough to suppress the effects of truncation errors. The plot suggests that our method produces results equal to the closed form formula in terms of the normalisation (8). Further comparison of individual matrix elements reveals that, on average, the approximate matrix matches (6) up to 14 decimal places with the worst difference matching only up to 11 decimal places.

Figure 3

Normalisation of displaced Fock states . The matrix was constructed using the closed form formula (black bullets) and approximated with TAME (red, solid) where we set and . In both plots the dimension  for TAME was determined via Algorithm 1.

What remains to be determined is the proper choice, or rather the methodology of choosing a sufficiently large working dimension given the target dimension . In the subsequent paragraphs we are going to show that it is practical to set the dimension as small as possible. The de facto standard scaling and squaring matrix exponentiation algorithm[73,74] relies on matrix multiplication with the actual number of matrix products depending on the binary logarithm of the 1–norm[75] of the exponentiated matrix.

The 1–norm[75] of the argument inside the matrix exponential within (9) readswhere the final approximation holds asymptotically. Therefore the asymptotic computational complexity of the matrix exponential in (9) scales as  in terms of matrix products. The complexity of each matrix multiplication, specified in terms of FP operations, depends on the algorithm it utilizes. A naive textbook implementation scales as poorly as , whereas the more sophisticated Strassen algorithm[76] scales approximately as . Consequently the computational complexity of (9) scales as under optimal conditions. It is therefore imperative to keep the dimension as low as possible.

We propose a simple iterative algorithm for finding optimal . Suppose a sufficiently sized matrix is correct on some region spanning where . Suppose the matrix exponential () algorithm is also consistent: for a differently sized matrix with dimension the computed matrix exponential is correct on a region of at least the same size. Given these assumptions, which are upheld by the standard implementation[73,74], we introduce the Algorithm 1 as follows. First we take the desired dimension of the correct region and set an equality tolerance for small numbers: our condition with proclaims two numbers identical if they match up to their twelfth decimal place. Then we search for a pair of larger matrices such that their regions match. The search process is significantly simplified by taking the dimension of the second larger matrix to be constantly shifted from the first larger matrix. To improve its speed we always recycle one of the matrices in the next iteration instead of recalculating it every time. The depth of the search is specified by the factor h. In our experience the dimension is found somewhere well below  in the case of displacement, hence we set the depth h above that. Once the search algorithm finishes successfully, we obtain .

Verification of approximated matrices

In general, we can not verify the matrix (9) constructed via TAME simply by comparing its elements against some exact solution for the obvious reason: if we knew the exact solution we would not be in this situation in the first place.

We have used normalisation (8), or more precisely the implied necessary condition of unitarity , to detect outright incorrect matrices in Fig. 2, but alas, necessary conditions alone can not be used to prove the matrix correct. In Fig. 2 we determined that employing the recurrent formula[53] in construction of was ill-advised due to accumulation and consequent amplification of rounding errors over the course of the computation. While we can not safely use the recurrent formula to construct an arbitrary truncated displacement matrix, we can use it to determine whether a candidate matrix, for example one constructed via TAME (9), possesses appropriate structure as the formulae define relations between neighboring matrix elements.

We can repurpose the relations Eq. (56–58) from[53] to construct an error matrixfor a given candidate matrix G. The rounding errors are not amplified in computation of the error matrix as there is no recursion. Its elements give the difference between the actual elements of the candidate matrix and the values they should have been based on their neighbors, and , and the structural constraints given in[53].

In Fig. 4 we compare the decadic logarithm of the difference for approximated using (a) TAME (, ) and (b) the plain matrix exponential (). In each plot we display the row-wise value with blue line. The surrounding light-blue area stretches one standard deviation from the mean. The maximal difference within each column is represented by the red line. Finally the dashed black horizontal line (at ) roughly corresponds to the double precision unit round-off[60].

Figure 4

Verification of the matrix approximated using (a) TAME ( and ) and (b) plain matrix exponential (). Blue lines mark the row-wise values, light-blue region stretches a standard deviation away from the mean. The maximal difference within each column is represented by the red line. The dashed horizontal line corresponds to the unit round-off in double precision floating point number representation. (a) The matrix is structurally correct. The average differences are negligible, their values falling below the unit round-off. The maximal differences match up to 11 decimal places. (b) The matrix maintains correct structure in its first third. The truncation errors manifest in the rest of the matrix as an exponential explosion in the maximal difference (around the 100th column) and a steady rise in the mean value.

In Fig. 4 (a) the matrix is structurally correct, with the maximal difference still matching up to 11 decimal places. On average the differences fall below the unit round-off, essentially making the errors negligible. In Fig. 4 (b) the matrix constructed using the plain matrix exponential maintains the correct structure in the first third of its columns, however, the truncation errors begin to manifest at that point. This can be observed as an exponential explosion in the maximal difference (around the 75th column) and a steady rise in the mean value. We saw a similar manifestation of truncation errors in Fig. 2 where the columns incorrectly retained their normalization as if the truncated matrix remained unitary.

Numerical simulation of the preparation circuit

The CV nature of the preparation scheme in Fig. 1, described with relations (1) and (2), makes its exact numerical simulation not only impractical, but outright impossible. We can, however, perform an approximate numerical simulation of the formulae on a TFS. We have already proposed TAME as the method for approximating the truncated displacement operator. We have yet to ascertain a key ingredient of the simulation. We must determine the optimal dimension of the TFS, which should be large enough to support all the quantum states occurring in the simulation.

Following the Fig. 1, we begin with the TMSV state. One of its modes is attenuated by the channel. This only reduces its energy and, as a consequence, the required support shrinks in size. We can therefore safely disregard the attenuating channel and simplify the expression for the marginal state (1) into . We then require that both the initial and the displaced TMSV states are faithfully approximated on the dimensional TFS for all the possible values of and . By taking the largest displacement and squeezing rate considered in the simulation, we can iteratively determine as the least dimension such that the cutoff error (3) falls below some threshold . This condition readsfor the displaced TMSV state. While the coefficients of the TMSV state are determined trivially, the matrix elements of the displacement operator can not be, in general, obtained analytically with (6) and we must employ alternate means such as TAME.

In the simulation we consider , corresponding to roughly squeezing, and , hence we set while searching for . Once the dimension is found, we determine its respective using the Algorithm 1. With the thresholds we get and for . Note that for this particular , the TAME matrices constructed on and dimensional TFS are identical in double precision FP arithmetic.

In the following sections we use the numerical methodology we developed to determine the benefits of using PNRD, APD, and APD cascades in a pair of applications of the preparation circuit. First we discuss preparation of non-linearly squeezed states (Section “Engineering non-linearly squeezed states”) and then follow with construction of well defined non-classical superpositions of Fock states (Section “Preparation of high fidelity qubit in Fock basis”). In both applications the figures of merit are functions depending on the resulting density matrix (1) and the associated probability of success (2). We approach the analysis with a rudimentary grid based exploratory strategy for optimization. We divide the region into equidistant grid of points and evaluate the numerically approximated relations (1) and (2) for each point and each experimental scenario defined by the overall efficiency of the setup and expected measurement outcome respective to the POVM elements. These entailrepresenting the click of the ideal APD and the first six PNRD outcomes relevant in our preparation scheme, as well as PNRD approximations employing APD cascadeswhere . The POVM elements represent outcomes where exactly n detectors click within APD cascade comprising M detectors[52]. The individual probabilities readThis way we procure an assortment of probabilities P(i, j) and density matrices corresponding to the and sequences. We then utilize these values in objective and constraint functions, that will be discussed in detail in the following sections, to analyze the performance of the preparation scheme in particular applications and its response to different experimental configurations.

The numerical simulation and the analysis of its results was implemented using a number of open source software libraries[77-83] in Python.

Engineering non-linearly squeezed states

Nonlinear squeezing was originally introduced[51] as a measure quantifying the quality of approximate cubic states suitable for optical measurement-induced quantum gates[37,51]. It has been shown to apply to higher ordered phase squeezing gates as well[37] and was recently discussed in detail[67]. The ideal cubic operation facilitates unitary evolution with interaction Hamiltonian proportional to . When approximatively implemented in the measurement induced fashion[51] its action in the Heisenberg picture can be represented by the operator transformationwhere the operators correspond to the signal state and to some ancillary mode. The first terms of both relations correspond to the ideal cubic interaction . The additional term () represents the nonlinear quadrature of the ancillary mode and embodies the undesirable noisy contribution. It can be suppressed by choosing an appropriate ancillary state with the right structure. Effects of this contribution, or more precisely its variance and mean, vanish for the ideal cubic state. In general, neither the variance nor the mean vanish for physical approximations of the ideal cubic state. Good approximations, however, minimize their values and consequently the variance of this contribution may be used to quantify the quality of these approximate cubic states.

The preparation scheme presented in Fig. 1 can be utilized for production of quantum states approximating the ideal cubic state. We have discussed the methodology of the simulation in detail in Section Numerical simulation of the preparation circuit. In essence we search for optimal values of squeezing and displacement that lead to high quality cubic state approximations while maximizing the probability of successful preparation. The optimization is performed for various experimental scenarios involving different detectors and taking a range of overall losses into account.

To measure the approximation quality we adapted the nonlinear quadrature and the concept of nonlinear squeezing discussed in[51] to fit our simulation. We employ the nonlinear quadrature and use its variance to measure the nonlinear squeezing of arbitrary states . Potential effects of Gaussian squeezing[51,67] on are eliminated by minimizing over the parameter . Consequently we base our analysis on the minimized quantity normalized with respect to the minimal variance achievable by Gaussian states [67].

The numerical simulation yields density matrices along with the P(i, j) probabilities of success corresponding to different experimental parameters. We then compute the individual moments required in the calculation of from the elements of density matrices . We avoid the matrix representation of the operators in the computation to avert truncation errors and employ closed form formulae instead. The minimization with respect to within is solved analytically.

We thus obtain M(i, j) values for their respective matrices and P(i, j) probabilities. We then divide the dataset corresponding to each experimental scenario into bins based on values of the variance M(i, j) and find the maximal attainable probability P(i, j) within each bin.

In Fig. 5 we present a comparison of the attainable variance as a function of success probability P. We examine different detection outcomes, in particular the PNRD projection onto (red line) and its three approximations realized through an APD cascade[52] where three APD detectors out of four (dashed magenta), five (magenta) and ten (blue) click. Their respective POVM elements , and are given by the relation (14). We consider a single APD detector (black line) as well. The plots show (a) , (b) , and (c) transmission efficiency . The results are normalized with respect to the minimal variance achievable by Gaussian states. The optimal cubic state approximations[51] constructed on v dimensional TFS are marked with dashed horizontal lines. These states were found by searching for pure states spanning the first v Fock states that would minimize the variance of the non-linear quadrature[51]. Their inclusion makes it possible for qualitative comparison with the states produced by our scheme.

Figure 5

A comparison of attainable variance as a function of success probability. The variance is normalized with respect to the minimal variance achievable by Gaussian states. We use the same vertical and horizontal axes in the plots to show the contrast between the almost ideal (a) and lossy (b, c) scenarios with , and transmission efficiencies. Horizontal dashed lines are used to mark the optimal cubic state approximations constructed on low-dimensional TFS. We encode the information about the POVM elements as follows: APD click with solid black line, PNRD projection onto with solid red, APD cascades comprising four (, dashed magenta), five (, magenta) and ten (, blue) detectors where three detectors click. Overall, utilizing the PNRD (solid red) produces states with lowest non-linear variance, therefore producing comparatively better approximations of the cubic state. In both (b) and (c) a single APD outperforms the APD cascades comprising five and four detectors for probabilities greater than . In this regime the cascade comprising ten detectors still offers advantage over single APD. In (c) a single APD outperforms APD cascades comprising either four, five or ten detectors for success probabilities larger than roughly .

In general using PNRD yields the best results. In the idealized scenario with efficiency the PNRD projecting onto  approaches the variance of the optimal non-linearly squeezed state . It also attains the best values consistently across the considered transmission efficiencies, therefore producing comparatively better approximations of the cubic state than either the APD cascades or a single APD. In the and regimes, the APD cascade comprising ten detectors promises better performance than a single APD or any other cascade configuration for that matter. In the low-efficiency mode () we can see that a single APD outperforms APD cascades for probabilities of success greater than . This can be attributed to the imperfections inherent to APD cascades[52]. Their flaws become emphasized with increased loss, rendering a single APD to be the better choice.

In conclusion, unless a PNRD capable of distinguishing at least three photons is available, it is advantageous to use a single APD in any practical scenario with non-ideal transmission efficiency as long as success probabilities larger than approximately are desired. The advantage of a single APD can be offset by using an exorbitant number of detectors within APD cascade.

Preparation of high fidelity qubit in Fock basis

The qubit-like superposition represents one of the simplest non-Gaussian quantum states of light. It serves an important role in quantum information processing[84] and is one of the resources available in contemporary experimental quantum optics[85-87]. As such it is has been employed in experimental demonstrations of various theoretical concepts including witnessing of non-Gaussianity[88-91] and hybrid entanglement[44,92] in quantum communication.

This family of quantum states can be produced with the preparation scheme we have previously introduced in Fig. 1. We can search for the optimal squeezing and displacement parameters to obtain a given target state with sufficient fidelity and maximal performance in terms of success probability.

We compute the corresponding values for the P(i, j) probabilities and density matrices obtained in the simulation described in detail in Section Numerical simulation of the preparation circuit. We then divide the dataset for each experimental scenario into bins comprising subsets of data satisfying where specifies a moving fidelity threshold. The maximal attainable probability P(i, j) is then found for each subset.

In Fig. 6 we demonstrate the relative improvement in probability of successfully engineering states by employing different detectors instead of a single APD. We consider a pair of target states, and , both evaluated for and transmission efficiencies. These target states were chosen to probe the improvement for unbalanced superpositions biased either towards or states. In the plot we show the result obtained for projection onto (red line) realized by PNRD and the results obtained with its approximations realized through APD cascades where a single detector out of ten (, magenta), five (, blue) and two (, black) clicks. The POVM elements of the cascades were defined in (14). The figure of merit is defined as with respective to individual detection outcomes.

Figure 6

Benchmarking the performance of PNRD and its approximations using APD cascades relative to a single APD detector in preparation of particular superpositions parametrized with . The PNRD projection onto is represented by red line, whereas the magenta line corresponds to APD cascade comprising ten detectors where a single detector clicks (), blue line to cascade of five detectors () and black line depicts the case with two detectors (). The plots demonstrate preparation of two distinct states while considering different transmission efficiencies. In (a) and (c) we aim to prepare . In (b) and (d) we target . In plots (a) and (b) we consider transmission efficiency, while in (c) and (d) we consider mere . The horizontal dashed line marks a twofold improvement in each plot. (a) In the high-efficiency regime we obtain roughly tenfold improvement in the high-fidelity preparation of the state. The PNRD approximations using more than two detectors offer a significant improvement as well. (b) While the advantage of PNRD is reduced when targeting the state biased towards , it still offers roughly four times better performance. (c) The PNRD detector and its approximations offer 2-3x higher success probability even in the lower-efficiency scenario. (d) The PNRD detector and its approximations offer roughly twofold improvement.

In general, conditioning on the PNRD detection outcome yields the best results. In the high-fidelity regime with efficiency the relative improvement is roughly tenfold for and roughly four times better for . The APD cascades comprising ten and five detectors follow. The relative lead of the PNRD diminishes in the efficiency regime. Its advantage also dwindles when we consider target states closer to , such as the state. While the cascade comprising two detectors falls short in every case, it still outperforms a single APD, albeit not by a lot.

Conclusions

We have analyzed the numerical accuracy of several currently available methods[53,66] used in construction of the truncated coherent displacement operator, an essential ingredient of state preparation in quantum optics[14,30,32,36] and many protocols used in quantum information processing[5-7,13,14,66]. We have proposed an alternative approach promising a better accuracy. Our method is based on the standard matrix exponential[73,74] with truncated argument. We compute the matrix exponential on a higher-dimensional space and truncate the resulting matrix to the target dimension, thus stripping erroneous matrix elements away from the truncated displacement operator. To avoid negatively impacting computational performance, the higher dimension should be ideally kept as low as possible. To this end we provide an off-line search algorithm that can be used to determine its optimal value. To ascertain the accuracy of the resulting matrix we complement the construction method with a verification strategy based on the recurrent formulae discussed in[53].

We have used our construction method for analysis of non-Gaussian state preparation scheme based on suitable manipulation of a two mode squeezed vacuum with subsequent photon counting measurement[14,39,63,64] in the context of engineering non-linearly squeezed cubic states[14,38,51] for measurement induced cubic gates[14,37,51] and construction of qubit-like superpositions between vacuum and single photon states. The latter application can be verified experimentally with currently available technology. We have compared the effects of different detection mechanisms, including APD, PNRD and its approximations using APD cascades[52] with varying number of APD detectors, to determine practical approach towards state preparation. In our analysis we have optimized the free parameters of the prepraration scheme, the initial squeezing and the displacement, to attain optimal results in both applications. This analysis also provides additional metric which can be used to quantify the quality of APD cascades. We have found that in practical applications when PNRD is not available, using a single APD to engineer non-linearly squeezed states offers better performance compared to employing APD cascades comprising small number of detectors. We attribute this counter-intuitive result to the imperfections inherent to APD cascades[52] which are exaggerated with increased loss; these flaws became significant for overall loss. The primary cause of this behaviour lies within the employed avalanche detectors as a single click may be triggered by multiple photons. While this is a critical issue for multi-photon state engineering, it is not as significant for single-photon states. We have determined that using APD cascade, even if one comprising only a pair of APD detectors, improves upon using a single APD in preparation of high-fidelity non-Gaussian qubit-like superpositions.

Our circuit variant can be extended to utilize multiple displacements and detectors. Similarly the proposed method for numerical construction of truncated unitary operators is not limited to displacement only can be applied to, for example, squeezing or cubic phase-shift operators. Furthermore, the method could be employed in preparation of a wider variety of quantum states with practical applications, such as GKP states[14].

Supplementary Information.