Quantum Simulations for Strong-Field QED

Quantum field theory in the presence of strong background fields contains interesting problems where quantum computers may someday provide a valuable computational resource. In the NISQ era it is useful to consider simpler benchmark problems in order to develop feasible approaches, identify critical limitations of current hardware, and build new simulation tools. Here we perform quantum simulations of strong-field QED (SFQED) in $3+1$ dimensions, using real-time nonlinear Breit-Wheeler pair-production as a prototypical process. The strong-field QED Hamiltonian is derived and truncated in the Furry-Volkov mode expansion, and the interactions relevant for Breit-Wheeler are transformed into a quantum circuit. Quantum simulations of a"null double slit"experiment are found to agree well with classical simulations following the application of various error mitigation strategies, including an asymmetric depolarization algorithm which we develop and adapt to the case of Trotterization with a time-dependent Hamiltonian. We also discuss longer-term goals for the quantum simulation of SFQED.


Introduction
Quantum computing technology is rapidly developing.At present, real quantum devices do not offer useful advantage over classical computers, but in the not terribly distant future, quantum computing may mature into a valuable tool with diverse applications.In high energy physics (HEP), it is hoped that quantum computers will be able to simulate phenomena in models where classical simulation techniques are limited by sign problems or the need to explore a large Hilbert space [1][2][3].In the near term, it is important to develop the tools to maximize the power of NISQ devices on simulation problems, and to flesh out the space of target problems for which systematic improvements in quantum computing power would lead to clear advantage over classical simulations [4,5].
The Schwinger process, e + e − pair creation in 1+1 dimensions in a strong background electric field [6], is a widely used benchmark for quantum simulations in HEP (see, e.g., [7][8][9][10][11][12][13][14][15][16][17].)As a nonperturbative, dynamical process in a confining gauge theory, it captures some of the essential physics of four-dimensional QCD-like theories in a simpler setting more tractable for NISQ-era devices.The potential applications of quantum computing to QCD are well known [18][19][20][21][22][23].Perhaps less widely appreciated are the potential applications of quantum simulation to QED, particularly in the presence of strong background electromagnetic fields (commonly referred to as strong-field quantum electrodynamics, or SFQED -for a recent review, see [24].)Indeed the Schwinger process [25][26][27] is a simple example of an SFQED phenomenon, but the richness and complexity of the theory arises in the presence of dynamical photons.Remarkably, although QED near the perturbative vacuum has been tested with unparalleled precision, there are still aspects of QED -most dramatically in the regime of ultra-strong fields -where the behavior is not fully understood (see e.g.[24,28]).Qualitatively, ultra-strong fields should yield explosive particle production and backreaction, leading rapidly to a complicated quantum state.Thus quantum computers may one day provide a unique probe of the most extreme regimes of QED.Similarly, quantum computers could provide a valuable tool to study phenomena in moderately intense fields.In this regime, which will be probed in upcoming experiments [29], analytic techniques are available for studying few-particle scattering amplitudes, at low loop order, in idealized background fields.As in QCD, quantum simulations could provide complementary access to real-time phenomena in more complex states.
These are future goals.Near-term quantum computers will be limited by size, connectivity, and noise.Therefore the present task is to design general simulation frameworks and test them on simple benchmark processes, tractable on present-day hardware.
Currently, quantum noise is unavoidable in simulations, and error mitigation is a powerful tool to "fit and subtract" various sources of noise.Fortunately, the last decade has seen rapid progress in the development of a suite of effective mitigation techniques [48][49][50][51][52][53].Here we employ measurement mitigation, Pauli twirling [50], and depolarization mitigation.For depolarization mitigation, we improve on so-called "self-mitigation" [19,54,55] by relaxing the symmetric depolarization assumption.We find that at present these techniques are essential to obtain acceptable results.Without the combination of these mitigation techniques, our quantum simulations would yield unusable data.
This paper is organized as follows.In Sec. 2 we derive the momentum-space SFQED Hamiltonian and describe the truncation approach we use for quantum simulations.In Sec. 3 we translate the nonlinear Breit-Wheeler process into quantum circuits.In Sec. 4 we describe the specific benchmark process of interest and develop a theoretical understanding of the relevant interference phenomena using time-dependent perturbation theory.We also perform real-time classical simulations in order to compare with quantum simulation results.In Sec. 5 we show raw quantum data and its comparison with Qiskit noisy simulations.In Sec. 6 we describe the error mitigation strategies used to obtain the final mitigated results, which are found to be in good agreement with the classical simulations.Finally, in Sec.7 we describe some directions for future work.

Conventions
We work in natural units with ℏ = c = 1 and report values in MeV.m = 0.511 MeV is the mass of the electron and e = 0.303 is the electric charge.
It is convenient to work in light-front coordinates [33,56,57].The coordinate fourvector is x µ = (x + , x − , x 1 , x 2 ) = (x + , x − , x ⊥ ), where x ± ≡ x 0 ± x 3 .Similar notation holds for other four-vectors.The Minkowski metric and inverse are To alleviate the subsequent bookkeeping of factors of two or one-half, we work mainly with raised indices, so that x ∓ → x ± 2 and x i → −x i (for i = 1, 2).In other words, the covariant coordinate form is ).An exception to this convention is with the four-gradient, which will be written as We take x + to be the light-front time coordinate.The spatial 3-vector and kinetic momentum are We define their inner product as px p − is light-front energy, and free particles satisfy the dispersion relation In the light-front, p + and hence p − are positive semidefinite.We will work in a Fock space constructed from a momentum lattice, and the zero mode p + = 0 does not propagate.Therefore we will take p + > 0.

