Quantum process tomography with coherent states

We develop an enhanced technique for characterizing quantum optical processes based on probing unknown quantum processes only with coherent states. Our method substantially improves the original proposal [M. Lobino et al., Science 322, 563 (2008)], which uses a filtered Glauber-Sudarshan decomposition to determine the effect of the process on an arbitrary state. We introduce a new relation between the action of a general quantum process on coherent state inputs and its action on an arbitrary quantum state. This relation eliminates the need to invoke the Glauber-Sudarshan representation for states; hence it dramatically simplifies the task of process identification and removes a potential source of error. The new relation also enables straightforward extensions of the method to multi-mode and non-trace-preserving processes. We illustrate our formalism with several examples, in which we derive analytic representations of several fundamental quantum optical processes in the Fock basis. In particular, we introduce photon-number cutoff as a reasonable physical resource limitation and address resource vs accuracy trade-off in practical applications. We show that the accuracy of process estimation scales inversely with the square root of photon-number cutoff.


Introduction
Assembling a complex quantum optical information processor requires precise knowledge of the properties of each of its components, i.e., the ability to predict the effect of the components on an arbitrary input state. This gives rise to a quantum version of the famous "black box problem", which is addressed by means of quantum process tomography (QPT) [1,2,3]. In QPT, a set of probe states is sent into the black box (here an unknown completely-positive, linear quantum process E over the set of bounded operators B(H) on a Hilbert space H) and the output states are measured. From the effect of the process on the probe states it is possible to predict its effect on any other state within the same Hilbert space.
QPT exploits linearity of quantum process over its density operators. If the effect of the process E(ρ i ) is known for a set of density operators {ρ i }, its effect on any linear combination ρ = β i ρ i equals E(ρ) = β i E(ρ i ). Therefore, if {ρ i } forms a spanning set within the space L(H) of linear operators over a particular Hilbert space H, knowledge of {E(ρ i )} is sufficient to extract complete information about the quantum process.
However, practical implementations of QPT become demanding especially for systems with large Hilbert spaces. For dim(H) = d, dim L(H) = d 2 , which implies that at least d 2 unknown operators {E(ρ i )}, each with d 2 unknown parameters, must be estimated. This procedure requires preparation of at least {ρ i } d 2 i=1 states, subjecting each to the unknown process E, and determining each element of {E(ρ i )} d 2 i=1 through measurement (each with d 2 unknown elements), thereby inferring an overall number of d 4 parameters. Furthermore, in order to build up sufficient statistics for reliable estimates of the output states, each measurement should be performed many times on multiple copies of the inputs. Thus, a large number of realizations and measurements are required for complete tomography of E.
An additional complication, especially for QPT of quantum optical processes, is associated with preparation of the probe states. Typical optical QPT implementations deal with systems consisting of one or more dual-rail qubits [4,5,6], which implies that the probe states are highly nonclassical, hence difficult to generate.
These difficulties have been partially alleviated in the recently proposed scheme of "coherent-state quantum process tomography" (csQPT) [7]. This scheme is based on the observation that the density operator ρ of a generic quantum state of every electromagnetic mode can be expressed in the Glauber-Sudarshan representation [8,9], where P ρ (α) is a quasi-probability distribution referred to as the quantum state's "P -function" and integrated over the entire complex plane [10]. Linearity hence implies that measuring i.e., determining the effect of the unknown process on all coherent states enables a prediction of its effect upon any generic state ρ according to The implementation of csQPT is advantageous because (i) coherent states are readily available from lasers, (ii) coherent states of different amplitudes and phases can be produced without changing the layout of the experimental apparatus, and (iii) output-state characterization can be performed using optical homodyne tomography [11], which obviates the need for postselection and provides full information about the process in question.
Moreover, csQPT has been tested experimentally on simple single-mode processes, such as the identity, attenuation, and phase shift operations [7]. Furthermore, csQPT has been used to characterize quantum memory for light based on electromagnetically-induced transparency [12].
An apparent obstacle to csQPT, however, is that the P function for many nonclassical optical states exists only in terms of a highly singular generalized function [13,14]. A remedy therefor is provided by Klauder's theorem [15], which states that any trace-class operator ρ can be approximated, to arbitrary accuracy, by a bounded operator ρ L ∈ B(H) whose Glauber-Sudarshan function P L is in the Schwartz class [16], so integration (3) can be performed. The Klauder approximation can be constructed by low-pass filtering of the P function, i.e., by multiplying its Fourier transform with an appropriate regularizing function equal to 1 over a square domain of size L × L and rapidly dropping to zero outside this domain. Ref. [7] employs this method to implement csQPT.
Practical implementation of Klauder's procedure is however complicated, because it requires finding the characteristic function of the input state and subsequently its regularized P function. This function features high-frequency, high-amplitude oscillations that limit the precision in calculating the output state (3). Furthermore, Klauder's approximation is ambiguous in the choice of the particular filtering function as well as the cutoff parameter L [7].
Here we improve csQPT to overcome the above problems. Specifically, we develop a new method for csQPT that eliminates the explicit use of the Glauber-Sudarshan representation and thus removes the inherent ambiguity associated with employing Klauder's approximation for csQPT. In Sec. 2.1, we obtain an expression for the process tensor in the Fock (photon number) basis that can be directly calculated from the experimental data. Using this tensor, the process output for an arbitrary input can be calculated by simple matrix multiplication rather than requiring integration and high-frequency cut-offs. In this way, transformations between the Fock and Glauber-Sudarshan representations, which were necessary in Ref. [7], can be sidestepped. Using our new approach, we easily extend csQPT from its restrictive single-mode applicability to multi-mode processes and even to non-tracepreserving conditional processes. These extensions are particularly relevant for quantum information processing circuits, whose basic components are inherently multi-mode and conditional [17].
Process tomography is successful if, for every input state, the estimate for the process output closely approximates the actual process output state, and the worst-case error of this estimate, given by a distance between the actual and estimated process outputs, is less than a given tolerance. For states over infinite-dimensional Hilbert spaces, this concept of error is however not meaningful because the finiteness of sampling implies that the process is necessarily under-sampled, hence cannot be determined with bounded error. Instead we could consider the process estimation restricted to a finite-dimensional subset of B(H). This version of process tomography can always be successful with a sufficiently large amount of sampling.
Of particular practical interest is the subspace B(H) defined by an energy cut-off, i.e., estimating the process without accessing any information about its high-energy behavior. This restriction is naturally consistent with our choice to work in the Fock basis, because then the resulting process tensor is of finite size and with many practical settings (e.g. quantuminformation processing with photonic qubits). In Sec. 2.2, we provide process error estimates for several input state subsets that extend beyond B(H).
Many interesting processes are phase symmetric; that is, an optical phase shift of the input state results in the same phase shift of the output. This property dramatically simplifies the experiment because one needs to collect data only for coherent states whose amplitudes lie on the real axis rather than the entire complex plane. This prompts us to discuss, in Sec. 3, how to obtain the process tensor for phase-symmetric processes, which we test on the experimental data from Ref. [12]. Next, in Sec. 4, we illustrate our method by analytically deriving the superoperators for certain fundamental quantum optical processes using the Fock basis. The paper is concluded in Sec. 5 and is supplemented with two appendices.

