A quantum algorithm for high energy physics simulations

Particles produced in high energy collisions that are charged under one of the fundamental forces will radiate proportionally to their charge, such as photon radiation from electrons in quantum electrodynamics. At sufficiently high energies, this radiation pattern is enhanced collinear to the initiating particle, resulting in a complex, many-body quantum system. Classical Markov Chain Monte Carlo simulation approaches work well to capture many of the salient features of the shower of radiation, but cannot capture all quantum effects. We show how quantum algorithms are well-suited for describing the quantum properties of final state radiation. In particular, we develop a polynomial time quantum final state shower that accurately models the effects of intermediate spin states similar to those present in high energy electroweak showers. The algorithm is explicitly demonstrated for a simplified quantum field theory on a quantum computer.

Particles produced in high energy collisions that are charged under one of the fundamental forces will radiate proportionally to their charge, such as photon radiation from electrons in quantum electrodynamics. At sufficiently high energies, this radiation pattern is enhanced collinear to the initiating particle, resulting in a complex, many-body quantum system. Classical Markov Chain Monte Carlo simulation approaches work well to capture many of the salient features of the shower of radiation, but cannot capture all quantum effects. We show how quantum algorithms are wellsuited for describing the quantum properties of final state radiation. In particular, we develop a polynomial time quantum final state shower that accurately models the effects of intermediate spin states similar to those present in high energy electroweak showers. The algorithm is explicitly demonstrated for a simplified quantum field theory on a quantum computer.
While quantum computers hold great promise for efficiently solving classical problems such as querying databases [1] or factoring integers into primes [2], their most natural application is to describe inherently quantum physical systems [3]. The most direct connection between quantum systems and quantum computers occurs for analog circuits that try to mimic the evolution of a Hamiltonian as closely as possible [4]. Analog circuits for a generic quantum field theory are impractical, as the closest analogs to high energy scattering processes are other high energy scattering process.
A promising alternative to analog circuits are digital quantum circuits, which use quantum algorithms to describe inherently quantum physical systems without directly implementing the system's Hamiltonian. However, many physical systems are too complex or have too many degrees of freedom to fully model with a quantum circuit using near-future noisy intermediate scale devices [5]. This is true for a generic quantum field theory, where there are both continuous quantum numbers as well as an infinite number of degrees of freedom. While tools have been developed to model quantum field theories by discretizing spacetime [6] and even including continuous quantum numbers [7], the number of quantum bits (or their continuous analog) required to compute any relevant scattering amplitude is impractically large. Results with simplified quantum field theories on a lattice are promising [8], but the full dynamics of high energy scattering processes are too complex for both lattice methods as well as traditional perturbation theory when the number of final state particles becomes too large.
A successful classical approach for simulating these dynamics is known as the parton shower [9], which relies on reorganizing the traditional perturbative series about a fundamental coupling constant to instead expand around the collinear and soft limit of emissions [10,11]. This leads to different series expansions where each term in the new series includes infinitely many terms from the original series expansion and is the basis of parton shower Monte Carlo (MC) programs [12][13][14][15], which are a key component of high energy quark and gluon scattering simulations.
Parton shower models are implemented using classical Markov Chain MC (MCMC) algorithms to efficiently generate high multiplicity radiation patterns. This reliance on classical MCMC algorithms implies that several quantum interference effects need to be neglected. For example, showers describing emissions in the strong interaction can only be implemented in the limit of a large number of colors (N C = 3 → ∞). While an impressive research effort to include subleading color effects exists [16][17][18], there is a fundamental limitation in the ability of MCMC methods to efficiently capture these physical effects. For showers describing the electroweak interactions [19], interference effects can arise because physically distinct particles can have related interactions, such that amplitudes which differ in their intermediate particles can interfere with one another. An important example is the interference of amplitudes involving intermediate Z bosons and photons.
Our primary motivation is to develop a quantum circuit describing the quantum properties of parton showers. In this work, we consider showers with interferences from different intermediate particles, using a simplified model that captures these effects without having to introduce the full complexity of the Standard Model (SM). The variable describing the scale of the shower evolution is discretized and at each step an emission can occur or not. We will show that a classical MCMC is not able to capture the important quantum interference effects in this model, and that a full classical calculation scales ex-ponentially with the number of steps 1 . The proposed quantum algorithm will be able to sample from the full probability distribution in polynomial time. After describing the physics of the simplified model, we will introduce the quantum circuit and show empirical results with both a simulated and a real quantum computer. The article ends with outlook towards a full SM quantum parton shower algorithm.
To begin, consider a simple quantum field theory, with two types of fermion fields, f 1 and f 2 , interacting with one scalar boson φ governed by the following Lagrangian: In such a theory, the scalar field φ can couple to either fermion via the coupling constants g 1 or g 2 , or to one fermion of each type, with coupling constant g 12 . The collinear dynamics of the theory are that fermions can radiate scalars (f → f φ) and scalars can split into pairs of fermions (φ → ff ). The couplings of fermions to scalar bosons occur in the Higgs sector of the SM, and it has been demonstrated that the final state collinear radiation at high energy can be written in terms of a parton shower [19,24].
We will revisit the connection to the full SM later; for now, we will consider final state radiation governed by Eq. (1) with generic couplings. This model can contain important interference effects when all couplings are non-zero, since the unobserved intermediate state of the fermions can be a superposition of f i for i ∈ {1, 2}. The observable final state is a set of fermions and bosons with their corresponding energies and locations inside the 'jet' of particles. Ignoring the φ → ff splitting for now, the jet is specified by the number and kinematic properties of the emitted bosons. For the amplitude there are n − 1 internal fermions and thus a total of 2 n−1 unobservable configurations. For example, to leading order in the coupling constants, where p i is the momentum of boson i andÂ n (p 1 , . . . , p n ) denotes the contribution to the amplitude which does not involve the coupling constants and is the same for all configurations. Already with two emissions, there is a nontrivial interference between the f 1 swapping twice and not swapping at all. For the emission of more bosons, the combinatorics required to obtain the coupling constant rapidly becomes combinatorially hard. Note that in the limit that g 12 → 0, A i→i n is simply ∝ g n i . To modelÂ n , one needs to choose a physical scale to order emissions. One common choice is to evolve based on the angle of emissions θ, where the emissions are ordered down to a collinear cutoff θ = > 0 below which further emissions cannot be resolved. In the strongly ordered limit that applies to parton showers, θ 0 θ 1 · · · θ n , the kinematic part of the amplitude factorizes as followŝ A n (θ 1 , . . . , θ n ) =Â(θ 0 |θ 1 )Â(θ 2 |θ 1 ) . . .Â(θ n |θ n−1 ) , (5) whereÂ(θ n |θ n−1 ) denotes the amplitude to emit one particle at angle θ n given the angle of the previous emission. In general, AÂ also depends on the momentum fraction of the emitted particle; in our setup, this is completely factorized from the angular dependence and henceforth ignored.
In the limit g 12 → 0, Eq. (5) allows for an efficient MCMC method for calculating high-multiplicity cross sections. This is performed by introducing two splitting functions and a no-branching probability (Sudakov factor) whereP (θ) encodes the scale-dependence of the emission probability. The Sudakov factor encapsulates the virtual (and unresolved real) contributions and is responsible for the resummation mentioned above. The Sudakov factor and splitting function satisfy the unitarity relation A classical parton shower would predict the cross-sections to be One can efficiently sample from Eq. (9) using a Markov Chain algorithm by generating one emission at a time, conditioned on the last emission. While this will correctly reproduce the physics of a theory with g 12 = 0, it does not reproduce the interference arising in the full theory given by Eq. (1) (still excluding φ → ff ), where the fermion can change in the emission. The resulting interference effects can only be included by working with the amplitudes directly. A single emission that changes the type of fermion can be treated using a density matrix formalism [19], where each splitting function is represented through a splitting matrix as In the limit of g 12 → 0 we have P i→jφ (θ) → δ i,j g 2 iP (θ), but for non-zero g 12 the full matrix structure of the splitting function needs to be retained. The complexity of taking this into account to all orders, reduces to the full amplitude calculation.
In what follows, we construct a quantum algorithm to sample from the full amplitude, including all interference effects. We consider the complete case, including φ → ff , which still follows the Markov Chain of amplitudes in Eq. (5). The core idea of the quantum algorithm is to encode the particles as qubits (Appendix A) and first rotate to a particle basis where there is no mixing between fermion states (Appendix B). In this superposition basis, emissions between states are uncorrelated. Sudakov factors can then be used to govern the no emission probability of the uncorrelated fermions. The bulk of the quantum circuitry will then be dedicated to book-keeping, to encode the emission history and decide which fermions/bosons radiate/split at a given step in the shower. Figure 1 is the quantum circuit implementing the quantum final state radiation algorithm for one of N steps. The circuit calls for six registers, which are are detailed in Appendix A and summarized in Table I. The initial state consists of n I particles (which can be fermions or bosons) in the f 1/2 basis. One starts by rotating this initial particle state from the f 1/2 basis to the f a/b basis, using a simple unitary R operation discussed in Appendix B. Then, a series of operations evolving the particles states are applied: the number of particles of each type are counted (U count ), Sudakov factors are used to determine if an emission occurred (U e ), given an emission, a particular particle is chosen to radiate/branch (U h ), and the resulting particle state is updated (U (m) p ). Finally, the state is rotated back to the f 1/2 basis through the R † operation. This process is repeated for all of the N steps. The rotation needs to be performed separately at each step because in general the matrix R depends on θ through the running of the couplings. At each step, there are four operations, which are summarized in Table II. More details can be found in the appendices.
Performing the evolution in the f a/b basis and then rotating to the f 1/2 basis, creates interferences between equivalent final states which had different intermediate fermions. One event is generated by measuring all of the qubits after the final rotation back to the f 1/2 basis. By repeating the entire process, we can generate a large number of events which we can then use to compute physical observables for our theory. The number of standard quantum gates (single qubit and CNOT gates) required at each step is discussed in Appendix I and summarized in Table II.   The practical challenge with above circuit is that it requires more connected qubits and operations than are currently available in state-of-the-art hardware. In order to show an implementation of our algorithm, we therefore consider a special case that is amenable to measurement on existing technology. This special case ignores the φ → ff splitting (naturally suppressed in gauge theories, but not in the scalar-only theory), ignores the running coupling, and has only a single fermion (possibly in a superposition) as the initial state. This results in a much simpler circuit since there is only one fermion, but an arbitrary number of scalars (Appendix I). A decomposition of the resulting circuit into single qubit and CNOT gates requires n gates = 12 N + 2 (Appendix G). This model is however still sufficiently complex that the classical MCMC described earlier 2 fails to capture important quantum effects when g 12 = 0. Figure 2 presents the normalized differential cross sections of four examples from a class of observables, i θ α i , for both classical simulations/calculations, quantum simulators [25], and chip experiments of public and Hub Analytical (g 12 = 0) 24 step simulation (g 12 = 0) 4 step simulation (g 12 = 0) 24 step simulation (g 12 = 1) 4 step simulation (g 12 = 1) 4 steps IBMQ Tenerife (g 12 = 0) 4 steps IBMQ Tenerife (g 12 = 1) 2 step sim. with ff (g 12 = 0) 2 step sim. with ff (g 12 = 1) 14 24 step Classical MCMC 24 step simulation (g 12 = 0) 4 step simulation (g 12 = 0) 24 step simulation (g 12 = 1) 4 step simulation (g 12 = 1) 4 steps IBMQ Tenerife (g 12 = 0) 4 steps IBMQ Tenerife (g 12 = 1) 2 step sim. with ff (g 12 = 0) 2 step sim. with ff (g 12 = 1) Member quantum chips through cloud access on the IBM Quantum Experience. All cases are started from the initial state containing a single f 1 fermion. The data of experimental measurements shown in Figure 2 were collected on the IBM Q 5 Tenerife chip. This quantum computer has five qubits, so N = 4 is the maximum number of steps that can be modeled. In addition to presenting the simplified model with both quantum hardware and simulations, Figure 2 also shows a simulation with the full model (including φ → ff ) for 2 steps.
For the top left plot, the histogram of the 4 step quantum simulation agrees exactly with the 24 step simula-  tion. For the other plots, the observable is introducing a dependence on the number of steps, which disappears as N → ∞. When interference effects are turned off (g 12 = 0), we find excellent agreement for all observables between both the classical and quantum simulations as well as the quantum computer measurements. For g 12 = 1 the spectra are harder, leading to more emissions and at larger angles. For all quantum simulations the fraction of events with no emissions (first bin in the bottom right plot) agree separately for each value of g 12 . This is because the simulation is started with a single fermion state, where the splitting φ → ff is irrelevant. For a higher number of emissions, the φ → ff splitting affects the distribution, and in particular lowers the fraction of events with a single emission. The data points for the quantum computer are much closer together for g 12 = 0 and g 12 = 1 compared with the simulations. Additional simulations with a noise model indicate that this is likely due to the noise present in the existing hardware. There are clear differences between the g 12 = 0 and g 12 = 1 case also for the full model, but the granularity is not yet sufficient to resolve the structure of these differences.
With improved quantum hardware beyond the current noisy intermediate scale devices [5], our algorithms will be able to produce calculations that are currently not possible with classical devices. On the algorithmic side, near term development will include modifying Eq. 1 to be the electroweak sector of the SM. Capturing the full interference effects of collinear electroweak radiation already has important practical implications for ultra-high energy dark matter indirect detection, where PeV-scale particles charged under the electroweak force would radiate copious electroweak bosons not captured properly by classical MCMC methods [26]. Circuit simplification and quantum error correction may result in useful results for near-term hardware. The further goal of this line of research is to also capture interference and entanglement from soft radiation, which is the key challenge in modeling strong force final state radiation. There has been significant progress with classical algorithms to include some of these effects as corrections [27][28][29][30][31][32], but a complete calculation may only be possible with pure quantum or quantum-classical hybrid algorithms. The richness of quantum phenomena in high energy physics makes them an excellent testbed for studying the power of quantum algorithms. By focusing on final state radiation, quantum algorithms may be able to provide key insight into the dynamics of quantum field theories underlying the laws of nature.