SFQED Hamiltonian
Given a quantum state |ψ⟩, the generator of light-front time x + translations is P + , i∂ + |ψ⟩ = P + |ψ⟩.Equivalently, ∂ + |ψ⟩ = − i 2 P − |ψ⟩ where we identify P − ≡ H as the Hamiltonian.The Hamiltonian density is given by the Legendre transform where L is the Lagrangian density Within the covariant derivative D µ = ∂ µ + ieA µ + ieA µ we make explicit the electromagnetic gauge field A µ and the classical1 background field A µ .ψ is the Dirac field and the Dirac adjoint is defined as usual by ψ = ψ † γ 0 .We work in the gauge A + = 0 for the photon field and assume the Lorenz gauge ∂ µ A µ = 0 for the background field.For plane waves, A µ is a function of a light-front wavevector κ µ : A µ = A µ (κ µ x µ ).If the background field is a null field propagating in the −ẑ direction, κ µ = (0, 2ω, 0, 0), then the Lorenz gauge reduces to A + = 0.Here ω is the frequency of the background wave.
Under these gauge conditions the Euler-Lagrange equations for A µ and ψ reveal two constraint equations.For the electromagnetic field, we have For the Dirac field, we obtain the Dirac equation (i / D − m)ψ = 0, but defining the projectors and applying Λ − from the left, we obtain the constraint where ψ ± = Λ ± ψ and ψ + + ψ − = ψ.We will solve the constraints exactly for ψ + and A − , Note that γ 0 γ + = 2Λ − and Λ 2 − = Λ − , resulting in the appearance of ψ − in the solution for A − .
In defining the Green functions appearing in Eq. (11), it is customary to take anti-symmetric boundary conditions for fields at the longitudinal boundaries [32,34,58].As a consequence the zero mode is omitted from the spectral decomposition.
Having fixed the gauge and solved the constraints exactly, the Hilbert space we construct is a physical one, with no additional non-gauge-invariant sectors.This approach minimizes the number of qubits needed to represent the degrees of freedom, at the expense of introducing additional interactions.
Inserting Eqs.(11) into Eq.( 6) we obtain a nonlocal expression for the Hamiltonian density in terms of the dynamical fields A i and ψ − .The full Hamiltonian density can be separated into several terms: the fermion energy the photon energy the fermion-background energy the fermion-photon-background interation the fermion-photon interaction the four-fermion interaction and the double fermion-photon interaction

Momentum-Space Hamiltonian and Discretization
Most simulations of quantum field theories on quantum computers use a real-space lattice discretization.We instead use a mode decomposition and construct the corresponding Fock space in the usual way.This approach has the advantage that it makes extracting physical quantities much easier with noisy, low-qubit-number devices.In terms of asymptotic scaling, the mode expansion results in a Hamiltonian with O(n 4 ) terms, compared with O(n 3 ) terms for a spatial lattice.We define creation and annihilation operators a † /a, b † /b, and c † /c for the electron, positron, and photon, respectively.These respect the (anti)commutation relations Here, s, r = ± 1 2 are helicity indices for the fermions and λ = 1, 2 are polarization indices for the photon.The Fock states are labeled by occupation numbers and require an ordering, particularly for the fermions.We write the states in the following order: where the superscript is the helicity/polarization, the subscript is the four-momentum, and n is an index for all allowed four-momenta for a given particle and helicity/polarization.(For a given three-momentum, p − is determined by the dispersion relation.)We note that the fermion creation and annihilation operators induce a phase factor of (−1) ζ , where ζ is the sum of the number of fermions to the left of the operand as written in Eq. (19).The fermion states can be at most singly occupied, while the photon states can be arbitrarily highly occupied.The mode expansion of the fields also requires helicity bispinors for the Dirac field and polarization vectors for the electromagnetic field.Since ψ − is projected with Λ − , it is convenient to define basis bispinors for each helicity as We define linear polarization four-vectors In the free limit e = 0 the constraint equation determines ϵ − j = 2k j k + , and the photon is transversely polarized, ϵ µ j k µ = 0.With these definitions, the Schrödinger picture mode expansions are These fields can then be substituted into the Hamiltonian densities to obtain the Schrödingerpicture Hamiltonian: Subsequent integrals may be evaluated by using After constructing the Hamiltonian in this way one takes the additional step of normal-ordering the creation and annihilation operators.(See App.A for sample calculations.)This procedure renormalizes away the simplest "ear diagram" divergences, and it clarifies that the Fock vacuum is the exact ground state in the light-front.
Simulating the theory requires discretizing the momenta, truncating the set of singleparticle states, and truncating the photon state occupation numbers.We use a momentum lattice with lattice spacing 2π L .The creation and annihilation operators then obey discrete (anti)commutation relations: As a result, they are scaled by a factor of √ L 3 relative to the infinite volume case.For example, a s p → √ L 3 a s p in passing from the continuum to discretized single-particle state space.