Formalism: determining the quantum process matrix
We study general quantum optical processes E acting on quantum states of light and begin with the simplest case for which only a single electromagnetic field mode is involved. An arbitrary quantum state ρ can be expressed in the Fock representation as Subjecting this state to an unknown process E, and imposing linearity, yields where E mn jk := j|E(|m n|)|k is a rank-4 tensor, hereafter referred to as the "process tensor" (superoperator). Thus, by expressing input and output states in the Fock basis, a quantum process can be uniquely represented and characterized by its rank-4 tensor, which relates the matrix elements of the output and input states according to Below we show how to estimate process tensor elements E(|m n|) for m, n over a finite domain. Because α| (|m n|) |α = e −|α| 2 α nᾱm √ m!n! (8) is in the Schwartz class, the Glauber-Sudarshan P representation is guaranteed to exist for any operator |m n| (m, n ∈ N 0 ) [18]. The P function is for ∂ m α := ∂ m /∂α m and α and its complex conjugateᾱ treated as independent variables, and δ 2 (α) ≡ δ Re(α) δ Im(α) . By inserting representation (9) into Eq. (6), and exploiting linearity of the process, we obtain the process tensor This expression can be simplified by using Eq. (10) and performing integration by parts: Thus we have eliminated the need to make use of the Glauber-Sudarshan representation for quantum states. The process tensor is found by taking partial derivatives (with respect to α andᾱ) of the matrix elements of ̺ E (α), which are estimated from experimental data and evaluated at α = 0.
The mathematical procedure defined by Eq. (12) is simpler and computationally faster (see Sec. 3) than employing Eq. (11) with a regularized version of P L,mn (α) replacing the tempered distribution P mn (α) described in Refs. [7,12]. Equation (12) has been used to determine the fidelity of quantum teleportation of a single-rail optical qubit based on measurements performed on coherent states (see supplementary material in Ref. [19]).
Generalization to the multi-mode case is straightforward. In the M-mode case, let us introduce the notation |n := |n 1 , n 2 , . . . , n M (with n ∈ N M 0 ) for multi-mode Fock states and |α := |α 1 , α 2 , . . . , α M (with α ∈ C M ) for multi-mode coherent states. Then the matrix elements of the output and input states with respect to the Fock basis are related to one another by the rank-4 M tensor where Similarly to the single-mode case, we employ the Glauber-Sudarshan P representation for the multi-mode operator |m n|, with the overall P function being a product of the P functions for the constituent modes: Multiple integration by parts yields where Equations (12) and (16) complete our coherent-state tomography formalism and show that coherent states provide a complete set of probe states for characterizing quantum optical processes, insofar as the expression for ̺ E (α) completely determines the process tensor. The above formalism is not restricted to trace-preserving quantum processes. Indeed, trace preservation was not required in the derivation of our results. Thus, our method is applicable to all quantum optical processes that are mathematically described by completelypositive maps, but may be trace-preserving, trace-reducing or even trace-increasing. Tracenonpreserving quantum processes are either conditional processes or part of a larger process E = E 1 + E 2 , which is trace-preserving as a whole, but whose components E 1 and E 2 may increase or decrease the trace, respectively. A conditional process is a process that is conditioned on a certain probabilistic event; it may be heralded if the event is observed. One of the most notable examples of such a process is a probabilistic conditional-NOT gate (CNOT), which forms the basis for the Knill-Laflamme-Milburn linear-optical quantum computing scheme [17]. Other examples are photon-addition and photon-subtraction processes, whose superoperators are derived in Sec. 4.
In experimental csQPT, states ̺ E (α) are determined using homodyne tomography [11]. It is important to remember, however, that this procedure reconstructs a density matrix normalized to unity trace: When measuring non-trace-preserving processes, one must recover the trace information contained in ̺ E (α). This is done by measuring the probability p α (E) = Tr [̺ E (α)] of the process heralding event for all α's for which the measurements are performed. The state to be used in Eqs. (12) and (16) An interesting feature of Eqs. (12) and (16) is that complete information about a quantum optical process is contained in its action on an infinitesimally small compact set of all probe coherent states in the immediate vicinity of the vacuum state. From a mathematical point of view, this feature can be understood by realizing that, for any j, k ∈ N 0 , the matrix element j|̺ E (α)|k is an entire function (see Appendix A), i.e., a complex-valued function in the variables α,ᾱ that is holomorphic over the whole complex plane, and so is its product with the exponential e |α| 2 . Hence, each term e |α| 2 j|̺ E (α)|k is infinitely differentiable over the whole complex plane and is identical to its Taylor series expansion in any point of C. Moreover, Eq. (12) implies that the process tensor is determined by the corresponding Taylor coefficients at α = 0. The same conclusion applies to the multi-mode case, in which we deal with entire functions on C M .

