The numerical emulation of quantum physics and quantum chemistry often involves an intractable number of degrees of freedom and admits no known approximation in the general form. In practice, representing quantum-mechanical states using available numerical methods becomes exponentially more challenging with increasing system size. Recently, quantum algorithms implemented as variational models have been proposed to accelerate such simulations. Here, we study the effect of noise on the quantum phase transition in the Schwinger model within a variational framework. The experiments are built using a free space optical scheme to realize a pair of polarization qubits and enable any two-qubit state to be experimentally prepared up to machine tolerance. We specifically exploit the possibility to engineer noise and decoherence for polarization qubits to explore the limits of variational algorithms for noisy intermediate-scale quantum architectures in identifying and quantifying quantum phase transitions with noisy qubits. We find that despite the presence of noise, one can detect the phase transition of the Schwinger Hamiltonian even for a two-qubit system using variational quantum algorithms.

The numerical emulation of quantum systems underpins a wide assortment of science and engineering and touches on fields ranging from statistical and quantum physics to biology and even to life-1,2 and behavioral-sciences.3–5 A physical simulator bootstraps one physical system to emulate the properties of another. While the time and memory required in the simulation of physical systems, particularly strongly correlated many-body quantum systems, using traditional computers often scales exponentially in the system size, the same is not always true for the physics-based quantum simulator. Indeed, Richard Feynman first speculated that instead of viewing the simulation of quantum systems using classical computers as a no-go zone due to its apparent computational difficulty, Feynman argued6 that physical systems themselves naturally posse computational capacity to be harnessed and used.

Variational approaches to optimization and simulation of eigenstates7–14 have been used recently to port ideas from machine learning15 to enhance algorithms with quantum processors.13,15–18 These approaches rely on an iterative quantum-to-classical variational procedure. Proven to be a universal model of quantum computation in Ref. 19—where the ansatz circuits are proven to be universal in Ref. 20—the variational approach to quantum computation arose naturally as the pathway between a static simulator and a fully programable gate-based quantum information processor. The variational model of quantum computation is the algorithmic workhorse of the current NISQ (noisy intermediate-scale quantum) technology era.

Recent experiments realize variational algorithms on different quantum hardware including superconducting qubits,10,11 trapped atoms,8,9,12 and photonic quantum processors.7,21 The most common application of quantum variational algorithms includes quantum chemistry applications.22,23 The original purpose of the algorithm was finding ground eigenvalues and eigenvectors; Ref. 24 shows that variational techniques can also find excited states, and various other proposals further expand the limits of applicability.25 

The variational quantum eigensolver (VQE)7 performs classical optimization to minimize an expected Hamiltonian value. The purpose of this algorithm is to determine the eigenvalues of a particular Hamiltonian, which describes a physical system, for example, the interaction of spins or electronic systems.8,9 A classical computer initially sets a vector of parameters θ={θi} for i, and an experimental setup prepares a parameterized quantum state |ψ(θ). After that, the state is measured, and the evaluation of the mean Hamiltonian value occurs. The parameters θ are adjusted to find the ground-state energy

Emin(θ)=minθψ(θ)|H|ψ(θ).
(1)

Therefore, the problem consists of using classical optimization algorithms to select optimal parameters θ corresponding to the (ideally) minimal value of energy (1).

Here, we report an experimental implementation of VQE in a photonic system. We target the exploration of a quantum phase transition in the Schwinger model. We specifically exploit the possibility to engineer noise and decoherence26 for polarization qubits to explore the limits of variational algorithms for NISQ architectures in identifying and quantifying quantum phase transitions with noisy qubits.

We implement VQE with polarization-encoded qubits using the experimental scheme shown in Fig. 1. Initial state preparation is carried out by a two-photon source based on a spontaneous parametric downconversion process (SPDC) in the Sagnac interferometer.27 A 405-nm laser diode beam is divided by a polarization beam splitter (PBS), which makes it possible to pump a 30-mm long periodically poled potassium titanyl phosphate nonlinear (PPKTP) crystal in two opposite directions. As a result of a type-II SPDC, pairs of signal and idler photons with orthogonal polarizations are generated in both directions. Then each photon pair is divided on the PBS and sent to different arms of the scheme. Thus, at the output of the two-photon source, we have the following entangled state:

|ψin=α(θ1,θ2)|HV+β(θ1,θ2)|VH,
(2)

where the coefficients α and β depend on the angular positions θ1 and θ2 of waveplates QWP1 and HWP1, which are placed in the pump beam. By rotating QWP1 and HWP1, we can alter the degree of entanglement of the initial state. The photon pairs are coupled to single-mode fibers and transferred to the measurement part of the setup. Motorized quarter-wave (QWP2, QWP3) and half-wave (HWP3, HWP4) plates are placed in each arm after the single-mode fiber channel, allowing to obtain any polarization state at the output. Finally, the Wollaston prism spatially separates the vertical and the horizontal polarizations to detect the prepared states using single-photon detectors in each of the arms. According to the measurement results, the classical algorithm transfers the new parameter values to the motorized plates until the optimal set of parameters is obtained.

FIG. 1.

Experimental setup implementing the variational quantum eigensolver algorithm using a pair of polarization qubits. Half-wave plates (HWP3, HWP4) and quarter-wave plates (QWP2, QWP3) in each channel prepare the desired variational state |ψ(θ). Wollaston prisms implement projective measurements in the computational basis (horizontal and vertical polarizations). HWP1 and QWP1 control the polarization of the pump beam, while HWP2 rotates the polarization by 90° to ensure the same polarization of the pump for the two pumping directions for the nonlinear crystal (PPKTP) inside the Sagnac interferometer. Liquid crystal variable retarders (LCVRs) are used to artificially introduce dephasing noise on the first qubit, when necessary.

FIG. 1.

Experimental setup implementing the variational quantum eigensolver algorithm using a pair of polarization qubits. Half-wave plates (HWP3, HWP4) and quarter-wave plates (QWP2, QWP3) in each channel prepare the desired variational state |ψ(θ). Wollaston prisms implement projective measurements in the computational basis (horizontal and vertical polarizations). HWP1 and QWP1 control the polarization of the pump beam, while HWP2 rotates the polarization by 90° to ensure the same polarization of the pump for the two pumping directions for the nonlinear crystal (PPKTP) inside the Sagnac interferometer. Liquid crystal variable retarders (LCVRs) are used to artificially introduce dephasing noise on the first qubit, when necessary.

Close modal

We should note that estimation of a single mean value of a Hamiltonian requires projective measurements in several bases, while the Wollaston prism projects only onto |H and |V states. To change the basis, one may use an additional pair of QWPs and HWPs, mounted just before the Wollaston prism. However, we chose a more economic setup, where the local unitary transformation of the initial state and the transformation of the measurement basis are compiled together, e.g.,

H|BUHWP(θ4)UQWP(θ3)=H|UHWP(θ4)UQWP(θ3).
(3)

Here, U(θ3,4) is a transformation of a corresponding waveplate with an axis angles θ3,4, and B is a unitary matrix that changes the basis. New angles θ3,4 are calculated automatically in our algorithm to perform measurements in desired bases.

By mapping experimental optical elements to the gate model, we arrive at the ansatz preparation circuit with six tunable parameters θi,i=1,,6, which is presented in Fig. 2. The parameters θi physically correspond to the waveplates' rotation angles. A general waveplate with a phase shift δ and an axis position θ performs the transformation U(δ,θ)

U(δ,θ)=V(θ)D(δ)V(θ),V(θ)=1cos(θ)ıσysin(θ),D(δ)=e+ıδ|11|.
(4)

A controlled-X gate corresponds to SPDC in the nonlinear crystal.

FIG. 2.

Schematic of the VQE algorithm. The initial state is prepared by three single-qubit gates and a Controlled-X gate. UQWP(θ1) and UHWP(θ2) are used to control the initial state in experiment. Other four single-qubit transformations serve both for the ansatz state preparation and the measurement basis change.

FIG. 2.

Schematic of the VQE algorithm. The initial state is prepared by three single-qubit gates and a Controlled-X gate. UQWP(θ1) and UHWP(θ2) are used to control the initial state in experiment. Other four single-qubit transformations serve both for the ansatz state preparation and the measurement basis change.