Interaction Picture
The SFQED Schrödinger-picture Hamiltonian can be split into a "free" piece, quadratic in the fluctuating fields, and an "interaction" piece, accounting for interactions of order e and e 2 .The free Hamiltonian H 0 = H ψ +H A +H ψA is composed of number operators, so it is diagonal in the Fock basis, whereas the interaction Hamiltonian V = H ψAA +H ψA +H 4ψ +H 2ψA is not.As a result it is convenient to work in an interaction picture.Since even the Schrödingerpicture Hamiltonian is time-dependent in the present context, let us review the derivation of the interaction picture time evolution operator.
The Schrödinger-picture time evolution operator is where H S (x + ) is the complete time-dependent Schrödinger-picture Hamiltonian.The lower limit of integration is set to zero but in general denotes some initial fiducial time at which the Schrödinger-, interaction-, and Heisenberg-picture states are the same.Schrödinger-picture states are evolved with U , |ψ S (x + )⟩ = U |ψ(0)⟩, and satisfy the Schrödinger equation with H S .Now define the free evolution operator and the interaction picture states In the last line we have defined the interaction-picture time evolution operator U int = U † 0 U .It can be shown to satisfy where the second line defines the interaction-picture interaction Hamiltonian H int .(Essentially it amounts to inserting the Volkov mode solutions [59] into the Schrödinger-picture interaction Hamiltonian V -see Appendix A.) The solution to Eq. ( 29) is The Schrödinger picture time evolution operator can be written as U = U 0 U int , where U 0 is entirely diagonal and amounts to phases when acting on free-particle basis states, so the typical probability between these states can be calculated just using U int .

