Fourier Quantum Process Tomography

The characterization of a quantum device is a crucial step in the development of quantum experiments. This is accomplished via Quantum Process Tomography, which combines the outcomes of different projective measurements to deliver a possible reconstruction of the underlying process. The tomography is typically performed by processing an overcomplete set of measurements and extracting the process matrix from maximum-likelihood estimation. Here, we introduce a new technique, referred to as Fourier Quantum Process Tomography, which requires a reduced number of measurements, and benchmark its performance against the standard maximum-likelihood approach. Fourier Quantum Process Tomography is based on measuring probability distributions in two conjugate spaces for different state preparations and projections. Exploiting the concept of phase retrieval, our scheme achieves a complete and robust characterization of the setup by processing a near-minimal set of measurements. We experimentally test the technique on different space-dependent polarization transformations, reporting average fidelities higher than 90% and significant computational advantage.

In principle, one could extract the analytical relations between the operator parameters and the outcomes of suitable projective measurements [19].However, this proves to be often incompatible with realistic experimental noise, typically yielding nonphysical reconstructions.This inconvenience can be overcome by formulating the process tomography as an optimization problem, as first proposed for the tomography of quantum states [20].
In this framework, the most elementary scenario is the characterization of an SU(2) gate Û acting on a twolevel quantum system (qubit).Polarization of photons provides a natural way of encoding qubits, with Û implemented via one or multiple birefringent waveplates.Accordingly, the characterization of devices acting on light polarization can be accomplished via QPT [21].
Here, we address the more challenging scenario of characterizing optical SU (2) gates that are dependent on some d-dimensional degree of freedom, hereafter referred to as lattice.We introduce a new technique, named Fourier Quantum Process Tomography (FQPT), that allows retrieving all the parameters of the unknown transformation by processing only three sets of projective measurements collected in 2 conjugate planes.This method applies to SU(2 × d) transformations, which can be decomposed in a 2 × 2 block-diagonal form.* Correspondence email address: francesco.dicolandrea@uottawa.ca FQPT is validated experimentally on complex polarization transformations realized via liquid-crystal metasurfaces (LCMSs) patterned with high spatial frequencies [22], and its performance is compared with a standard maximum-likelihood (ML) approach.In this experiment, the measurements can be conveniently chosen to be performed in two conjugate planes, namely the near and far field, wherein the light distributions are directly connected via a Fourier transform.If the near field is associated with an intermediate plane or two intermediate planes are selected, then the Fourier transform is replaced by a paraxial Fresnel propagator [23].FQPT can also be implemented in other platforms.For instance, integrated photonic technologies [24] can support additional chips specifically implementing the Quantum Fourier Transform (QFT) algorithm [25][26][27].At the same time, SU(2) operations could be implemented either in the polarization [28][29][30] or path encoding.In the latter case, the waveguide array would simulate a composite lattice.Similar schemes have also been reported in several non-photonic platforms [31][32][33].