Energy cutoff and estimation of the error of approximation
As discussed in Sec. 1, the incompleteness of the information acquired in the experiment is accommodated in csQPT by evaluating the process tensor over a restricted finite-dimensional subspace H of the Hilbert space H with a fixed maximum number N of photons. The incurred expense is that, through this reduced tomography, only approximate information about the process will be inferred: for a given input state ρ, the predicted output is not E(ρ), but rather is the trace-normalized projection of ρ onto B( H) and is the predicted output of the reconstructed process for input stateρ. In Eqs. (18) and (19), Π is the projection operator onto H. If the input state ρ is outside B( H), the process output estimation error E(ρ) −Ẽ(ρ ) 1 (where ρ 1 = Tr ρ † ρ is the trace norm) is generally unbounded. However, it is possible to bound the error for certain practically important classes of input states and processes.
For example, all linear-optical processes involving only linear-optical elements (interferometers, attenuators, conditional measurements) do not generate additional photons, and thus map B( H) onto itself, soẼ(ρ) = E(ρ). For such processes, the error for a particular input ρ can be estimated according to E(ρ) −E( ρ) 1 ≤ E ρ− ρ 1 , with the superoperator norm defined as E ≡ sup{ E(B) 1 :B ∈ B(H) , B 1 ≤ 1} [20]. If the process is known to be trace-nonincreasing, we have E ≤ 1 [21] so the error is bounded from above by Note that the above result is not sufficient for evaluating the error for a general process, because this error is given by the deviation of E(ρ) from E( ρ) rather than from E( ρ) [Fig. 1].
A further error bound can be obtained for the class of trace-preserving processes that do not increase the mean energy, acting on a set of input states whose mean energy does not exceed a certain value [22]. We illustrate this for a single optical modeâ with frequency ω and Hamilton operatorĤ = ω(â †â + 1/2) whose eigenvalues are denoted by h n = (n + 1/2)ω. Suppose that the quantum states ρ of interest satisfy Tr[ρĤ] ≤ U. According to Ref. [22], if we choose the cutoff dimension dim( H) = N + 1 such that U/h N +1 ≤ γ for some (small) γ > 0, the reconstructed process output errors are bounded from above as where Conversely, if we want to achieve a certain upper bound ǫ on the error of approximation (which corresponds to a lower bound on the desired accuracy of the process characterization), we first solve Eq. (22) for γ = γ(ǫ), and then find the minimum N γ ∈ N such that U/h Nγ +1 ≤ γ. Any cutoff dimension N + 1 > N γ is then sufficient for our purpose. For γ ≪ 1, ǫ ≈ 2 √ γ, which yields This implies that the error of approximation scales as 1/ √ N with the cutoff dimension N + 1. For example, in order to achieve a 10% error in Eq. (21), we need ǫ = 0.05 and thus γ ≈ 0.0006. For the input mean energy bound corresponding to one photon (U = 3/2ω), the required cutoff is N ≈ U/γ ≈ 250. This calculation shows that the above error estimate is extremely conservative.