Quantum Circuit Design
There are many ways to build the same unitary as a quantum circuit, but in the NISQ era we must be careful to minimize the use of noisy gates, particularly CNOTs and SWAPs.At present we cannot afford to simulate the most scientifically interesting SFQED processes with many degrees of freedom.However, as a first step toward the long-term goal and as a proof of principle, we will show that a simple example process, a tree-level three-body process with exclusive final states, is reliably simulated.In this section we build three-qubit circuits that optimally simulate the parts of the SFQED Hamiltonian responsible for the tree-level nonlinear Breit-Wheeler process, shown in Fig. 1.We encode each particle in a qubit q.For the processes considered here it is sufficient to truncate boson occupation numbers at one.If q is in the |0⟩ state, then the Fock state is unoccupied; if q is in the |1⟩ state, the Fock state is occupied.For a truncated Hilbert space built from one single-particle state each for an electron, positron, and photon, the encoding is given by a bitstring of three qubits, such as |001⟩, written as a Fock state according to Eq. (19).Thus the state |001⟩ is a single photon and |110⟩ is an electron-positron pair.We follow Qiskit's little-endian notation, so that the qubits are numbered |q 2 , q 1 , q 0 ⟩.That is, the photon is qubit zero, q 0 .The terms in the Hamiltonian that describe the transition between these states are those with the operators a † b † c and c † ba.We transcribe these two operator monomials into This "decay" proceeds at tree level in SFQED in the Furry expansion [60] (see Eq. ( 77)).Applying the Jordan-Wigner transformation to represent three one-particle states by qubits, we translate the time evolution operator into a quantum circuit.
quantum gates with the Jordan-Wigner transformation.For example: where X i , Y i , and Z i are the Pauli operators acting on the ith qubit.Z gates account for the fermion ζ factor.By using Pauli identities, we can write the monomials as a linear combination of Pauli strings, each written purely with X and Y gates.In the case of three qubits, there are 2 3 = 8 different terms in the linear combination, since each qubit can be acted on with either an X gate or a Y gate.In other words, We emphasize that because the interactions in the Hamiltonian are limited to threeand four-body, the Pauli string scaling is not exponential but rather polynomial in the number of qubits.(The nonlocality in momentum space increases the connectivity and therefore the power of the polynomial; however, by solving the constraints explicitly, even a position-space lattice implementation would have higher connectivity.) Although the true time evolution operator is given in Eq. (30), this operator must be discretized as well to digitally simulate time evolution.The typical prescription is to use the Lie-Trotter approximation, also referred to as first-order Trotterization.Although higher order Trotterizations such as the second-order Suzuki-Trotter formula lead to better accuracy for longer time steps, the circuit length is increased per time step, leading to more quantum error.It is a problem-dependent question whether it is worthwhile to implement the second-order formula, and in our case, we found that first-order is preferable.Additionally, we wish to see fine-grained details not found using long time steps.Therefore we write where n denotes the number of time steps used in the approximation, so x + n = ∆x + and x + n = n∆x + .As a result, we must build a quantum circuit representing the matrix exponential of Eq. ( 32).Here again we opt for a first-order as opposed to second-order Trotterization for the same reasons discussed above, so that Each exponential factor now forms a subcircuit.Importantly, the order of these factors does not matter in the approximation, and it is very useful to shift certain subcircuits around and apply quantum gate identities to minimize the number of gates, particularly CNOTs.There are numerous ways to simulate each subcircuit.We employ the GHZ diagonalization used in previous literature [19,61].Subcircuits for Pauli strings with an odd number of X (Y ) gates can be diagonalized with a GHZ transformation circuit.The transformation circuits are shown in Fig. 2 and the diagonalization relations are shown in Tab. 1.We note that the qubit that is acted on by X-or Y -basis transformation gates is the qubit that will always have a Z gate upon diagonalization.Therefore, to reduce the use of long-ranged CNOT gates, it is ideal to place these transformation gates at a middle qubit, making the middle qubit the control qubit.Odd Xs S X (string Since S X and S Y are unitary, we may write e S † X AS X = S † X e A S X .Then by rearranging the matrix exponentials in Eq. ( 34), we can write the Trotterization with one half diagonalized by S X and the other half diagonalized by S Y .This way Eq.( 34) becomes Each matrix exponential can now be simulated solely with CNOTs and R Z gates.The four relevant circuits are shown in Fig. 3 and the final circuit is shown in Fig. 4. In practice, IBM quantum computers use a limited set of native quantum gates: X, R Z (θ), CNOT, I, and √ X.Thus, the final step in circuit creation is to translate the circuit into these gates.For example, the Hadamard gate becomes with an overall phase of e i π 4 , which falls out of the probability.
4 Nonlinear Breit-Wheeler and the null double slit experiment In the previous section we constructed a gate implementation for a truncation of the SFQED Hamiltonian.The truncation was drastic, but designed to capture the tree-level Breit-Wheeler interaction of Fig. 1 with fixed arbitrary initial and final states, and arbitrary null background fields.With this simple first-order interaction, it is possible to realize interesting quantum dynamics.An elegant example is a null version of the double-slit interference pattern [43,44].We will use this process as a benchmark for quantum simulation of the truncated SFQED Hamiltonian, computing the real-time probabilities and comparing with classical simulation.First, to build intuition, we describe the process using first order perturbation theory.Actually, although the truncated Hamiltonian operates on an eight dimensional Fock space, for the purposes of the Breit-Wheeler process, we need only consider a two-level subsystem corresponding to the basis states |e − e + ⟩ = (1, 0) and |γ⟩ = (0, 1).We are interested in the amplitude ⟨e − e + |U |γ⟩.
The simplest Hamiltonian truncation corresponds to the Hilbert subspace where the electron and positron have the same helicity with equal and opposite transverse momenta.Let the helicities be s = 1  2 , and let the photon polarization be λ = 1.We take the photon momentum to be K µ = (2p + , 0, 0, 0), the electron momentum to be P µ = (p + , p − , p 1 , 0), and the positron momentum to be Q µ = (p + , p − , −p 1 , 0).In the latter two cases the mass shell condition fixes p − = (p 1 ) 2 +m 2 p + .The background field we will consider corresponds to a linearly-polarized plane wave A µ (κ µ x µ ) = (0, 0, A 1 (κ µ x µ ), 0).For a null field traveling in the −ẑ-direction, we have κ µ x µ = ωx + .The electric field strength E can be nondimensionalized as ξ = eE mω , where e is the gauge coupling [62].A simple case to consider is a very short-duration pulse.An example is given by the gauge potential Analytical treatment becomes straightforward if we take the limit of ω → ∞, creating a delta-function pulse and a Heaviside-function potential: In the null double-slit process of [43], we allow the photon to collide with two such pulses separated in light-front time.Let one pulse appear at x + = 0 and the other appear at x + = T .Then our background field is Fig. 5 shows a spacetime diagram of this nonlinear Breit-Wheeler process.

Perturbation theory
Since the process is relatively simple, we can describe the dynamics with first-order perturbation theory.This analysis is somewhat tangent to the main purpose of the paper, for which we could make do merely with a quantum simulation and a classical exact diagonalization numerical treatment for comparison.However, it is helpful to have an analytic description in order to understand the physics of the results and the limitations of perturbation theory.
To calculate the probability |⟨e − e + |U (x + )|γ⟩| 2 we expand the interaction-picture time evolution operator to first order, The interaction-picture Hamiltonian consists of Here Using y + 0 θ(z + − T ) n dz + = (y + − T )θ(y + − T ), we can evaluate the phase to be Now let us define the parameter combinations Rewriting e aθ(y + ) = (e a − 1)θ(y + ) + 1, e i(f (P )+g(Q)) is expressed compactly as Note that θ(y + ) → 1 since we are considering y + ≥ 0. Putting together the pieces we find Integrating over light-front time, we have the probability amplitude We can find separate probabilities for x + ≷ T .If x + < T , then the probability of pairproduction is If x + > T , then the probability is Let us now inspect Eq. ( 47), the probability after the first pulse has arrived.Time-averaging over a cycle, the probability is which is maximized with respect to the background field strength when α is minimized.This corresponds to mξ = p 1 , or α = − (p 1 ) 2 p + and β = −α.Consequently, if we tune p 1 near to mξ we maximize the conversion probability.Furthermore, we note that since α is parabolic in mξ, there will always be pairs of momenta modes that reach the same level of enhancement.
Restoring the oscillating part of Eq. ( 47), we also see that the maximum probability is achieved when x + = (2n+1)π α+p − for n ∈ Z + .Let us now inspect Eq. (48).Taking β = −α and time-averaging, we find the following ratio with Eq. ( 47) To maximize this ratio, we recognize that α < 0, so we need T = (2n+1)π α+p − .If we further take the high-energy limit p 1 ≫ m, then α = −p − and we obtain This gives the factor of enhancement from the first time-averaged probability to the next.

Classical Simulation
In the near-term the accuracy of quantum computations must be benchmarked against classical simulations.We begin by configuring the momentum lattice.The lattice spacing is 2π L and to approximate the continuum limit, we need this spacing to be less than the typical momentum in the processes of interest.We will take a large L = 50π m .Let us now investigate the pair-production of electron-positron pairs with differing transverse momenta.As discussed above the electron four-momentum is given by P µ = (p + , p − , p 1 , 0), the positron four-momentum is Q µ = (q + , q − , q 1 , 0), and the photon fourmomentum K µ = (k + , 0, 0, 0).Here p + = q + , so k + = 2p + , and q 1 = −p 1 .Due to equal and opposite transverse momenta, we can classify an electron-positron pair by the produced electron transverse momentum p 1 .
For concreteness we will take ξ = 6, so that the maximum possible enhancement is seen for p 1 = 6m.In addition to simulating this momentum mode, we simulate also p 1 = 0 and p 1 = 12m modes, which should have equal enhancements to each other.p + is arbitrary and we fix it to 7.2m.Correspondingly, the light-front energies are k − = 0 for the photon, and p − + q − ≈ (0.3, 10, 40)m for p 1 = (0, 6, 12)m, respectively.Physically, the sign of the momentum results from a negative transverse electric field E ⊥ = −∂ + A ⊥ .Positive particles will experience a force in the −x 1 direction while negative particles will be accelerated in the +x 1 direction.
In addition to momentum, the Fock states possess helicity and/or polarization quantum numbers.We let the initial photon be an equal superposition of both polarizations.We will simulate the production of electron-positron pairs with all four different helicity configurations, and sum their probabilities to obtain the final pair-production probability.Accordingly, we keep two photon polarization states for a single momentum mode and four electron-positron helicity states for three electron-positron momentum modes.This creates a 14-dimensional Hilbert space.
Lastly, let us discuss state preparation.The Fock states are not eigenstates of the interacting Hamiltonian.During time evolution, Fock states evolve into superpositions of the true eigenstates.We use adiabatic state preparation to map the initial and final states back and forth between these bases, adiabatically turning on/off the electric coupling e. Adiabatic turn-on/off is performed linearly over 1000 time units, compared to the physical evolution time of 100 units during the two-pulse encounter, where the coupling is constant.However, the effects of adiabatic turn-on are negligible for the finite p 1 modes due to kinematic suppression (these modes are of high light-front energy).Since adding thousands of Trotter steps is prohibitive in the NISQ era, we will focus on the p 1 = 6m mode in the quantum simulations, for which adiabatic turn-on can be neglected2 .Furthermore we measure probabilities just after the collision with the second pulse, neglecting adiabatic turn-off, which is also a good approximation for this mode.
Fig. 6 shows the classical real-time Hamiltonian simulation.As expected, the pairproduction probability for the p 1 = 6m mode can experience constructive or destructive interference depending on the light-front time delay T separating the pulses.The p 1 = 0 and p 1 = 12m modes experience equal enhancements asymptotically.The various oscillatory behaviors appearing before and after the second pulse are captured by the perturbative formulae.After the first pulse, Eq. (47) shows that both the frequency and amplitude of oscillations are sensitive to the combination α + p − , which is quadratic in p 1 .If p 1 = mξ, then α + p − is minimized and we see relatively shorter frequencies and larger amplitudes.If p 1 ̸ = mξ, then we see longer frequency and shorter amplitude oscillations (as in the p 1 = 0 and p 1 = 12m modes).From Eq. (48) we see that the oscillation frequencies after the second pulse depend on the combination α + β + p − .Although still quadratic in p 1 , this expression is minimized and equal to m 2 /p + , independent of ξ, when p 1 = 2mξ.This is why the p 1 = 12m mode has a lower frequency than the p 1 = 0 and p 1 = 6m modes.The amplitudes of the modes now have a more complicated dependence in p 1 .Roughly speaking, a factor of 1 (α+β+p − ) 2 leads to large amplitudes for p 1 ≈ 2mξ.Away from this minimum, however, the amplitudes are suppressed by O(1/ξ 2 ).Correspondingly for the p 1 = 0 and p 1 = 6m modes we see much smaller amplitudes: the second pulse effectively "stops" the oscillations.For these modes far from 2mξ, adiabatic turn-off is therefore less significant.
The simulations allow us to visualize interference effects in real time, which are most important in the p 1 = mξ mode.The first pulse gives a kick to the p 1 = mξ mode and its probability starts to oscillate.The second pulse stabilizes the oscillation of this mode near its maximum.At the same time, the application of two kicks also yields an amplitude in the p 1 = 2mξ mode.As a side note, if we introduce a third field with strength −mξ, we might expect to obtain further constructive interference in the p 1 = mξ mode.In other words, if we vary the net field between mξ with consecutive pulses, we expect to see multiple constructive interference enhancements in the p 1 = mξ mode pair-production probabilities.Fig. 7 exhibits this qualitative behavior3 using a background field of the form The light-front time differences T odd − T even and T even − T odd were kept constant.Subsequent enhancements exhibit somewhat shorter frequency oscillations, so fine-tuning the pulse separations could lead to greater net constructive interference.

Quantum Results
Near-term quantum devices are limited by quantum noise.This means we cannot precisely perform the same classical simulation on a quantum device.The equivalent number of qubits needed for the Breit-Wheeler simulation of Fig. 6 is n = 14. 4 Based on the discussion below Eq. ( 32), we would then expect of order 3(2 × 2)(2 × 2)(2 × 2) = 192 distinct Pauli strings in the Hamiltonian.For a rough lower bound estimate, a single Trotter step could be implemented with ≳ 600 gates.Moreover, the majority of Trotter steps would require long-reaching CNOT gates which would introduce many SWAP gates and hence more noise.For this reason, we will return to the three qubit quantum simulation discussed in Sec. 3, achieved by truncating to fewer helicities and polarizations and momenta.
Beyond gate errors, the probabilities we can actually measure are limited by shot noise.This is a problem particularly for real theories characterized by a weak coupling, like QED, where interesting processes may have quite low probabilities.To avoid requiring an infeasible number of shots, we artificially boost the coupling constant e = 0.303 → 60 and decrease the lattice parameter L → 6π m .For the electron momentum, we let p + = p 1 = 8m 3 .(It still is the case that the positron momentum has q + = p + and q 1 = −p 1 .For the photon, k + = 2p + as before.)With this choices we are still effectively near the continuum limit, since the momentum lattice spacing 2π L = m 3 is an eighth of the typical momentum magnitude.Although taking e = 60 appears decidedly non-perturbative, what our simulations compute are exclusive probabilities in small cells of Fock space, or effectively differential probabilities times small phase space factors.By retaining only a few states in the full Fock space, each with a small volume ( 2π L ) 3 , the truncated quantum mechanical model is still perturbative despite the large value of e. Tree-level results in the full QFT may then be obtained by rescaling.We emphasize however that this approach will break down if many more Fock states are kept in the simulation (increasing the amplitudes of second-order and higher transitions in the truncated model.)Expanding the truncated theory in this way is of course essential to match with continuum SFQED beyond the leading order and to see more interesting dynamics, so care will have to be taken to design simulations that can operate at the physical coupling without requiring prohibitively large shot counts.In any event, for our purposes an upscaling of the coupling is sufficient.Fig. 8 shows that for e = 60 we still obtain qualitatively reasonable agreement between Trotterization and firstorder perturbation theory in the truncated quantum mechanics (or the tree level prediction of full SFQED) in this case.
The single Trotter step circuit in Fig. 4, once expressed in native IBM gates, is 40 gates long.While the classical simulation discussed previously was performed over 21000 timesteps, most of which were adiabatic turn-on/off steps, the quantum simulation must be limited to a few hundred gates total in order to avoid total decoherence.To this end, we omit adiabatic turn-on/off in the quantum simulations.Instead, we fix the initial state to the free photon state |001⟩, which as we have discussed previously is a good approximation for the physics of interest.We also begin the time evolution at the (light-front) time of collision between the photon and the first pulse and end it just after the collision with the second pulse (e.g. the quantum simulation is performed only between the dark gray lines of Fig. 6).Within this region, the pair-production probability oscillates.Due to NISQ limitations, we can only afford to capture a small portion of these oscillations, such as half a period.Using 10 Trotter steps, we show how the probability changes from minimum to maximum.To fully capture the noisy dynamics, we also record the γ → γ probabilities.
It is convenient to estimate the impact of noise first on a classical simulator, using a noise model based on historical measurements.Results using Qiskit's FakeNairobiV2 noisy simulator [63] are shown in Fig. 9a.The data suggest that we should expect at least 10% errors already in the initial timestep, increasing with time.At e = 60, perturbation theory still yields probabilities below one, and the two methods agree qualitatively.Perturbation theory provides a useful understanding of the process, despite the large coupling, because the phase space volume is small (both for the final state and the Hilbert space truncation as a whole, which prevents higher order effects from becoming large.)Fig. 9b shows raw quantum data from ibm_nairobi [64].The real noise level we observe is higher than seen in the noisy FakeNairobiV2 simulator.At x + = 4 the error in the γ → γ probability is already 60%.Such results are not uncommon from NISQ devices and can be substantially cleaned by mitigation algorithms.

Error Mitigation
Mitigating against specific sources of device noise is critical to extract useful information from near-term quantum simulations.We employ measurement mitigation and Pauli-twirling [50], as well as a modified version of self-mitigation [19,54,55] (mitigation of depolarization noise).The combination of these three error mitigation techniques significantly improved the quantum data as shown in Fig. 13b.For completeness we discuss each mitigation algorithm in some detail.

Measurement Mitigation
With measurement mitigation we account for readout errors of bitstrings.For example, there is the probability that upon readout, the measured bitstring is the computational basis state |000⟩ as opposed to |001⟩.For n qubits, the mitigation strategy is to run 2 n circuits that  prepare individual computational basis states and derive statistics for how many counts there are for the 2 n possible results.
Let {|ψ i ⟩} for i = 0, 1, . . ., 2 n − 1 be the computational basis.Prepare circuits C 0 , C 1 , . . ., C 2 n −1 that prepare each computational basis state as in Fig. 10.For each circuit q 0 X q 1 q 2 For a general quantum circuit, let #» c true be the column vector of noiseless measurement counts of each basis state.The product C #» c true distributes the true measurement counts among other basis states such that the calibration matrix models the effects of readout noise.Thus, define #» c meas = C #» c true .
In reality, noisy quantum measurements yield #» c meas , so in principle we can obtain the mitigated counts by matrix inversion: #» c true = C −1 #» c meas .However, simply inverting C can lead to negative counts and counts that do not sum to the number of shots N .Therefore, we utilize a least-squares minimization protocol.This method takes advantage of the notion that #» c meas − C #» c true = 0. Let #» x be an estimate for #» c true .Then the problem becomes finding subject to the constraints that 0 ≤ x i ≤ N and | #» x | = N .The minimization will yield the best estimate of #» c true as predicted by C. The mitigated probabilities are then 1 N #» c true .

Pauli Twirling
CNOT gates may have different errors depending on the state that is fed to them.Pauli twirling implements gate identities around CNOTs in order to feed them different states, thus "symmetrizing" the noise.Fig. 11 shows examples of CNOT identities.In total there are 16 ways to twirl a CNOT gate with extra Pauli gates.We randomly twirl all CNOTs with different identities in each run, then average over the twirled results in post-processing.

Depolarization Mitigation
The quantum channel often used to describe depolarizing noise [65] is Figure 11: Sample Pauli twirls of the CNOT gate.For error mitigation, each circuit was randomly compiled with every CNOT twirled using these equivalent gates.
Under this channel a state ρ may turn into the fully decohered state I 2 n with probability p, and with probability 1 − p the state remains unchanged.This type of depolarizing noise model is called symmetric depolarization and is the basis for many implementations of "selfmitigation" [19,54,55].A more general noise model could take ρ not just to I 2 n , but to other general mixed states, which is what we will now consider.
We generalize the depolarizing noise model by introducing a transfer matrix M : where ) is the true (noisy measurement) probability of computational basis state |ψ i ⟩.This implies the following quantum channel: Symmetric depolarizing noise is the case where M is a one-parameter matrix with entries M ij = δ ij (1 − 2 n ϵ) + ϵ. ϵ may be estimated by preparing circuits with similar noise elements as the circuit of interest.Likewise, to estimate the elements of M ij in our model, we construct circuits C 0 , C 1 , . . ., C 2 n −1 , similar to those for measurement mitigation.However, in addition to preparing each computational basis state, the circuit applies an identity operator with the same gates and the physics circuit, as shown in Fig. 12.Each noisy identity is The depolarizing mitigation circuit C 1 (x + 2 ) for the second time step, preparing |ψ 1 ⟩ = |001⟩ when n = 3 qubits.The noisy identity is the time evolution circuit used for the quantum simulation followed by its Hermitian conjugate.
the current time evolution operator followed by its inverse.This means that the mitigation circuits C i are time-dependent C i (x + ), allowing us to mitigate against time-dependent noise.The transfer matrix elements are then estimated as We emphasize that in both the symmetric case and our case, Pauli-twirling is essential to convert more general types of noise from the CNOTs into depolarization noise.Thus, depolarization mitigation is the last step of the error mitigation pipeline.Both the physics and mitigation circuits are twirled and the results are averaged before depolarization mitigation.
Furthermore, since the inserted identity operators are effectively "twice the physics circuit," each mitigation circuit is practically a measure for the noise "squared".This implies that we should be taking √ M as the calibration matrix. 5That is, As before, to determine #» c true , we solve an optimization problem subject to the same constraints on #» x as before.The weights W are taken to be W = diag( 1 We checked that symmetric depolarization noise mitigation following Pauli twirling did not perform well on our problem: the noise profile indicated some states "thermalize" with each other, but not all together.On the other hand, the construction of a complete calibration matrix for depolarization mitigation is computationally expensive and does not scale well with problem size.For the small size studied here it was not critical to optimize it further, but it would be interesting to explore the performance restricted classes of calibration transformations constructed from fewer parameters and mitigation circuit measurements.

