Variational Simulation of Schwinger's Hamiltonian with Polarisation Qubits

The numerical emulation of quantum physics and quantum chemistry often involves an intractable number of degrees of freedom and admit no known approximations in a general form. In practice, representing quantum-mechanical states using available numerical methods become 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 NISQ architectures in identifying and quantifying quantum phase transitions with noisy qubits. We find that noise does not impede the detection of the phase transition point in a large range of noise levels.

The numerical emulation of quantum physics and quantum chemistry often involves an intractable number of degrees of freedom and admit no known approximations in a general form. In practice, representing quantum-mechanical states using available numerical methods become 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 NISQ architectures in identifying and quantifying quantum phase transitions with noisy qubits. We find that noise does not impede the detection of the phase transition point in a large range of noise levels.

I. INTRODUCTION
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 behavioralsciences [3][4][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 manybody 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 argued [6] that physical systems themselves naturally posses the computational capacity to be harnessed and used.
Variational approaches to optimization and simulation of eigenstates [7][8][9][10][11][12] have been used recently to port ideas from machine learning [13] to enhance algorithms with quantum processors [13][14][15]. These approaches rely on an iterative quantum-to-classical variational procedure. Proven to be a universal model of quantum computation in [16]-where the ansatz circuits are proven to be universal in [17]-the variational approach to quantum computation arose naturally as the pathway between a static simulator and a fully programmable 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.
Up-to-date experimental works realize variational algorithms on different quantum hardware including superconducting qubits [10,11], trapped atoms [8,9,12] and photonic quantum processors [7,18]. The most common application of quantum variational algorithms includes quantum chemistry problems. Despite the fact that the main purpose of the algorithm is finding eigenvalues and eigenvectors, [19] show that variational techniques can also find excited states, and various other proposals further expand the limits of applicability [20].
The variational quantum eigensolver (VQE) [7] performs classical optimization to minimize a mean Hamiltonian expected value-found by quantum hardware. 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 ∈ N and an experimental setup prepares a quantum state |ψ(θ) parametrized by these control parameters. 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: Therefore, the problem consists of using classical optimization algorithms to select optimal parameters θ corresponding to the (ideally) minimal value of energy.
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 decoherence [21]  identifying and quantifying quantum phase transitions with noisy qubits.

II. ENCODING THE ANSATZ STATE IN POLARIZATION QUBITS
Ansatz state preparation is strongly connected with the particular experimental realization, because elements with tunable parameters used in the setup determine the parametrization of the ansatz. Therefore, we start from the experimental scheme to clarify the origin of the chosen ansatz.
We implement a VQE algorithm using polarizationencoded qubits. The experimental setup scheme is shown in Fig. 1. The initial state preparation is carried out by a two-photon source based on spontaneous parametric down-conversion process (SPDC) in the Sagnac interferometer [22]. 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 KTP nonlinear crystal (PPKTP) 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: 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 singlemode 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. 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 QWP and HWP retarders, 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 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.
measurement basis are compiled together, e. g., here U (θ) is a transformation of a corresponding waveplate with an axis angle θ and B is a unitary matrix that changes the basis. New angles θ 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 (δ, θ): A controlled-X gate corresponds to the SPDC of the pump photon in the nonlinear crystal. Taking into account 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 U 1 ⊗ U 2 is applied to get the probe state |ψ(θ) : Unitaries U 1 and U 2 are composed of the waveplate transformations: 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 reader is referred to Appendix A for details.

III. THE SCHWINGER MODEL
The Schwinger model describes interactions between Dirac fermions via photons in a two-dimensional space. 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 H N describes electronpositron pair creation and annihilation, their interaction and takes into account the particle mass: (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 where we set the background electric field parameter 0 to zero. The problem Hamiltonian can be encoded in the multiqubit system by using its decomposition into Pauli strings: where N denotes the number of qubits and h α ∈ R 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 The quantum phase transition manifests itself in the behavior of the order parameter For polarization-encoded pair of qubits the order parameter is simply a projector onto |V H state: Two-qubit Schwinger Hamiltonian has four nondegenerate eigenvalues E 1 , . . . , E 4 . Two intermediate eigenvalues, E 2 = 2 and E 3 = 1, are constant and do not depend on the mass m. The largest and the smallest eigenvalues, E 1,4 = 1/2 ± m 2 + 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 E min ≡ E 4 . The graph of its dependence on m is depicted in Fig. 3a. Also Fig. 3b shows the order parameter versus mass. The solid lines correspond to the exact analytical solutions, 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.
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 m = −1/2 the minimum has a lower value, E min < −3/2, while residing near the plateau with E = −3/2. So the landscape of E(θ) in the punctured neighbourhood of m = −1/2 becomes a flat valley. A valley landscape puzzles the gradient-based optimizers and significantly slows convergence [23]. Therefore, the algorithm terminates at the wrong value. A little step noticeable in Fig. 3b illustrates this situation. When m is far away from the phase transition point, the plateau does not strongly influence the results, because E min 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 |Ψ − : where W 1,2 are some single-qubit unitary matrices. Consequently, an arbitrary Bell state |ψ 0 is invariant under the following transformation: 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 10 5 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 Appendix B for details).

IV. EFFECTS OF NOISE
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 (LCVR) that 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 where E j 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 in our experiment. Our experimental setup allows to explore 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 versus 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).
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 equal 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.