ACKNOWLEDGMENTS
We thank A. Ryd for useful discussions on the treatment of spin correlations in classical MCMC programs. This work is supported by the DOE under contract DE-AC02-05CH11231. In particular, support comes from Quantum Information Science Enabled Discovery (Quan-tISED) for High Energy Physics (KA2401032). We acknowledge access to quantum chips and simulators through the IBM Quantum Experience and Q Hub Network.
Appendix A: The registers of the quantum circuit The quantum circuit introduced in this paper has a total of 6 registers. The first two registers are physical registers, holding the information created by the circuit. The final 4 registers are work registers, which means that they are reset to their original value after each step. Thus they hold no information after the circuit has been run, and the same work registers can be used for each step. As discussed in other appendices, additional work qubits will be necessary when actually implementing some of the more involved circuit operations.
The first register, |p , contains the flavor information about each particle. Each particle in the system can be in one of 6 states |0 , |φ , f a/b , and f a/b . To encode these 6 states one requires 3 qubits, and we choose the representation as where the third and fourth states are not used and one chooses f 1/2 and f a/b before and after the basis change discussed in Appendix B, respectively. Since there can be up to N + n I particles in the system (where n I is the initial number of particles and N is the number of steps), one needs a total of dim[|p ] = 3(N + n I ) (A2) qubits to encode this register. The second register, |h , holds the information about which particle emitted a particle at a given step. At the start of the m th step (where the first step has m = 0), there are up to m + n I particles that can have emitted the extra particle, and at the m th step |h m needs to be able to hold the integers 0 . . . m + n I (where 0 denotes no particle having emitted something). When considering N steps, the register therefore needs to hold qubits.
The third register, |e temporarily holds the information whether an emission has occurred in the current step. This is binary information, and therefore requires a single qubit, giving dim[|e ] = 1 . (A5) The remaining three registers are count registers, which temporarily hold the information about how many bosons, fermions of type a and fermions of type b (counting both f andf ) are in the current state. Since the count registers are used for every step, they have to hold the integers 0, . . . , N + n I . We again choose the binary representation to hold these integers, and one needs qubits.
The summary of these registers was already shown in Table I.
At the start of the circuit, all work registers |e , |n φ , |n a , and |n b are initialized to |0 , where for the count registers |0 refers to the integer 0 in binary notation. For the physical registers, all history registers |h m as well as the particle registers |p m>n I are initialized to zero. The only non-zero registers are |p m≤n I , which are initialized to the initial particle content (possibly in a superposition).