Mitigated Results
For each time step, we performed 10 Pauli-twirls on the quantum simulation circuit as well as each of the 2 3 = 8 mitigation circuits.Each time step was sent as a different set of circuits to ibm_nairobi at different times, so in addition we ran 8 measurement mitigation circuits each time step.Therefore, 98 circuits were sent to ibm_nairobi at each time step.
As before we can compare with Qiskit's noisy FakeNairobiV2 simulator.Fig. 13a shows the mitigated noisy simulator results.The mitigation is substantial and suggests that the real quantum simulation should also see improved results vs the raw data.
Fig. 13b shows mitigated quantum data.As expected, the error mitigation protocol greatly improves the simulation results.The largest percent error we see is about 15% of the classical simulation.However, substantial errors do not occur until the last two time steps.This suggests that had we performed a longer or finer time evolution, decoherence may have been too strong to be well-mitigated.

Conclusion
Quantum field theory in the presence of strong background fields is a rich arena where quantum computers may provide valuable new access.We have seen that for a simple strong-field QED process -nonlinear Breit-Wheeler pair-production in a sequence of laser pulses -the noisy quantum computers of today are already able to simulate the dynamics reliably with the help of error mitigation routines.Due to the polynomial scaling in circuit depth, the simulation of SFQED with many particles and strong backreaction presents an interesting, plausible long-term goal.
Capturing more complex electromagnetic cascade processes requires retaining a much larger Hilbert space.With a less drastic truncation one could also study helicity and polarization effects in cascades, as well as begin to tackle the regime of ultra-strong fields, where few other tools exist.The light-front Fock space formulation used in this work can be applied directly to more complicated processes with bigger Hilbert spaces, and the operator-circuit map can be straightforwardly extended.
It is also possible to simulate O(30-50) qubit quantum computers with state-of-the-art GPUs and supercomputers [66].Larger-scale SFQED quantum circuits could be designed on these simulators with adjustable noise.This capability could be useful to develop efficient circuits for future quantum computers, and to establish noise targets for interesting physics processes.Investigations in this direction are ongoing.
(Here repeated indices, although all raised, are understood to be summed over.)It is convenient to define the following coefficients R: R s,j ≡ w s † γ j w −s R + 1 2 ,j = 1 j = 1 i j = 2 R − 1 2 ,j = −1 j = 1 i j = 2 (66) Integrating over light-front space we obtain