II. THEORY
A qubit rotation of an angle 2E around the axis n = (n 1 , n 2 , n 3 ), with 0 ≤ E < π and |n| = 1, is described by an SU(2) operator where σ 0 is the 2×2 identity matrix and σ = (σ 1 , σ 2 , σ 3 ) is the vector of the three Pauli matrices.The characterization of an optical SU(2) gate is typically performed by processing an overcomplete set of 16 projective measurements of the form where |a⟩ and |b⟩ are extracted from the three sets of states forming the Mutually Unbiased Bases (MUB) of arXiv:2312.13458v1[quant-ph] 20 Dec 2023 [1,8,10].The process tomography of the gate is then accomplished via an ML approach, i.e., by minimizing a cost function expressing the distance between the experimental outcomes I exp ab and the corresponding theoretical predictions I th ab [10,21]: This approach may become unfavourably timeconsuming and less accurate in the case of transformations acting on high-dimensional Hilbert spaces.
Here, we consider the case of unitaries, which depend on some additional parameter, e.g., a spatial variable or a lattice position.More precisely, we assume the parameters E, n 1 , n 2 , n 3 to be functions of r, where r can be any set of either discrete or continuous variables.In this current study, we assume position r as a continuous variable.However, the discrete case can be obtained by replacing integrals with summations, i.e., f (x) dx → x f (x).
The unknown unitary process acts on quantum states whose Hilbert space is the tensor product of a qubit space H i , associated with an internal degree of freedom, and a high-dimensional space H r , associated with the lattice.
We consider an input state uniformly distributed along r, with an internal state prepared as one of the MUB states |a⟩.|a⟩ is assumed to be the positive eigenvector of σ 1 , σ 2 or σ 3 .The unknown unitary U (r) acts on this state, and then a projection onto the same internal state is performed.The probability distribution along r is I aa (r) = | ⟨a|U (r)|a⟩| 2 .Another measurement is performed to retrieve the probability distribution Ĩaa (k) in the reciprocal space of r.This is enabled by a QFT of the final state: One can extract the wave function phase in each plane from the two probability distributions by applying phaseretrieval techniques.We specifically employed the Gerchberg-Saxton (GS) algorithm [34].In particular, we can retrieve the amplitude A a (r) = I aa (r) and the phase α a (r) = arg( ⟨a|U (r)|a⟩).However, the latter is determined up to an unknown constant ξ a .The amplitude and phase are related to the process parameters according to where a = 1, 2, 3, depending on which Pauli matrix is considered (σ a |a⟩ = |a⟩).Thus, we obtain Equations ( 6a)-(6b) show that the extracted parameters depend on the global phase shifts ξ a , which cannot be estimated from the phase retrieval method.Indeed, any phase that differs from α a (r) by a constant global shift yields the same measured amplitude in the direct and reciprocal space.Considering the ambiguity due to the global phase shift, we list all the possible energy modulations compatible with the measurements.In practice, we select N values of the global phase shift ξ a,j = 2πj/N , with j = 0, 1, ..., N − 1.We specifically set N = 64.Only one of the candidates can best describe the process under investigation.To find it, we perform an additional measurement in the reciprocal space, obtained by evolving any input state without projection on the internal degree of freedom, e.g., where |b⟩ can be chosen arbitrarily.Crucially, this last measurement also provides the normalization factor for all the data.In this way, by numerically simulating the far field obtained from each of the N possible SU(2) evolutions, we can isolate the realization associated with the physical setup by sifting the one that minimizes the distance with the measurement of Eq. ( 7).We finally remark that, in principle, the method also works with only 2 sets of measurements.For instance, one could process I 11 (r) and I 22 (r), together with Ĩ11 (k) and Ĩ22 (k), to retrieve n 1 (r) and n 2 (r).In this case, the third component could be directly computed from the normalization condition: Combined with the final measurement in the reciprocal space, these add up to a minimal set of 5 total measurements, in perfect agreement with the argument provided in Ref. [18].However, we observed that experimental noise may yield nonphysical results, i.e., n 3 features a nonzero imaginary part at some location.This suggested integrating the minimal setting with an additional set of measurements [18,20,21].
When accounting for the noise, the proposed technique thus requires only 7 measurements instead of the conventional 16: 3 measurements in the lattice space, 3 in the reciprocal space, and a last one, still in the reciprocal space, to fix normalization and remove all the ambiguities on the parameters of the unitary.