Appendix B: Diagonalizing the splitting matrix
In this appendix we discuss the rotation required to go from the basis with fermions f 1/2 to a new basis with f a/b . The splitting matrix in Eq. (10) can be written in terms of the coupling constants g 1 , g 2 and g 12 as The coupling matrix G can be diagonalized as with where The matrix U in Eq. (B2) is given by with u = (g 1 − g 2 + g ) 2g . (B6) Since each particle is represented by a 3-qubit state, the operation R that rotates a single particle from the f 1/2 basis to the f a/b basis is represented by a 8 × 8 unitary matrix R, and it must be applied to all of these 3-qubit particle states. It is defined in terms of the matrix U , introduced in Eq. (B5). For the representation of the particles given in Appendix A, one has where I denotes the 2×2 identity matrix. The rotation R correctly mixes the fermion states, while it leaves alone the |φ and |0 states. Because of the running of the coupling constants the matrix U , and in turn the matrix R, will be different at each step in the evolution.
Appendix C: Populating the register for counting the particles In this section we give details about the operation which counts the number of each particle type in the current state and stores these numbers in the count registers |n φ , |n a and |n b . As discussed, at the beginning of each step they are in the state |0 . To perform this counting we apply the controlled U (m) count gate, which is broken down in Figure 3. For each particle in the state |p we apply the unitary operation U + to the appropriate count register. The operation U + is defined on a set of integer states ranging from 0 . . . N + n I as or in matrix form, (U + ) ij = 1 if j = i + 1 mod (N + n I ) and 0 otherwise. This is a simple operation, and the gate decomposition of the U + operator can be found in Appendix G.
The circuit operation for counting the particles.