Figure 1 :
Figure 1: Nonlinear Breit-Wheeler pair-production.Double lines indicate e + e − propagation on the background field.This "decay" proceeds at tree level in SFQED in the Furry expansion[60] (see Eq. (77)).Applying the Jordan-Wigner transformation to represent three one-particle states by qubits, we translate the time evolution operator into a quantum circuit.

Figure 2 :
Figure 2: GHZ diagonalization circuits for Pauli strings with an odd number of Xs (left) or odd number of Y s (right).Diagonalization facilitates implementation of the time evolution subcircuits with minimal gate counts.

Figure 4 :
Figure 4: First-order Trotterization circuit for a single time evolution step.By alternating the S X and S Y halves of the circuit in subsequent time steps, cancellations can be used to minimize the gate count.The first step has 16 CNOT gates, but only 12 CNOT gates are added with each succeeding step.

Figure 5 :
Figure 5: Two delta-function laser pulses propagate in the −ẑ direction, separated in time.A photon collides head-on with the first pulse, leading to a superposition of the |e − e + ⟩ and |γ⟩ states.Collision with the second pulse may cause attenuation or enhancement of the |e − e + ⟩ state, depending on the pulse separation.Sketched is a case of total destructive interference in pairproduction probability for a given electron-positron momentum mode.Light-front coordinates are convenient because all particles interact with the laser pulses at the same light-front time x + .