III. EXPERIMENTAL RESULTS
In photonic setups, a qubit can be encoded into photon polarization, which is typically manipulated via optical waveplates.In the circular polarization basis, where |L⟩ = (1, 0) T and |R⟩ = (0, 1) T are left and right circular polarization states, respectively, a waveplate L δ,θ having the birefringence δ and the optic axis oriented at θ with respect to the horizontal direction can be expressed in the matrix form A single waveplate thus implements a rotation of an angle −δ/2 around the equatorial axis n = (cos 2θ, sin 2θ, 0).Nevertheless, one can cascade multiple waveplates to implement more general operations [35,36].
We apply FQPT to periodic polarization transformations induced by complex LCMSs [22].These can be modeled as optical waveplates having patterned opticaxis modulation θ = θ(x, y) and fixed, but tunable, birefringence [37].In particular, we test our method with LCMSs featuring high spatial-frequency modulations along the x and y directions.This scenario is experimentally more challenging than the one addressed in Ref. [18], where simple combinations of polarization gratings were considered.We benchmark our technique against the standard ML approach processing a whole set of 16 polarimetric measurements (all taken in the near field), taking into account both the timing and the accuracy of the final reconstruction.For the minimization, we employed the NMinimize routine from Wolfram Mathematica [38].
The experiments are realized with the setup sketched in Fig. 1a).A Ti:Sa laser (wavelength λ = 810 nm) is coupled to a single-mode fiber.The output Gaussian mode is magnified with a telescope lens system, f 1 = 125 mm and f 2 = 200 mm (not shown in the figure).The beam waist is measured to be 2.6 ± 0.1 mm.In this way, the overall beam size is larger than the largest periodicity on the plates, that is Λ = 2.5 mm.A combination of a linear polarizer (P 1 ), a half-wave plate (HWP 1 ) and a quarter-wave plate (QWP 1 ) is needed to prepare any input polarization state.The beam then propagates through one or multiple LCMSs, engineering a space-dependent SU(2) optical operator.Another set QWP 2 -HWP 2 -P 2 is adjusted to project onto any state.The last element (P 2 ) is removed when performing the measurement of Eq. (7).Each polarimetric measurement is collected on a CCD camera, placed either after a 4f system (f = 150 mm) or in the focal plane of a lens (f ′ = 250 mm), depending on if the measurement is realized in the near field or in the far field, respectively.Recall that the far-field light distribution is pro- portional to the transverse momentum distribution, i.e., to the Fourier transform of the input field.The first experiment is realized with a single LCMS displaying a complex periodic modulation along the x axis (see Fig. 1b)).The optical retardation is set at δ = π.The periodicity of the process simplifies the analysis, as only a discrete spectrum of Fourier components is expected to appear in the far field for all the measurements.Furthermore, in this one-dimensional (1D) realization, the intensity modulations can also be integrated along the y direction to mitigate experimental imperfections.
Figure 2a) shows the far-field intensity distributions recorded for the ⟨H|U |H⟩ configuration, with |H⟩ = (|L⟩ + |R⟩)/ √ 2. As discussed above, the farfield distribution can be discretized and a normalized power spectrum P (m) is extracted.Figures 2b)-c) illustrate the comparison between the experimentally measured total far-field distribution for a |L⟩ input state (see Eq. ( 7)) and the reconstruction performed via FQPT and ML, respectively.The agreement with the experimental observation is quantified in terms of the similarity estimator s = ( m P exp (m)P rec (m)) 2 , where P exp (m) and P rec (m) are the (normalized) experimental and reconstructed far-field distributions.We obtain s FQPT = 97.2% and s ML = 97.0%.
Another metric that can be considered is the absolute distance between the two distributions, computed as ∆ = ( m |P exp (m) − P rec (m)|) 2 .We report ∆ FQPT = 0.088 and ∆ ML = 0.083.In both figures, the insets show the same comparison for the 1D near-field amplitude I HH (x).Reconstructions of individual parameters are plotted in Fig. 2d).The agreement between the predictions of the two methods is quantified in terms of the fidelity F = Tr U † ML U FQPT /n, where n = 2 is the dimension of the internal degree of freedom [39].An excellent average fidelity is obtained, F = 95.7%,where F denotes the average fidelity computed over all the pixels.These results prove that both methods provide reliable reconstructions.It must be noted that FQPT achieves satisfactory reconstructions by only processing a near-optimal set of measurements.Moreover, a brute-force minimization approach tends to jump between the parameters associated with processes U and e iπ U = −U , as these both generate the same experimental outcomes [18,[40][41][42].For this reason, a continuity constraint must be enforced between consecutive pixels.Conversely, our method is assumption-free as it does not rely on any a-priori hypothesis on the process parameters.The technique is also extremely efficient on the computational level.The total times required for a complete reconstruction are t FQPT ≈ 1 min and t ML ≈ 30 min.
Recall that the GS algorithm is executed to retrieve the unknown phases.This algorithm is based on an iterative strategy that presents intrinsic limitations [43].The convergence speed is sensitive to the initial guess of the phase, and the convergence to a global minimum is not guaranteed.Moreover, the presence of noise in the input data can severely affect the accuracy of phase retrieval.To overcome these limitations, for each set of polarimetric measurements, the algorithm is run N 1D T = 100 times for N 1D I = 1000 iterations, randomly initializing the phase guess at each trial.In doing so, the best phase reconstruction can be selected (up to global shifts) as the one minimizing the total distance between the reconstructed and measured near-field amplitudes.The second experiment is realized by cascading the previous LCMS with a second one, patterned along the y direction (see Fig. 1b)).The sequence thus implements a two-dimensional (2D) periodic SU(2) transformation.The birefringence setting is δ 1 = π/2 and δ 2 = π.Figure 3a) shows the experimental far-field distribution ĨHH (k x , k y ), from which a 2D power spectrum P (m x , m y ) is extracted.The total far-field reconstructions for a |L⟩ input state are plotted in Fig. 3b), where the comparison with the experimental result is also provided.We obtain s FQPT = 88.5% and s ML = 96.7%, with a total error of ∆ FQPT = 0.080 and ∆ ML = 0.013.The slight deterioration in the final prediction can be ascribed to the reduced performance of the GS algorithm in two spatial dimensions [44].Figure 3c) shows the near-field amplitude I HH (x, y), as reconstructed from FQPT and ML, compared with the experimental measurement.In this case, the GS algorithm runs N 2D T = 50 times for N 2D I = 500 iterations.In this realization, we compress the experimental images by integrating light intensity over 11 pixel × 11 pixel regions equally distributed on the camera.This allows for both minimizing the errors due to local intensity fluctuations in the image area and keeping the computation time within the same range as the 1D experiment.The reconstructions of individual parameters are shown in Fig. 3d).An alternative visualization of the reconstructed process in terms of the local eigenvector and rotation angle is provided in Fig. 1c).A good agreement is observed between the two predictions, with average fidelity F = 91.4%.