Appendix D: Sudakov factors in the quantum circuit
In this section we discuss how we implement the second operation required in each step, which decides whether an emission happens or not. In the a/b basis the splitting can not change the flavor of the emitting fermion, and the evolution can therefore be described in terms of individual splitting functions and Sudakov factors, just as in a usual MCMC. For the fermions there are 2 different splitting functions where i ∈ {a, b}. The splitting of the bosons is given by Using these splitting functions, one can define Sudakov factors, which describe the probability to have no emission from a given particle in a given step m. One finds where and The probability to have no emission from a state containing n φ bosons and n a/b fermions of type a/b, is then given by From this one can derive the probability to have a branching at a given step, which is given by One therefore finds the unitarity condition This splitting probability can be encoded in the quantum circuit through the rotation U Appendix E: Selecting a particle to radiate or split The second operation discussed in the previous appendix only decided if an emission happened or not, but did not have any information about which of the existing particles the emission originated from. This decision happens in the third operation, which creates the emission history by deciding if an emission happened, which particle might have emitted and assigning the correct amplitude for the given splitting. In order to select the histories one needs to "loop" over all particles in the register, up to the (m + n I ) th particle for step m, which can be written in terms of sub-operations as shown in  in the particle register, and the final operation ensures that the emission qubit is back in the |0 state after the operation. The sub-operation which controls on the k th particle is shown in Figure 6. Here a control on |p k means that the controlled unitary operation U h depends on the flavor of particle p k , while φ k , a k and b k are true or false if |p k is either a boson, an a-fermion, or a b-fermion, respectively. 6: Circuit for the k th sub-operation from the second operation in each step.
For each particle in |p we apply U if the emission has occurred in the given step. U (m,k) h is a 2 × 2 unitary sub-matrix which always acts between the states |0 and |k of |h m . Defining where P a , P b and P φ are given coefficients, the mentioned 2 × 2 submatrix is given by The coefficients of the matrix U (m,k) h depend on n φ , n a and n b , which is why we control on the count registers. Thus, if an emission has occurred, in each sub-operation the controlled U (m,k) h gate rotates between the states |0 and |k in the |h register. This is done recursively in a way that builds up the correct amplitudes for each possible emission history.
After each application of U (m,k) h , the count register is reduced, changing the value of P (n φ , n a , n b )(θ m ) in the next step. For example, if it was a f a which emitted, the count n a will go to to n a − 1. This means in particular that in the last sub-operation one has such that the last of the 2 × 2 sub-matrix is always of form As a result, in the last sub-operation the amplitude of the |0 state of |h m is fully transferred to the |m state.
In the history register, this operation generates a superposition of states corresponding to all the possible emissions which could have happened. The amplitude for the emission to be associated with a particle p k is given by but this procedure includes the interference from all possible flavors each particle can have. Notice that at the end of the step, the U − gates have been applied conditionally on all of the particles in |p , which is exactly the inverse of the first operation, where we counted the particles. As a result the three count registers will be back to the initial state |0 at the end of each step, ready to be used again. Furthermore, the emission register is reset back to |0 by the last controlled X operation.
Appendix F: Adjusting the particle state The final operation in each step adjusts the particle flavors according to the emissions that happened. For example, if a boson splits into a f afa pair, we must remove a φ from |p and add an f a and anf a . In general, if it is a fermion that emits we simply have to add a boson to the |p register, while if it is a boson which emits we must add fermion-antifermion pair to |p and remove the boson which emitted. The schematic of the circuit which performs this operation is shown in Figure  7. Each sub-register in |p , made up of three qubits, cor- 7: Circuit for the operation at step m which fixes the particle register after the emission has happened. As before, M = m + nI − 1. Notice that if we control on |h being in the |k state, we apply Up to the k th sub-register |p k and the (M + 1) th sub-register |p M +1 .
responds to one particle state and can be in any of six possible states: 0, φ, f a , f b ,f a andf b . The sub-register |p M +1 will encode the new particle which has just been emitted and it always starts out in the 0 state, while the registers below encode the previous particle states.
The operation labeled U p , conditional on the k th state in |h , is a map from the k th and (M + 1) th particle states before the emission, and the same particle states after the emission. Notice that this operation is controlled on the history states, which specify which particle emitted, though they do not hold the information of what kind of particle that was. That information is provided by the k th particle state. The U p gate is always the same and we want it to take .