Phase-invariant processes
Many practically relevant processes, including the single-mode processes studied in Sec. 4, exhibit phase invariance. If two input states are identical up to a shift by an optical phase φ, the process outputs for these states differ by the same phase shift: For such processes, it is convenient to express the probe coherent states in polar coordinates: |α = re iθ = e inθ |r . Specifically, in these coordinates, we have [9] P mn (r, θ) = √ m!n! (m + n)! e r 2 +iθ(n−m) (−1) m+n d m+n dr m+n δ(r), and accordingly Hence Eq. (24) can be expressed as and the superoperator E [Eq. (26)] has the following explicit representation: In experimental tomography of phase-invariant processes [7,12], it is sufficient to measure the process output for a discrete set of coherent states {|r i } on the real axis of the phase space. The matrix elements of the output states can then be interpolated as polynomial functions where Q is the degree of the polynomial (which depends on the dimension of the truncated Hilbert space) and C l (j, k) are its coefficients. Furthermore, from Eq. (A.3) together with Eq. (28), it follows that, for phase-symmetric processes, when j − k is even or odd, j| E(|r r|) |k and its analytic extension to negative values of r are even or odd functions of r, respectively. By taking into account the symmetric or antisymmetric property of this function, we have additional information to be used in the interpolation procedure; the constructed polynomial has to contain only even or odd powers of r, respectively. In this way the precision of process estimation from the experimental data is substantially increased.
With the knowledge of the coefficients C l (j, k), Eq. (28) is further simplified to: The last result is significant in that one can obtain the process tensor directly from the experimentally reconstructed output states through simple summation. Moreover, if the dimension of the truncated Hilbert space is d = N + 1, from Eq. (30) it follows that only terms of power l ≤ 2N of the interpolation polynomial (29) contribute to the process tensor.
We have tested this procedure on experimental data [12] and calculated the process tensor in a few microseconds, which is a dramatic improvement in comparison to several hours required for the original procedure [7,12].