Close modal

Taking into account the ansatz preparation scheme, our VQE algorithm implementation consists of the four main steps:

  1. SPDC source emits the initial entangled state |ψin (2).

  2. Once the initial state has been prepared, a local unitary transformation U1U2 is applied to get the probe state |ψ(θ)
    |ψ(θ)=U1(θ3,θ4)U2(θ5,θ6)|ψin(θ1,θ2).
    (5)

    Unitaries U1 and U2 are composed of the waveplate transformations: U1=UHWP(θ4)UQWP(θ3),U2=UHWP(θ6)UQWP(θ5).

  3. The cost function E(θ)=ψ(θ)|H|ψ(θ) is calculated by summing up measurement results with coefficients depending on the problem Hamiltonian. Since usually Hamiltonian is expressed as a linear combination of Pauli observables and our setup allows only projective measurements, we first should decompose the Hamiltonian as a linear combination of projectors onto eigenbases of Pauli matrices. Change of basis is carried out according to the rule (3).

  4. The value of E(θ) is minimized as a function of the parameters θ using a classical optimizer routine. In particular, we use simultaneous perturbation stochastic approximation (SPSA) algorithm.

The Schwinger model describes interactions between Dirac fermions via photons in a two-dimensional space.33,34 In Ref. 12, the authors map the model to the lattice model of an electron-positron array. The Schwinger Hamiltonian exhibits a quantum phase transition: the signature of which (in finite dimensions) allows us to determine new features in VQE behavior and clarify its robustness to noise.

The Schwinger Hamiltonian HN describes electron-positron pair creation and annihilation, their interaction and takes into account the particle mass

HN=wj=1N1[σj+σj+1+H.c.]+m2j=1N(1)jσjz+gj=1NLj2.
(6)

It consists of the three terms: the first one is responsible for the interaction of an electron and a positron, the second depends on bare mass m of the particles, and the third stands for the energy of the electric field. We assume the coefficients w=g=1 and only consider the dependence of the Hamiltonian ground energy on the bare mass. The operators in the third term are given by

Lj=ε012l=1j[σlz+(1)l],
(7)

where we set the background electric field parameter ε0 to zero.

The problem Hamiltonian can be encoded in the multi-qubit system by using its decomposition into Pauli strings: Pα=σ1α1σ2α2σNαN with single-qubit Pauli operators σiαi{I,σix,σiy,σiz} as

H=ahαPα,
(8)

where N denotes the number of qubits, and hα are real coefficients. In further consideration, we will use this representation. We carried out numerical simulations and experiments for the case of two qubits, for which the Schwinger Hamiltonian takes the form

H2=1+σ1xσ2x+σ1yσ2y12σ1z+12σ1zσ2z+m2(σ2zσ1z).
(9)

The quantum phase transition manifests itself in the behavior of the order parameter

O=12N(N1)j>i(1+(1)iσiz)(1+(1)jσjz).
(10)

For polarization-encoded pair of qubits, the order parameter is simply a projector onto |VH state

O=14(1σ1z)(1+σ2z)=|VHVH|.
(11)

Two-qubit Schwinger Hamiltonian has four non-degenerate eigenvalues E1,,E4. Two intermediate eigenvalues, E2=2 and E3=1, are constant and do not depend on the mass m. The largest and the smallest eigenvalues, E1,4=1/2±m2+m+17/4, vary with m in a symmetric manner.

We are interested in the ground energy of the Hamiltonian that corresponds to the minimal eigenvalue EminE4. The graph of its dependence on m is depicted in Fig. 3(a), and Fig. 3(b) shows the order parameter vs mass. The solid lines correspond to the exact analytical solutions, and dots represent the results of simulations and experiment. A phase transition signifies itself in the rapid change of the order parameter from one to zero, and it is expected near the point m=1/2, where O=1/2.

FIG. 3.

Eigenvalue (a) and order parameter (b) vs bare mass m. Solid lines—analytical solution, cyan points (violin plot)—simulations, and red points—experimental result. The probability distribution for each m in simulations is obtained using 30 trials. For the experiment, all attempts are shown as distinct points.