(F2)
Here we used the vector representation of the particle states given in Eq. (A1). We can write this transformation as a single unitary operator as follows: Since the particle states of different flavors are orthogonal to one another, this transformation is unitary. The circuit decomposition of Eq. (F3) can be found in Appendix G.

Appendix G: Circuit Decomposition
We now explain in some detail how to break down the operations in our general quantum circuit from Figure 1 (including φ → ff and the running coupling) into standard quantum gates (single qubit gates and CNOT gates), so that we can run the circuit on a simulator and eventually on an actual testbed. While every effort was made to find an efficient breakdown of the circuit, we anticipate that a reduction in the number of standard quantum gates is still possible. The following discussion gives the number of gates required for each step 0 ≤ m < N − 1 in the evolution. Table II gives the number of gates needed after summing over all steps.

The first sub-operation, Ucount
We start with the counting operation shown in Figure 3. We store integers in the counting registers using the conventional bit representation, then the U + gate (see Appendix C) can be implemented as shown in Figure  8. A general integer a has the form |q ...q 3 q 2 q 1 , where = log 2 (a) is the number of bits necessary to store the integer (we round up to the nearest integer). Therefore, in our circuit the number of gates needed to implement a specific U + gate depends on the maximum integer we might have to store. As shown in Figure 3, the U + gate is  controlled on the particle state |p being a type a fermion, a type b fermion or a φ boson. As illustrated in Fig. 9, the first two cases require controlling on two qubits from |p , while the latter case requires controlling on all three the qubits from |p . These controls are applied to all of the operations shown in Figure 8, yielding many instances of an X-gate being controlled on multiple qubits. It is a known result (see e.g. Ref. [33]) how to decompose a C (n) (U ) operation, requiring n−1 work qubits, 2×(n−1) Toffoli gates plus a C(U ) operation. A Toffoli gate requires 16 standard gates while a C(U ) operation where U is real requires 5 standard gates in general (although if U = X it is simply a CNOT gate). For n > 2 and controlling on all qubits being in the |1 state, we then need standard gates. To this count we add 2 X-gates for each time we control on a qubit being in the |0 state instead of the |1 state. Using these results, the total number of standard gates necessary for the counting operation when simulating the m th step is: The above number includes many pairs of adjacent X gates (coming from controlling on a |0 , rather than |1 ) that cancel. Ignoring all such X gates gives c count (m, n I ) = 873 log 2 (m + n I ) − 968 .
The true answer lies in between (G2) and (G3); the effect is small and henceforth we ignore the difference arising from controlling on |0 versus |1 . We therefore write the final answer as It is possible to rearrange the particle representation given in (A1) to use only 2 controls for all, but subsequent operations become more complicated in this case.