Examples: superoperators of important quantum optical processes
In this section, we illustrate our new method by applying it to some fundamental quantum optical processes, whose effects on coherent states are known. Specifically, using Eqs. (12) or (16), we analytically derive corresponding superoperator tensors E mn jk in the Fock basis. The results are summarized in Table 1.

Identity
For the identity process (E id ), ̺ E id (α) = |α α|, the matrix elements of the output states are Inserting these elements into Eq. (12) yields E mn jk = δ mj δ nk , as expected.

Attenuation and lossy channel
For attenuation of light fields (E att ), the process's effect on single-mode coherent states is given by ̺ Eatt (α) = |ηα ηα|, where 0 ≤ η < 1. The matrix elements in the Fock basis are From Eq. (12), we obtain which depends explicitly on η.

Photon subtraction and addition
Photon subtraction is defined as a process that removes a single photon from the light field, whereas photon addition adds a single photon. Photon subtraction has been used by Ourjoumtsev et al. [23] to generate optical Schrödinger kittens (coherent superpositions of low-amplitude coherent states) from squeezed vacuum states for the purpose of quantum information processing. Single-photon-added coherent states can be regarded as the result of the most elementary amplification process of classical light fields by a single quantum of excitation; being intermediate between single-photon Fock states (fully quantum-mechanical) and coherent (classical) ones, these states have been demonstrated to be suited for the study of smooth transition between the particle-like and the wavelike behavior of light [24].
Here we discuss idealized single-mode photon subtraction and photon addition. Both processes are non-trace-preserving. For example, photon subtraction can be approximately realized [23] by a highly-transmissive beam splitter, whose reflected mode is directed to a detector and whose transmitted mode constitutes the output, respectively, as illustrated in Fig. 2a. Any click in a detector implies extraction of photon(s) from the input mode by the beam splitter. As the beam splitter has low reflectivity, here single-photon extraction events are more likely than multi-photon events. An approximate experimental realization of photon addition is illustrated in Fig. 2b. The input quantum state ρ enters the signal channel of a parametric down-conversion setup. Provided that detector dark counts are neglected, a photon detection in the idler mode heralds photon addition to the signal mode, which contains the output state of the process.

Schrödinger cat generation
The unitary evolution according toÛ Kerr (χ) ≡ exp −iχ â †â 2 for χ = π/2, if applied to coherent states, generates Schrödinger cat states (hereafter denoted as E cat ) [25,26] with matrix elements The superoperator tensor for this non-Gaussian unitary process obtained via Eq. (12) is Interestingly, this process does not change the total particle number of any input state.

Beam splitter
Now let us consider the beam splitter as an example of a two-mode process. The unitary beam splitter transformation is given by [27] where Θ is the parameter identifying how the beam splitter transmits or reflects beams. Specifically, its action on coherent state inputs |α 1 and |α 2 is given as with T ≡ cos(Θ/2) and R ≡ sin(−Θ/2) being the transmissivity and reflectivity, respectively. By knowing the effect of the process on two-mode coherent states, we can calculate the corresponding tensor using Eq. (16), which yields E m 1 m 2 n 1 n 2 as an explicit function of T and R.