IV. DISCUSSION AND CONCLUSIONS
In this work, a new technique for fast and accurate Quantum Process Tomography is demonstrated.This is accomplished via a non-interferometric scheme requiring no a-priori information on the unknown operator.In the case of complex polarization transformations, our method achieves performances very close to the standard tomography based on an overcomplete set of measure-ments.It offers experimental advantage, only requiring a near-minimal set of measurements, and computational speed-up, outperforming the standard approach by at least one order of magnitude.It appears naturally suitable for all experimental setups that allow for easy access to conjugate domains.For this reason, this method can also be implemented on other physical platforms, such as quantum circuits for atoms [45,46] and electron beams [47,48].
Nevertheless, we believe that the performance of the method can be further improved by adopting optimized strategies for phase retrieval.For instance, Convolutional Neural Networks represent a promising solution [49][50][51].We also expect that FQPT can detect non-unitary evolutions [52] if equipped with some minor modifications, such as including multiple intermediate planes in the final analysis.At the same time, it would be interesting to adapt the present method to retrieve multiphoton gates and complex operations in high-dimensional Hilbert spaces.
In principle, this procedure can also be applied to processes acting on a m × d dimensional space, having the irreducible form d i=1 U (i) , with U (i) ∈ SU(m).In this case, one must consider the decomposition of U (i) in the generators of SU(m), the generalized Gell-Mann matrices.Although the analytical relations between the process parameters and the measurement outcomes become significantly more complicated, the number of measurements required by FQPT will still be optimal.This will be investigated in successive works.

Figure 1 .
Figure 1.Fourier Quantum Process Tomography.a) A space-dependent polarization transformation U (x, y) is implemented via LCMSs.The process tomography is performed by preparing and projecting onto MUB states.The resulting intensity distributions are measured in the far field, Ĩaa(kx, ky), and image plane, Iaa(x, y), of the optical operator.b) Plots of the optic-axis pattern of two plates used to test the technique.These are patterned along orthogonal directions to create a complex 2D modulation.c) Reconstruction of the process U (x, y) = Tx(x)Ty(y) using FQPT and the traditional ML technique.The arrows represent the local eigenvectors n(x, y) and their color is associated with the local rotation angle E(x, y) (see also Fig. 3d) for details).

Figure 2 .
Figure 2. One-dimensional complex polarization transformations.a) Far-field intensity distribution ĨHH (kx) for a 1D periodic SU(2) transformation.A power spectrum P (m) is extracted after discretizing the reciprocal space.b)-c) Comparison between the experimental total far-field distribution for a |L⟩ input state and the reconstruction obtained from FQPT b) andML c).The insets show the reconstructed near-field amplitude IHH (x), compared with the experimental measurement.d) Reconstructions of the process parameters across one period performed by FQPT (red curves) and ML (blue curves).The expected periodicity is captured by both approaches.

Figure 3 .
Figure 3. Two-dimensional complex polarization transformations.a) Far-field intensity distribution ĨHH (kx, ky) for a 2D periodic SU(2) transformation.A power spectrum P (mx, my) is extracted.b) Comparison between the total far-field distribution for a |L⟩ input state and the reconstruction obtained from FQPT and ML.c) Reconstructed near-field amplitude IHH (x, y), compared with the experimental measurement.d) Reconstructions of the process parameters across one period performed by FQPT and ML.The expected periodicity is captured by both approaches.