The second sub-operation, U (m) e
Let's now look at the operation in which we determine whether or not we had an emission, whose circuit is shown in Figure 4. If we are at the m th step, the largest number of particles we can have is m + n I , while the minimum is n I . This means that we have to apply U e gates controlled on all the possible combinations of three integers, ranging from 0 to m + n I , whose sum is in the range [n I , m + n I ]. There are c(m, n I ) = m + 1 6 (m 2 + 3 m n I + 5m + 3n 2 I + 9n I + 6) (G5) such such combinations. For each of these we run a C (3 log 2 (m+n I ) (U e ) operation, where the U e gates are R Y (θ) rotations. Using the results from above about C (n) (U ) operations, the total number of standard gates necessary for the emission operation is N sub2 (m, n I ) = c(m, n I ) (96 log 2 (m + n I ) − 27) .
of the |0 state. Therefore, the number of standard gates necessary to apply the controlled-U − operations when simulating the m th step is c count (m, n I ) given in (G3).
In the first part of the sub-operation circuit, the gate U h has the same controls on the count registers that we found in the emission operation, plus having controls on the particle state and on the emission qubit (which has to be in the |1 state for an emission to have happened). As we have seen the particle state can be in any of six states specified by three qubits; however, the emission probability is the same for particle and anti-particle, and no U h is the identity if the particle is in state |0 . Therefore the number of possible combinations of particle state and count states in sub-operation j is 3 c(m, n I ), which is the number of times we must apply U h gates controlled on 4 + 3 log 2 (m + n I − j) qubits. The history register contains states labeled by integers from 0 to m + n I and we use the standard bit representation to encode these integers in qubits. For the j th sub-operation the matrix U h is an a × a unitary matrix, where a = 2 b with b = log 2 (j) . Lastly, we must break down the operation in which we adjust the particle states in |p . Given the transformation in Eq. (F3), we can implement U p efficiently as shown in Figure 10. In the circuit H is the Hadamard gate and U r is given by Here, k = m + nI .
The operation U p is controlled on the possible states in |h . There are m + n I such states, each requiring log 2 (m + n I ) controls. Thus, for each of the m + n I occurrences of U p one adds log 2 (m + n I ) controls to each operation in Fig. 10. This gives N sub4 (m, n I ) = (m + n I ) (224 log 2 (m + n I ) + 143) .