Parametric down-conversion
Another two-mode process of interest is parametric down-conversion (PDC). In PDC, a crystal with an appreciably large second-order non-linearity is pumped by a laser field. Each of the pump photons can spontaneously decay into a pair of identical (degenerate PDC) or nonidentical photons (nondegenerate PDC). Here we consider a nondegenerate PDC process E PDC induced by the transformation [27] The effect of this unitary process on a two-mode coherent state is given by In Appendix B, we derive the process tensor in the Fock basis. The result can be expressed as: the hypergeometric function, (x) n := Γ(x + n)/Γ(x) the Pochhammer symbol and Γ(·) the Gamma function [28].

Conclusions
Coherent states are easily generated probe states for tomography of unknown quantumoptical processes. Here, we have presented a new, more efficient data processing technique for estimating a quantum process from similar experimental procedure of Ref. [7]. The original formulation was based on regularization and filtering of the Glauber-Sudarshan representations for quantum states, which are cumbersome to implement numerically. Furthermore, Ref. [7] introduces additional errors associated with regularization of the P function. In contrast, our new method to determine the process superoperator [Eq. (12) or Eq. (16)] is mathematically simpler, computationally faster and unique up to the choice of the energy cutoff. Moreover, we presented straightforward generalizations of coherent state quantum process tomography to multi-mode and non-trace-preserving conditional processes.
We have illustrated the new framework through several examples (summarized in Table 1). We have shown that it is straightforward to derive analytically exact and unique closed-form expressions for the superoperators for quantum optical processes whose effect on coherent states is known. For phase-invariant unknown processes, the formula to find the process tensor reduces to a simple summation of coefficients of a polynomial obtained from the experimentally reconstructed output states via interpolation.
An interesting consequence implied by our formulation [in particular, Eqs. (12) and (16)] is that complete information about a quantum optical process is entirely captured by its effect on a compact set of all coherent states |α in the immediate vicinity of the vacuum state. This is due to the entireness property of the image of processes on coherent states. It thus appears sufficient to perform tomography experiments only for a range of coherent states whose mean photon number is much smaller than that required for the method of Ref. [7] (see the suppl. material therein). However, coherent state quantum process tomography relies on the ability to approximately determine all the derivatives of a function which is obtained by interpolation from measured experimental data. Minimization of errors associated with this calculation imposes a lower bound on the phase space region over which the measurements need to be performed. For the time being, we have provided an evaluation of the error in the process estimation by introducing a truncation of the Fock space. For the class of processes respecting a certain energy constraint (which includes all processes that do not amplify the energy), we have determined (i) the cutoff dimension that is sufficient in order to achieve a certain degree of approximation accuracy, as well as (ii) the upper bound on the error of estimation for a given cutoff dimension.

Acknowledgments:
We acknowledge financial support by NSERC, iCORE, MITACS, QuantumWorks and General Dynamics Canada. AIL is a CIFAR Scholar, and BCS is a CIFAR Fellow. We would also like to thank Connor Kupchak for helpful discussions.
Appendix A. Proof that j| ̺ E (α) |k is an entire function According to Eq. (12), by knowing the complex-valued function j| ̺ E (α) |k (of the variable α) for any j and k, one can determine the process tensor E mn jk . Here we show that this function is an entire function so it can be represented as a power series that converges uniformly on any compact domain.
As a completely-positive quantum operation, E possesses a Kraus decomposition E(ρ) = L i=1Ki ρK † i , where L ≤ dim(H) 2 andK i are some Kraus operators on H (whose explicit form is not needed for our purpose). Hence we can rewrite the matrix elements of the output state as is the dual or adjoint map [29]. The complex-valued function α|Â |α (referred to as Husimi function ifÂ is a density operator)-whereÂ is any bounded operator on H-is an entire function of the two variables α andᾱ [13,30]. Hence the right hand side of Eq. (A.1) implies that the function j| ̺ E (α) |k is an entire function. By representing the coherent states in Eq. (A.1) in the Fock basis and using Eq. (6), we obtain which is a power series of the complex variables α andᾱ, hence convergent everywhere [13,30].
which can also be expressed in terms of a product of values of the hypergeometric function 2 F 1 , as given by Eq. (45).