FIG. 3.

Eigenvalue (a) and order parameter (b) vs bare mass m. Solid lines—analytical solution, cyan points (violin plot)—simulations, and red points—experimental result. The probability distribution for each m in simulations is obtained using 30 trials. For the experiment, all attempts are shown as distinct points.

Close modal

Exactly at the vicinity of the phase transition point m=1/2, we found a discrepancy between analytical solutions and VQE simulations. The Hamiltonian H2(m=1/2) has the ground energy Emin=3/2 with the corresponding eigenvector (|01|10)/2, which is a maximally entangled singlet state |Ψ. A distinguishing feature of the singlet state is its invariance under local unitary rotations USU(2): |Ψ=(UU)|Ψ. Therefore, the target function E(θ) remains constant on some parameter manifold. Note that this plateau does not change with the mass m, because m:Ψ|H2(m)|Ψ=3/2.

For m=1/2, the existence of the plateau is not a problem for the optimization algorithm, since the global minimum is attained at any point of the plateau. But for any m1/2, the minimum has a lower value, Emin<3/2, while residing near the plateau with E=3/2. So the landscape of E(θ) in the punctured neighborhood of m=1/2 becomes a flat valley. A valley landscape puzzles the gradient-based optimizers and significantly slows convergence.28 Therefore, the algorithm terminates at the wrong value. A little step noticeable in Fig. 3(b) illustrates this situation. When m is far away from the phase transition point, the plateau does not strongly influence the results, because Emin is much lower than E=3/2.

In our particular case, slow convergence originated from the invariance of the singlet state |Ψ being the Hamiltonian eigenvector for m=1/2. A more general view on the cause of the convergence problem is that it appears any time, when the ansatz is general enough to perform arbitrary local unitary transformations, and the Hamiltonian ground state is close to some Bell state (not necessarily |Ψ). Indeed, all Bell states are equivalent under local transformations, so we can find a local map that brings a Bell state |ψ0 to a singlet one |Ψ

|Ψ=(W1W2)|ψ0,
(12)

where W1,2 are some single-qubit unitary matrices. Consequently, an arbitrary Bell state |ψ0 is invariant under the following transformation:

USU(2):|ψ0=(W1UW1W2UW2)|ψ0.
(13)

If ansatz circuit is general enough to prepare different transformations of the form (13), then the plateau in the landscape of E(θ) appears. Therefore, when the Hamiltonian ground state is close to the Bell state, the nearby plateau will create flat valley landscape.

The simplest opportunity to get around poor optimizer convergence is by a correct choice of the initial point. We gathered statistics for 105 random initial points θ for m=1/2, 0, 1/2, and 10 and found that near the phase transition the algorithm sticks to the plateau much frequently than to the proper minimum (see the supplementary material for details).

In order to clarify the issue with the accuracy, we used parameter δ from Ref. 29. This parameter characterizes the closeness of the obtained energy E to the exact ground level E0 compared with the distance to the next energy level E1: δ=EE0E1E0. For “good-enough” accuracy, the parameter should be much less than one, δ1. In our work, the maximum value of δ is 0.176 for the experiment and 0.02 for simulations.

Compared to other types of quantum computers, photon circuits have low intrinsic noise levels. This means that we can add noise to the system in a controlled manner and get the dependencies of the parameters of interest on the noise level. We took advantage of this to evaluate the effect of noise on the phase transition that we observed without the noise. We expect that as the degree of dephasing increases, the phase transition will blur until it disappears completely. This will allow us to estimate the acceptable noise level in the system implementing VQE to identify quantum phase transitions.

The origin of the noise model used is connected with our experimental implementation. We artificially introduce noise to the system with liquid crystal variable retarders (LCVRs) that are placed directly before Wollaston prisms adding noise to measurement. Placing them in the ansatz preparation part would require taking into account fiber transformations and would allow the algorithm to compensate the noise. LCVRs allow us to change the phase of the specific polarization component of the light field. If the phase shift δ varies during the data acquisition time, then this leads to effective decoherence of the system state. The noise channel E(ρ) is thus the transformation (4) averaged over δ taken from some interval (depending on the noise strength). The explicit action of the noise channel is