Summary
Adding all sub-operations together and summing over 0 < m < N − 1, one finds that the overall scaling of our circuit is N 5 ln N . Fig. 11 shows the number of gates as a function of N for N < 50.
Note that one can obtain a much shallower circuit requiring less qubits if one takes into account that in the end states with different history registers do not interfere with one another. This implies that one can measure the history register after the third operation in each step, and reset it back to zero. This collapses the quantum state to one with a definite history. Having a state with a definite history gives definite knowledge about the number of bosons n φ , as well as the total number of particles n tot . This is because the history allows us to infer how many emissions have happened, which means that the state has a definite number of particles; since one also knows at any step at which an emission happened if the emitting particle was a fermion or a boson, one knows the total number of bosons. Thus, instead of counting and keeping track of the 3 values n φ , n a and n b , it suffices to only keep track of n a and from that derive n b = n tot −n φ −n a . Following similar steps outlined in this section in this case, one can easily see that the scaling of the depth of the circuit is reduced significantly, with an overall scaling of N 3 ln N , instead of the N 5 ln N . Unfortunately, current quantum hardware does not allow for such repeated measurements. Appendix H: Results for two steps Figure 12 presents an overview of the possible outcomes from running the full circuit with two steps. In this case with N = 2, 24 total qubits were used. Several interesting features can be observed from this result. First, given that the initial state in a single f 1 fermion, for g 12 = 0 one can only end up in a state containing an odd number of f 1 fermions. In contrast, for g 12 = 1, even with no emissions the fermion can change type. This can be traced to virtual corrections involving fermion changing interactions. This figure also gives a sense how big the effect of φ → ff splittings is (first three sets of bars versus last three sets of bars) for the chosen parameter values. This simulation is also used to predict the distribution of observables shown in Fig. 2.
Appendix I: Circuit with no φ → ff Ignoring the φ → ff splittings, ignoring the running coupling, and starting with only one fermion (possibly in a superposition) as the initial state, allows us to drastically simplify our quantum circuit, since all one needs now is a single qubit which represents the fermion flavors, and a boson register, which keeps track of whether or not a boson was emitted at a given step. This boson register is the equivalent to the emission register plus the particle register in the general circuit. We no longer need a history register, since we know the fermion is the only particle which can emit, nor do we need the count reg-isters since in this limit the probability of a boson being emitted only depends on the flavor of the fermion. The full evolution can be carried out with the much simpler circuit shown in Figure 13.
The U and U † gates are the same as in Eq. (B5), while the U a/b i gates are given by the matrices which encode the amplitude for the fermion to emit or not emit a boson at a given step. These gates are controlled on the fermion state since the gate parameters depend on the flavor of the fermion. The circuit construction demonstrates that the scaling for generating a single event is linear with the number of steps. At each step the U a and U b gates are conditionally applied to a new qubit, but after that the qubit is left alone until the final measurement at the end of the evolution. Therefore, at each step one could measure the qubit on which the U a/b gates act on, store the result in a classical register, reset it to the initial |0 state and reuse it for the next step. Using this method of repeated measurements and resetting the measured qubits one can rewrite the circuit in terms of just two qubits as shown in Figure 14. At each step one records the measurement on the second qubit, and at the very end the first qubit is measured. The combination of these measurements makes up one event. Note that because this circuit can be implemented using just 2 qubits, one can in fact find an efficient quantum inspired classical algorithm. This is derived in Appendix K.
We now discuss how to implement this circuit on currently available hardware. Given that repeated measurements as used in Figure 14 is not possible on currently existing hardware, we use the circuit shown in Figure 13. To implement this, one needs to break down the controlled operations into standard gates, namely single qubit gates and CNOT gates. To achieve this, one first uses the well known result In our case the gate U consists of a R Y (θ) rotation gate. Furthermore, we use the fact that for an arbitrary controlled-U operation, one has To apply this to the controlled-R Y (θ) gate one chooses where α, β and ψ satisfy This gives gates A, B, C, P (where P is the trivial identity matrix) that satisfy all conditions. Using this information one finds that each step requires a total of 12 simple quantum gates (8 single qubit gates and four CNOT gates), and in addition two transformations are required at the beginning and end of the circuit which also consist of single qubit gates. Generating a single event therefore requires a total of n gates = 12 N + 2 (I6) single qubit and CNOT gates.
Appendix J: Renormalization of the theory When computing higher order corrections to the simple model given in Eq. (1) one encounters divergences as in any quantum field theory. These divergences can be removed using a renormalization procedure by including counter-terms which are then determined using a socalled renormalization prescription. This is accomplished by writing the bare fields, masses and couplings of the unrenormalized theory [denoted with a superscript (0)] in terms of renormalized quantities, with the relation given by renormalization factors. This gives Starting from the bare Lagrangian given in Eq. (1), one therefore finds the renormalized Lagrangian where repeated indices are summed over, and we have defined g ii ≡ g i and g ij = g ji . The divergences of the counter-terms given in the last 4 lines of Eq. (J2) need to cancel all divergences arising from higher-loop diagrams, but their finite terms need to be fixed by a renormalization scheme. The simplest renormalization scheme is the so-called minimal subtraction (MS) scheme, in which the finite contributions of the counter-terms are defined to vanish, however other schemes are possible (and widely used). A reader not too accustomed with quantum field theories might be worried that one could choose a scheme that eliminates all contributions giving rise to the mixing between the fermions considered in this work. However, one needs to keep in mind that the loop corrections in the theory give rise to a scale dependence, and that the interactions in the renormalized theory depend on a scale called the renormalization scale. The renormalization scheme provides a relation at a single scale only, which implies that it is not possible to eliminate mixing effects at a general renormalization scale.
However, the fact that all renormalized quantities depend on a renormalization scale implies that the basis change that eliminates the fermion mixing is scale dependent, since it depends on the values of the scale dependent coupling constants. The evolution of a parton shower can be viewed as an evolution in the renormalization scale, such that each step in our algorithm corresponds to a different value of the renormalization scale. This explains why the rotation matrix U in Eq. (B5) depends on the step.