Figure 6 :
Figure 6: Classical simulations of nonlinear Breit-Wheeler pair-production differential probabilities as a function of light-front time.The three curves represent three different electron-positron pairs, labeled by the p 1 electron transverse momentum.The vertical light gray lines enclose the physical evolution period when the coupling is held constant.Each dark gray line represents a collision of the incoming photon with the background laser pulse.Only the p 1 = mξ mode sees significant constructive or destructive interference.Top: Two delta pulses with ξ = 6, spaced to show constructive interference in the production e + e − pairs with p 1 = 6m.Bottom: The same, but spaced to show destructive interference.(In the quantum simulations shown below, evolution is performed only between the dark gray lines.)

Figure 7 :
Figure 7: Eight consecutive delta pulses with ξ = 6 are represented by dark gray vertical lines.The signs of the pulses appear in the order + + − − + + −−.A careful selection of pulse spacings can lead to substantial enhancement of the pair production probability in specific momentum modes.

Figure 8 :
Figure 8: Comparison of real-time pair-production probabilities computed with perturbation theory and Trotterization for artificially large coupling.An initial delta pulse is encountered at time zero, and then the probability oscillation amplitudes change when the second pulse is encountered.At e = 60, perturbation theory still yields probabilities below one, and the two methods agree qualitatively.Perturbation theory provides a useful understanding of the process, despite the large coupling, because the phase space volume is small (both for the final state and the Hilbert space truncation as a whole, which prevents higher order effects from becoming large.) (a) Results from Qiskit's noisy simulator reflect modest errors, more significant in γ → γ than γ → e − + e + probabilities.(b) Data from IBM's quantum hardware reflects strong noise effects, particularly in γ → γ probabilities.

Figure 9 :
Figure 9: Unmitigated quantum results for nonlinear Breit-Wheeler pair-production.Pulses with mξ = 8m 3 are encountered at x + = 1 and x + = 10.The coupling is held fixed at e = 60.The electron transverse momentum mode is p 1 = mξ.
(a) Mitigated data from Qiskit's fake simulator shows excellent agreement with classical simulation results.(b) Mitigated data from IBM's quantum hardware agrees well with classical results until x + = 9, suppressing high levels of errors seen in raw data.

Figure 13 :
Figure 13: Mitigated quantum results for nonlinear Breit-Wheeler pair-production.Pulses with mξ = 8m 3 are encountered at x + = 1 and x + = 10.Robust measurement and depolarization mitigation techniques suppress errors at early time steps, but become less effective at late time steps.

Table 1 :
Diagonalizations of Pauli strings consisting of X and Y operators by the operators of Fig. 2.