ρ=E(ρ)=j=12EjρEj,Ei=V(θ)Di(δ)V(θ),D1(δ)=2ε2(eiδ001),D2(δ)=ε2(eiδ001),
(14)

where Ej are the Krauss operators, θ is a LCVR axis angle, and δ is a mean retardance. Noise strength is controlled by the parameter ε, 0ε1. We set θ=π/4 and δ=2π in our experiment.

Our experimental setup allows exploring the effect of this noisy channel on one qubit or simultaneously on both. Primarily, we simulated these two cases for different noise levels ε ranging from 0 to 1 with a 0.1 step to obtain the eigenvalues and the values of the order parameter vs m. As expected, the presence of noise in the system prevents the algorithm from converging to the exact eigenvalue, and noise escalation leads to convergence deterioration (Fig. 4).

FIG. 4.

Noise simulations for dephasing of one qubit (a) and (c) and both qubits (b) and (d). Figures (a) and (b) show the minimal eigenvalue dependence on m, and figures (c) and (d)—the dependence of the order parameter. Red lines correspond to noiseless simulations, and the color blur corresponds to the increase in noise strength ε from 0.1 to 1 in 0.1 steps. Points—experimental results and solid areas—theoretical prediction for points with different noise levels.

FIG. 4.

Noise simulations for dephasing of one qubit (a) and (c) and both qubits (b) and (d). Figures (a) and (b) show the minimal eigenvalue dependence on m, and figures (c) and (d)—the dependence of the order parameter. Red lines correspond to noiseless simulations, and the color blur corresponds to the increase in noise strength ε from 0.1 to 1 in 0.1 steps. Points—experimental results and solid areas—theoretical prediction for points with different noise levels.

Close modal

Finding appropriate eigenvalue becomes challenging for the case of simultaneous dephasing in both channels, and full dephasing (ε = 1) leads to degeneracy—the algorithm converges to 1 for any m. The phase transition in the order parameter blurs with increasing noise and disappears for ε = 1. Full dephasing makes the order parameter constant and equals to O=1/4 for any m. In the case of a single noise channel, the phase transition remains visible even with ε = 1, while the maximum value of order parameter is halved.

Quantum phase transitions as metal-insulator transition and transition between quantum Hall liquid states can be predicted and inquired by quantum algorithms. As we experimentally demonstrated, noise does not impede the detection of the phase transition point in a large range of noise levels. Only completely dephasing channels acting on both qubits prevent finding it in our model. This result demonstrates the noise-tolerance of VQE not only from speed and quality of convergence perspective but also from a practical point of view of determining the parameters of the Hamiltonian corresponding to a quantum phase transition.

We observe slow VQE convergence near the phase transition point and connect this behavior with the Hamiltonian ground state's closeness to the two-qubit singlet state. It seems to be a common effect for a combination of sufficiently general ansatz circuits and Hamiltonians, where the ground state exhibits additional symmetry. This hypothesis should be verified in future research. Possible approaches to circumvent poor convergence may include QAOA (Quantum Approximate Optimization Algorithm),17 because it uses specific ansatz adjusted for the target Hamiltonian.

Scalability is a major challenge for all modern quantum computing platforms, with the photonic one not being an exception. The experimental approach taken here may be relatively straightforward scaled up to 6–10 photons, and experiments on such scales are feasible.30 When the system is scaled up to larger number of photons, polarization encoding and free-space implementation used here are probably not the best option, and one should aim at integrated photonic circuits. Here, we may note that large fully programable circuits are available, and they can be used to realize parametrized transformations for variational algorithms,21 The exact forms of optimal variational ansatz for such encoding are not yet known and are an area of future research. However, a standard practice of using dual-rail encoded qubits allows one to realize variational algorithms as described here at an expense of finite probability of multi-qubit operations, which may be brought close to unity with an addition of extra photons.31 

The major challenge on the experimental side for large scale experiments will be photon loss in the circuit, which dramatically reduces the count rate for multi-photon events. So on short time scales, one should aim for higher brightness single photon sources and low-loss integrated optics. On a longer timescale, a fully integrated modular architecture for photonic computing may be developed, which has intrinsic loss tolerance, since the path length of each photon becomes independent of the circuit size due to the modular structure of the processor.32 