V. CONCLUSION
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 quantum approximate optimization algorithm (QAOA) [15], because it uses specific ansatz adjusted for the target Hamiltonian. The authors declare no competing interests. Data and code availability: The data that supports this study are available within the article. The code for generating the data will be made available on GitHub after this paper is published.

Appendix A: Classical optimizer
The target function under minimization E(θ) is the mean Hamiltonian value, but in the experiment, only random samples of E(θ) obtained by repetitive measurements are available. So experimental VQE is a stochastic optimization problem. We use a simultaneous perturbation stochastic approximation (SPSA) algorithm [24] as a classical optimizer in our VQE implementation. It is useful for high-dimensional problems, where the gradient of the objective function is not directly available, because SPSA requires only two function evaluations per iteration for any number of parameters in the optimization problem.
Single SPSA iteration proceeds as follows: 1. Generate a random vector ∆ with elements being ±1 with equal probability.
2. Estimate a gradient g: 3. Move to the new point θ : Scalar variables a and b are called meta parameters. The parameter a describes the iteration step and b defines finite difference to calculate the gradient. They change with the number of iterations k according to schedule: Usually final values a f and b f are set to zero to ensure convergence in the limit k → ∞. However, we use nonzero a f and b f to track a slow drift of the experimentally prepared probe state |ψ(θ) over time [25]. The drift occurs mainly due to the instability of polarization transformation in optical fibers connecting the SPDC source and measurement part of the setup.
Moreover, we find out influence of mass parameter m on VQE convergence-closeness to phase transition makes it slowly. So we adjust meta parameters for each m as We carried out numerical simulations of the VQE algorithm for m = 0, −1/2, 1, 10 to investigate how the choice of an initial point θ 0 affects convergence and explore the set of obtained solutions. Recall that the Schwinger Hamiltonian H 2 (9) undergo phase transition of the order parameter at m = −1/2, so points m = 0 and m = 1 are nearby and symmetric w. r. t. phase transition and m = 10 is an example of a distant point. To collect statistics, we execute the VQE algorithm 10 5 times for each m starting from freshly generated random initial points θ 0 . The points are distributed uniformly in a sixdimensional hypercube with the side length equal to π, which coincides with the period of the target function E(θ).
Each VQE run results in the final point θ, the energy level E(θ) (eigenvalue), and the order parameter O .  As one can see, there is a sharp peak near a wrong value O = 1/2 for both histograms and obtuse peaks approaching true solutions O ≈ 0.38 and O ≈ 0.62 for m = 0 and m = −1, respectively. As it was said in the main text, O = 1/2 corresponds to the plateau in landscape of the target function, which acts as an attractor for the optimizer. Fig. 6 illustrates evolution of found-eigenvalue distribution for different m. As expected, the histogram for m = 10 is unimoal and centered around the exact eigenvalue for the corresponding Hamiltonian. When m approaches phase transition point m = −1/2, the second peak emerges around E = −3/2, which is precisely the Hamiltonian eigenvalue for m = −1/2. This erroneous peak is small for m = 1, but it becomes even higher than the true one for m = −1. Closeness to phase transition point changes convergence statistics dramatically-less than 1% of simulations reveal proper values for m = −1.
We tried to clarify the structure of obtained solutions in the space of tuned parameters θ. First, we bring all found points θ to a hypercube that corresponds to the target function period. After that, for graphic purposes, we decreased the dimensionality of obtained solutions θ from six to three using principal component analysis (PCA). PCA finds a lower-dimensional hyperplane in the original space, which has a minimum average squared distance from the points to the hyperplane. Then points are projected to the approximating hyperplane. PCA helps to keep the real structure of the initial space and find any clusters of points with the same values. Fig. 7 presents obtained PCA projections for m = −0.5, 0, and 1 in two views, which will be called "top" and "side" for convenience. Color shows target function values E(θ). Blue points correspond to erroneous eigenvalues for the given m. The overall structure is similar for different m, especially for −0.5 and 0. However, the majority of converged values are not correct for m = 0. For m = 1, the fraction of good solutions increases, blue areas slowly disappear. This suggests that the algorithm converges to the desired point with higher probability, which is in perfect agreement with the histogram in Fig.  6b. There are two classes of proper solutions: points from the first class are located "inside" the erroneous plateau and for the second lie "outside". However, it appears difficult to isolate areas of initial points θ 0 that can guarantee finding the true minimum or lead to one or another class of solutions, and additional research is required.