See the supplementary material for the detailed information about classical optimizer and the algorithm convergence statistical research for random initial point.

The Skoltech team acknowledges support from the research project, Leading Research Center on Quantum Computing (Agreement No. 014/20). The MSU team acknowledges financial support from the Russian Foundation for Basic Research (RFBR Project No. 19-32-80043 and RFBR Project No. 19-52-80034) and support under the Russian National Technological Initiative via MSU Quantum Technology Centre. The authors declare no competing interest.

The data that support the findings of this study are available within the article. The code for generating the data will be made available on GitHub after this paper is published.

1.
N.
Lambert
,
Y.-N.
Chen
,
Y.-C.
Cheng
,
C.-M.
Li
,
G.-Y.
Chen
, and
F.
Nori
,
Nat. Phys.
9
,
10
(
2013
).
2.
F.
Neukart
,
G.
Compostella
,
C.
Seidel
,
D.
Von Dollen
,
S.
Yarkoni
, and
B.
Parney
,
Front. ICT
4
,
29
(
2017
).
3.
T.
Werlang
,
C.
Trippe
,
G.
Ribeiro
, and
G.
Rigolin
,
Phys. Rev. Lett.
105
,
095702
(
2010
).
4.
J.
Abadie
,
B. P.
Abbott
,
R.
Abbott
,
T. D.
Abbott
,
M.
Abernathy
,
C.
Adams
,
R.
Adhikari
,
C.
Affeldt
,
B.
Allen
,
G.
Allen
 et al,
Nat. Phys.
7
,
962
(
2011
).
5.
M.
Lubasch
,
J.
Joo
,
P.
Moinier
,
M.
Kiffner
, and
D.
Jaksch
,
Phys. Rev. A
101
,
010301
(
2020
).
6.
R. P.
Feynman
,
Found. Phys.
16
,
507
(
1986
).
7.
A.
Peruzzo
,
J.
McClean
,
P.
Shadbolt
,
M.-H.
Yung
,
X.-Q.
Zhou
,
P. J.
Love
,
A.
Aspuru-Guzik
, and
J. L.
O'brien
,
Nat. Commun.
5
,
4213
(
2014
).
8.
M.-H.
Yung
,
J.
Casanova
,
A.
Mezzacapo
,
J.
Mcclean
,
L.
Lamata
,
A.
Aspuru-Guzik
, and
E.
Solano
,
Sci. Rep.
4
,
3589
(
2015
).
9.
Y.
Shen
,
X.
Zhang
,
S.
Zhang
,
J.-N.
Zhang
,
M.-H.
Yung
, and
K.
Kim
,
Phys. Rev. A
95
,
020501
(
2017
).
10.
P. J.
O'Malley
,
R.
Babbush
,
I. D.
Kivlichan
,
J.
Romero
,
J. R.
McClean
,
R.
Barends
,
J.
Kelly
,
P.
Roushan
,
A.
Tranter
,
N.
Ding
 et al,
Phys. Rev. X
6
,
031007
(
2016
).
11.
A.
Kandala
,
A.
Mezzacapo
,
K.
Temme
,
M.
Takita
,
M.
Brink
,
J. M.
Chow
, and
J. M.
Gambetta
,
Nature
549
,
242
(
2017
).
12.
C.
Kokail
,
C.
Maier
,
R.
van Bijnen
,
T.
Brydges
,
M. K.
Joshi
,
P.
Jurcevic
,
C. A.
Muschik
,
P.
Silvi
,
R.
Blatt
,
C. F.
Roos
 et al,
Nature
569
,
355
(
2019
).
13.
A. V.
Uvarov
,
A. S.
Kardashin
, and
J. D.
Biamonte
,
Phys. Rev. A
102
,
012415
(
2020
).
14.
D.
Wang
,
O.
Higgott
, and
S.
Brierley
,
Phys. Rev. Lett.
122
,
140504
(
2019
).
15.
J.
Biamonte
,
P.
Wittek
,
N.
Pancotti
,
P.
Rebentrost
,
N.
Wiebe
, and
S.
Lloyd
,
Nature
549
,
195
(
2017
); arXiv:1611.09347 [quant-ph].
16.
V.
Akshay
,
H.
Philathong
,
M.
Morales
, and
J.
Biamonte
,
Phys. Rev. Lett.
124
,
090504
(
2020
).
17.
E.
Farhi
,
J.
Goldstone
, and
S.
Gutmann
, “
A quantum approximate optimization algorithm
,” arXiv:1411.4028 (
2014
).
18.
K.
Mitarai
,
T.
Yan
, and
K.
Fujii
,
Phys. Rev. Appl.
11
,
044087
(
2019
).
19.
J.
Biamonte
, “
Universal variational quantum computation
,”
Phys. Rev. A
103
(
3
),
L030401
(
2019
).
20.
M. E. S.
Morales
,
J. D.
Biamonte
, and
Z.
Zimborás
,
Quantum Inf. Process.
19
,
291
(
2020
).
21.
J.
Carolan
,
M.
Mohseni
,
J. P.
Olson
,
M.
Prabhu
,
C.
Chen
,
D.
Bunandar
,
M. Y.
Niu
,
N. C.
Harris
,
F. N.
Wong
,
M.
Hochberg
 et al,
Nat. Phys.
16
,
322
(
2020
).
22.
I. G.
Ryabinkin
,
S. N.
Genin
, and
A. F.
Izmaylov
,
J. Chem. Theory Comput.
15
,
249
(
2019
).
23.
R. M.
Parrish
,
E. G.
Hohenstein
,
P. L.
McMahon
, and
T. J.
Martínez
,
Phys. Rev. Lett.
122
,
230401
(
2019
).
24.
J. R.
McClean
,
S.
Boixo
,
V. N.
Smelyanskiy
,
R.
Babbush
, and
H.
Neven
,
Nat. Commun.
9
,
4812
(
2018
).
25.
O.
Higgott
,
D.
Wang
, and
S.
Brierley
,
Quantum
3
,
156
(
2019
).
26.
27.
A.
Fedrizzi
,
T.
Herbst
,
A.
Poppe
,
T.
Jennewein
, and
A.
Zeilinger
,
Opt. Express
15
,
15377
(
2007
).
28.
H. H.
Rosenbrock
,
Comput. J.
3
,
175
(
1960
).
29.
C.
Bravo-Prieto
,
J.
Lumbreras-Zarapico
,
L.
Tagliacozzo
, and
J. I.
Latorre
,
Quantum
4
,
272
(
2020
).
30.
X.-L.
Wang
,
L.-K.
Chen
,
W.
Li
,
H.-L.
Huang
,
C.
Liu
,
C.
Chen
,
Y.-H.
Luo
,
Z.-E.
Su
,
D.
Wu
,
Z.-D.
Li
,
H.
Lu
,
Y.
Hu
,
X.
Jiang
,
C.-Z.
Peng
,
L.
Li
,
N.-L.
Liu
,
Y.-A.
Chen
,
C.-Y.
Lu
, and
J.-W.
Pan
,
Phys. Rev. Lett.
117
,
210502
(
2016
).
31.
P.
Kok
,
W. J.
Munro
,
K.
Nemoto
,
T. C.
Ralph
,
J. P.
Dowling
, and
G. J.
Milburn
,
Rev. Mod. Phys.
79
,
135
(
2007
).
32.
S.
Bartolucci
,
P.
Birchall
,
H.
Bombin
,
H.
Cable
,
C.
Dawson
,
M.
Gimeno-Segovia
,
E.
Johnston
,
K.
Kieling
,
N.
Nickerson
,
M.
Pant
 et al, arXiv:2101.09310 (
2021
).
33.
T. M. R.
Byrnes
,
P.
Sriganesh
,
R. J.
Bursill
, and
C. J.
Hamer
,
Phys. Rev. D
66
,
013002
(
2002
).
34.
T.
Byrner
and
Y.
Yamamoto
,
Phys. Rev. A
73
,
022328
(
2006
).

Supplementary Material