{'fraction_non_alphanumeric': 0.035977056962025314, 'fraction_numerical': 0.011867088607594937, 'mean_word_length': 4.578221535745808, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Kinetic Theory and Stochastic Simulation of Field Quanta

We develop quantum electrodynamics into a kinetic-theory-like evolution equation for electrons, positrons and photons. To keep the"collision rules"simple, we make use of longitudinal and temporal photons in addition to the usual transverse photons. For our explicitly time-dependent approach, we introduce proper time-correlation functions. We then develop a stochastic simulation technique for solving the resulting kinetic equation and we comment on the subtleties associated with gauge invariance, Lorentz invariance and divergent integrals. To illustrate the validity and power of the proposed ideas, we show how a very simple simulation of the dynamics of field quanta can be used to obtain the main contribution to the anomalous magnetic moment of the electron.


I. INTRODUCTION
The goal of this work is to introduce particle-like simulations for quantum field theories. Rather than with particles in a classical sense, of course, we actually deal with field quanta. In the course of time, these field quanta undergo collisions of the types that we commonly represent by Feynman diagrams. We develop and illustrate the ideas in the context of quantum electrodynamics (QED), where the collisions are given by (i) emission or absorption of photons by electrons or positrons and (ii) electronpositron pair production or annihilation. The simulations to be developed consist of a stochastic sequence of such collisions.
Our particle-based simulations are supposed to be developed into an alternative to the widely used fieldbased simulations, which go back to K. G. Wilson's famous formulation of lattice gauge theory [1]. Computer simulations of lattice gauge theories with dynamic fermions [2,3] have been established as a very successful tool in nonperturbative quantum field theory, but they are extremely demanding from a computational point of view (four-dimensional lattices, continuous evolution by hybrid-molecular-dynamics algorithms, conjugate-gradient calculations). It is hence interesting to explore complementary simulation ideas.
In Section II, we introduce the ideas and notation of the four-photon approach to free electromagnetic fields and Dirac's spinor representation of free electrons and positrons. Interactions are then introduced in Section III, where also the relativistic covariance and the gauge invariance of the theory are discussed. Proper correlation functions for our explicitly time-dependent approach are presented in Section IV; in particular, we consider propagators, vertex functions and the electron form factor leading to the magnetic moment of the electron. We then have all the tools to develop a kinetic-theory-like simulation method for the field quanta in Section V and we present a particularly simple simulation that allows us to obtain the main contribution to the anomalous magnetic * hco@mat.ethz.ch; http://www.polyphys.mat.ethz.ch/ moment of the electron. A brief summary and a detailed discussion conclude the paper (Section VI). Some useful definitions and relations are compiled in the Appendix.
Throughout this paper, we use natural units with = c = ǫ 0 = 1 where is the reduced Planck constant, c is the speed of light, and ǫ 0 is the electric permittivity. In these units, the electric charge is given by e 0 = √ 4πα ≈ 0.30282212, where α is the fine-structure constant. Only one further unit, naturally taken as mass or energy, remains to be specified.

II. FREE QUANTUM FIELDS
We begin our development by introducing the free electromagnetic and electron-positron spinor fields. On the one hand, we lay out the basic notation and constructions, on the other hand, we offer a short but selfcontained introduction to QED in a non-standard approach with longitudinal and temporal photons.

A. Vector potential
To avoid the need of considering Coulomb interactions explicitly, we use the four-photon quantization of the electromagnetic field. This elegant idea was originally developed in 1950 by Gupta [4] for free electromagnetic fields and by Bleuler [5] in the presence of charged matter.
The four-vector potential with components A µ is introduced as a Fourier transform in terms of polarization states and the corresponding creation and annihilation operators, (1) We denote the corresponding Fourier component in the integrand by The temporal unit vector, is actually independent of q. The three orthonormal spatial polarization vectors are chosen as and The polarization vectors n 1 q and n 2 q correspond to transverse photons, n 3 q corresponds to longitudinal photons. Note the symmetry property n α −q = (−1) α n α q .
The sign ε α (defined as +1 for transverse and longitudinal photons, −1 for temporal photons) leads to a nonself-adjoint nature of the vector potential (1). The spatial components are self-adjoint, whereas the temporal component A 0 is anti-self-adjoint. The latter statement can be expressed as where we have assumed a Minkowski metric η µν = η µν with signature (−, +, +, +), that is, η 00 = η 00 = −1. It is hence important to note that Eqs. (1) and (3) specify A 0 with a lower index. With these sign conventions, we further obtain a modified completeness relation for the polarization vectors, α ε α (n α q ) µ (n α q ) ν = η µν .
The four photon creation and annihilation operators, a α † q and a α q , introduced in Eq. (1) satisfy the usual commutation relations, We can hence build the four-photon Fock space in the usual way by multiple application of all the operators a α † q for all q ∈ R d on a ground state (which is annihilated by any a α q ). The full Hilbert space factorizes into spaces obtained by repeated application of a α † q for each mode q (see, for example, Secs. 1 and 2 of [6] or Secs. 12.1 and 12.2 of [7] for more details on the construction of such Fock spaces).
There is only one modification of the standard construction: a modification of the scalar product, which we present as the construction of bra-from ket-vectors. If a basis vector of the ket space contains n temporal photons, the corresponding standard bra-vector is multiplied by a factor (−1) n . This minus sign associated with each temporal photon in the bra-space is introduced such that the vector potential (1) becomes self-adjoint for the modified scalar product. The use of an indefinite scalar product may be alarming because it might endanger the probabilistic interpretation of the results as, for example, the norm of an odd-number temporal photon state is negative. However, no interpretation problems arise for physically admissible states, where admissibility is defined in terms of a proper version of the covariant Lorentz gauge condition. In particular, the Lorentz condition implies that physical states involve equal numbers of longitudinal and temporal photons having the same momentum (see Chapter 17 of [8]), which gives us an idea of how interpretational problems are avoided. The physical states have been listed in an elementary way in (2.16) of the pioneering work [5]. An elegant way of characterizing the physical states has been given in Section V.C.3 of [9] by transforming from longitudinal and temporal to left (gauge) and right (droite) photons, and Physical states can contain an arbitrary number of gphotons, but no d-photons. Finally, the free Hamiltonian for our massless photons is given by where all four photons are clearly treated on an equal footing. Their energy is given by the relativistic expression for massless particles. The free Hamiltonian (13) is self-adjoint both for the standard and for the indefinite scalar product. The observable magnetic field, as obtained from the rotation of the spatial part of the vector potential (1), is given by the spatial components of Note that the magnetic field involves only transverse photon operators. The observable Fourier components are associated with a 1 † q −a 1 −q and a 2 † q +a 2 −q which, in view of the symmetries (7), coincide with the transverse Fourier components of the vector potential in Eq. (2).

B. Dirac fields and spinor properties
According to Eq. (13.59) of [7], the Hamiltonian associated with the Dirac equation for the free electron/positron can be written as where E p = p 2 + m 2 gives the energy of a relativistic particle with mass m and momentum p, and a summation over the possible spin values σ = ±1/2 is implied by the same index occurring twice. The operators b σ † p and b σ p create and annihilate an electron of momentum p and spin σ. Similarly, the operators d σ † p and d σ p create and annihilate a positron of momentum p and spin σ. As we are now interested in fermions instead of the bosonic photons, we need to specify the fundamental anticommutation relations for the creation and annihilation operators. All particle and antiparticle creation operators anticommute among each other, and so do the annihilation operators. The only nontrivial anticommutation relations are and The four-component spinor field associated with annihilating an electron or creating a positron, each coming with the two values ±1/2 of the spin, has the Fourier representation (see Eq. (13.50) of [7]) with the adjoint field operator creating an electron or annihilating a positron, where an asterisk implies both complex conjugation and transposition of a vector or matrix. Equation (19) should be compared to Eq. (1); instead of the polarization vectors for photons, the spinors u and v occur for relativistic electrons and positrons. In the standard representation, they are given by and v −1/2 where we have introduced the notationp j = p j /(E p + m) for conveniently normalized momenta p j = p j .
In the context of a relativistic theory, it is useful to introduceψ = ψ † γ 0 , where γ 0 is a diagonal 4 × 4 matrix with the diagonal elements (1, 1, −1, −1), becauseψ has a more natural Lorentz transformation behavior than ψ † (see Section 2.2 of [10]). In particular,ψ(x)ψ(x) is a Lorentz scalar, whereas ψ † (x)ψ(x) is not. With the corresponding definitionsū = u * γ 0 andv = v * γ 0 , Eq. (19) can be multiplied by γ 0 and rewritten as (24) We introduce the Fourier components in Eqs. (18) and (24) as For p = 0, the column vectors (20)-(23) clearly form an orthogonal basis of the four-dimensional vector space characterizing the spins of electrons and positrons. For arbitrary p, we still have orthogonality and completeness relations. Actually, there are two versions, one for spinors with equal p and one for spinors with opposite p. For example, one form of the orthogonality relations reads To establish the corresponding completeness relation, we consider the 4 × 4 matrices obtained as sums over tensor products of spinors and write them in a compact form in terms of the Pauli matrices (see appendix) and 2 × 2 unit matrices. We find and where the Dirac 4 × 4 matrices γ j required for the compact reformulation in the last step are listed in the appendix. For p = 0, we realize that Λ e (p) and Λ p (p) may be regarded as projectors to the electron and positron degrees of freedom, respectively. In Eqs. (27)- (30), there is no integration over p and we indicate the summations over σ explicitly. Equations (29) and (30) imply the completeness relation and the additional property

III. INTERACTIONS
The interaction between charged leptons and photons is given by the Hamiltonian where we have made use of the electric charge density, and the electric current density, Equations (34) and (35) define a four-vector. Note that, in the standard scalar product, the Hamiltonian H (1) is not self-adjoint because A 0 (x) is anti-self-adjoint. However, in the indefinite scalar product, the Hamiltonian H (1) becomes self-adjoint. As suggested in Eq. (13.61) of [7], we actually use the normal-ordered versions of the charge and current densities (all creation operators are to the left of all annihilation operators in a product) to avoid irrelevant infinite contributions.

A. Collision amplitudes
The interaction term (33) involves two spinor fields and one vector potential, and hence such a ternary interaction corresponds to one photon and two lepton lines in a Feynman diagram. After inserting the Fourier transforms (1), (18) and (24), we write the interaction between Dirac and electromagnetic fields in terms of two contributions, pp . The first contribution describes the spin dependent emission and absorption of photons by electrons and positrons, the second term accounts for pair production and annihilation. These contributions are of the following form, and The complex collision amplitudes h σσ ′ µ pp ′ andh σσ ′ µ pp ′ are proportional to the electric charge e 0 and fully characterize electromagnetic interactions of photons with electrons and positrons. For the functional form of the collision amplitude h σσ ′ µ pp ′ characterizing photon emission, we find The collision amplitudeh σσ ′ µ pp ′ for pair production is given byh These collision amplitudes have the following symmetries, and which are useful for verifying the adjointness properties of H Let us consider the following (non-normalized) states consisting of n electrons,n positrons andñ photons, where n l=1 p l + n l=1p l + ñ l=1 q l is the total momentum of the multi-particle state. These states are the base vectors of the Fock space that includes electrons and positrons in addition to photons and hence properly generalizes the Fock space of Section II A. The states (42) are eigenstates of the energy of the noninteracting system with eigenvalues pp on a state (42). Although the calculation is straightforward and the outcome is quite lengthy, we give the explicit result because it clarifies why we refer to a "kinetic-theory-like" approach and to "collision amplitudes:" σ 1 · · · σ n σσ 1 · · ·σnσ α 1 · · · α l−1 α l+1 · · · αñ p 1 · · · p n pp 1 · · ·pn q l − p q 1 · · · q l−1 q l+1 · · · qñ − n j=1n k=1ñ l=1 Equation (43) summarizes the eight collision rules of QED for electrons, positrons and photons. The first three terms decrease the number of particles by one, either through absorption of a photon by an electron or positron or through annihilation of an electron-positron pair into a photon. These terms contain double sums corresponding to the selection of two particles (to be replaced by one) and do not involve integrations. The next three terms increase the number of particles by one, either through emission of a photon by an electron or positron or through creation of an electron-positron pair from a photon. The hallmark of these terms is a single sum corresponding to the selection of a disappearing particle (to be replaced by two) and a momentum integration for one of the two new particles. The last two terms decrease or increase the number of particles by three, an electron, a positron and a photon in either case. They contain a triple sum, corresponding to selecting three particles, and a double integral, corresponding to the freedom of choosing two momenta, respectively. Momentum conservation is respected by all of the eight terms, and each of the terms is proportional to a collision amplitude h orh. These "collision rules" are the core of the simulation of QED to be developed below.

B. Lorentz and gauge invariance
After defining the Hilbert space, the fields and the various contributions to the Hamiltonian, we should now ask whether the resulting theory has the physically required invariance properties. The proper relativistic form is not at all obvious because the Hamiltonian approach assigns a special role to time. In the Hamiltonians (13) and (15) for the free electromagnetic and electron/positron fields, the relativistic energy-momentum relations for massless and massive particles provide the proper ingredient. For the interaction Hamiltonian (33) to be Lorentz invariant, A µ and J µ have to be Lorentz four-vectors. For A µ , this is nicely achieved by the four-photon approach which puts temporal and spatial components on an equal footing, except for the sign of the Minkowski metric properly arising from the use of an indefinite scalar product. For J µ , Lorentz transformation behavior is achieved by construction of the spinors, which are chosen to be trivial for particles at rest and then need to have a suitable momentum dependence for moving particles. The correct relativistic transformation behavior is hence hidden in the collision amplitudes (38) and (39).
By introducing four types of photons, we are dealing with more than the physically observable degrees of freedom. In other words, our Hamiltonian dynamics takes place in an unphysically large space. The benefits are the simplicity and elegance of the equations in the larger space. The drawback is the need to identify the physical observables. As we rely on the approach of Gupta and Bleuler, this problem has been solved by introducing left and right photons, where the latter correspond to unphysical degrees of freedom. We here play safe and consider only correlation functions involving the transverse photons directly related to the observable magnetic field.
One might be worried about the hidden character of fundamental physical invariance properties. However, experience with the standard form of the Maxwell equations for electromagnetic fields and the usefulness of conveniently fixed potentials (fixed gauge) clearly underlines that such a procedure is very common and absolutely unproblematic.

IV. CORRELATION FUNCTIONS AND MAGNETIC MOMENT
As we do not treat time and space on an equal footing and as we work in a particular gauge, it is important to introduce the appropriate correlation functions carefully. We here introduce and discuss some useful time-ordered two-and three-time correlation functions.

A. Two-and three-time correlations
We first consider the two-time correlation function of two operators A and B. Actually, it is more convenient to study Fourier transforms. When time-ordering is properly taken into account, the correlation in the frequency domain is given by where we have introduced the evolution operator The plus or minus sign in the definition (44) depends on how A and B are related to commuting or anticommuting operators, the small parameter ǫ ensures the convergence of the time integral. Note that R ω involves the full Hamiltonian and possesses the perturbation expansion, where R (0) ω is defined as R ω in Eq. (45), but with the free Hamiltonian H (0) instead of H. The simplicity of this perturbation expansion for R ω actually provides the motivation for considering Fourier transforms. A full perturbation theory of the correlation functions (44) requires an additional expansion of the vacuum states. It is given by the following relations between the vacuum states of the free and interacting theories, To obtain this result, it is important to realize that H (1) is self-adjoint for the modified, indefinite scalar product.
In the limit of a vanishing regularization parameter, the operators R (0) ω are self-adjoint for both the modified and the standard scalar products. Note the similarity of the perturbation expansions of the vacuum state and of R ω for ω = 0.
For time-ordered three-time correlation functions, the proper Fourier transform can be expressed as where the signs have been chosen for the case that A and C are fermion operators, whereas B is a boson operator. Note that the frequency ω ′ occurs all the way between A and B, and ω between B and C. If the order of A and B (or B and C) is reversed, the frequency ω ′ (or ω) picks up a minus sign.

B. Propagators
The standard form of Feynman's fermion propagator S αβ in Fourier space can be introduced by choosing A and B in Eq. (44) as Fourier components (25) and (26) of the fermion field. More precisely, we define where, for fermion operators, the minus sign must be used in the definition (44). In space-time, this corresponds to the usual definition in terms of a time-ordered product (see, for example, Eq. (13.72) of [7] or p. 112 of [11]), In a more compact notation, Eq. (50) can be expressed as where x = (t, x) and 0 = (0, 0) are four-vectors and T is the time-ordering operation.
For the free theory, it is straightforward to evaluate the correlation function in Eq. (49) for the Fourier components given in Eqs. (25) and (26). We obtain the free propagator In the second step, we have used Eqs. (29) an (30) to obtain the standard form of the free fermion propagator. A photon propagator can be introduced in the same way as the fermion propagator (49), but now with the plus sign for boson operators in the definition (44), From a straightforward calculation, we obtain the free photon propagator This result illustrates that the four-photon approach corresponds to a particularly convenient gauge.

C. Vertex functions
As important further correlation functions, we introduce the vertex functions Γ ν through special cases of a three-time correlation function, where the Γ ν are 4 × 4 matrices in spinor space, just like the fermion propagators S. In the compact four-vector notation, we can write where we have made use of the following four-vectors: With the propagators S(ω, p), D(ω, q) and the vertex functions Γ ν (ω ′ , p ′ , ω, p) defined in Eqs. (49), (53) and (55), we have now introduced the most important correlation functions QED.

D. Perturbation theory
To gain a better understanding of the vertex functions Γ ν (ω ′ , p ′ , ω, p), we would like to calculate the lowestorder contributions by perturbation theory. As the definition (55) involves three fields, the vertex functions vanish for the free theory. Only odd orders of perturbation theory can contribute. Even the first-order result, is already tedious to get within our approach, which is highly inappropriate for constructing perturbation expansions. It is hence recommended that perturbation theory should be based on compact expressions of the form given in Eqs. (51) and (56), for which elegant and powerful methods based on Feynman diagrams are available in any textbook on quantum field theory. Following the standard rules (see, for example, Appendix B of [7] or Appendix C of [11] with the explicit result given in Eq. (III.6.8) on p. 197) we obtain where this contribution results from the Feynman diagram shown in Fig. 1. By inserting the first part of Eq. (52) as well as Eq. (54) and using Cauchy's integral formula to perform the frequency integration over q ′ 0 , we can rewrite the expression (58) as a three-dimensional integral, One still recognizes the importance of fermion propagators, but now split into electron and positron contributions.

E. Form factors and magnetic moment
Perturbation theory suggests that Γ ν can be regarded as an effective interaction vertex when photons and electron-positron pairs are around during the basic ternary interaction process described by γ ν . Whereas γ ν characterizes the electromagnetic interaction of a structureless point electron, the full interaction including a cloud of photons and electron-positron pairs around the electron during interaction with an external photon is characterized by Γ ν . The physics of the effective interactions can be described elegantly in terms of the electromagnetic form factors F and G, which occur in the most general representation of four-vectors (see, for example, Section 10.6 of [12] or Section III.6 of [11]), where the real functions F (x) and G(x) depend on They satisfy the normalization condition F (0) + G(0) = 1, which is consistent with the result of first-order perturbation theory, that is, with the constant form factors F (x) = 1 and G(x) = 0. In view of the definition (38), it is natural to define the effective collision amplitude The magnetic moment, which characterizes the interaction of a charged particle with an electromagnetic field, is given by the g-factor and can be expressed in terms of the form factor (see, for example, Section 10.6 of [12] or Section III.6 of [11]) In lowest-order perturbation theory, we find g = 2, whereas the experimental result is slightly larger, where significantly more accurate experimental results would be available (the relative experimental error for g is smaller than 10 −12 ). The celebrated approximate one-loop result of Schwinger [13] corresponding to the Feynman diagram in Fig. 1 is Our goal is to find the anomalous g-factor deviating from 2 by kinetic-theory-like simulations.

V. SIMULATIONS
After introducing the proper formulation of Hamiltonians and correlation functions in the four-photon approach, we are now in a position to develop a kinetictheory-like simulation methodology. The eight collision rules in Eq. (43), together with the collision amplitudes (38) and (39), provide the corner stones for the stochastic simulation technique to be developed here.

A. Basic ideas
The simulations proposed in the present paper are strongly inspired by the powerful stochastic simulation techniques for quantum master equations [14][15][16][17]. Instead of looking for a deterministic solution of the Schrödinger equation with H (0) = H pp , we consider a stochastic process |ψ t in Hilbert space such that the expectation E |ψ t satisfies Eq. (66). Let us consider an important example of a stochastic process |ψ t . In that example, the process starts from a deterministic initial condition |ψ 0 and continuously evolves according to the Schrödinger equation (66) for the free Hamiltonian in the absence of interactions, Instead of having a continuously acting Hamiltonian H (1) , we introduce a rate parameter r that determines the occurrence of random times at which the following jump occurs, The expectation of the stochastic process with continuous evolution (67) and random jumps (68) occurring with rate r solves the Schrödinger equation (66). The inverse rate 1/r plays the role of a time step that, however, occurs randomly according to an exponentially decaying distribution with average 1/r. For large r, we recover an almost continuous evolution with small jumps. For small r, we have rare but large jumps according to the prefactor 1/r in Eq. (68). The weaker the interaction, the smaller can the jump rate r be chosen, and the more efficient becomes the simulation. The proposed stochastic technique is expected to be particularly powerful when a problem is amenable to perturbation theory. By choosing a finite rate r, nonperturbative results can be obtained. According to Eq. (43), the action of the Hamiltonian H (1) in the jump (68) involves sums and integrals. In our stochastic approach, it is most natural to perform such sums and integrals by Monte Carlo methods. This introduces a second source of randomness into our Hilbert space process. After averaging, we still obtain a solution of the Schrödinger equation (66).
In n jumps of the type (68), the trajectory typically progresses by n times 1/r. In view of the exponential time distribution, however, all positive values of time increments are possible, with the following probability density for the n-jump contribution at time t, with the waiting times τ 1 . . . τ n for the jumps 1 . . . n to occur. The last factor is the probability for having no further jump between t ′ and t. In view of the δ function in Eq. (69), it is particular convenient to calculate the n-jump contribution to the Fourier transform (45) of the evolution operator. We then obtain For r = 0, this corresponds to the usual perturbation theory in Eq. (46), but the summed result is actually independent of r. So, we can eliminate the stochastic aspect of the jump rate in calculating Fourier transformed time correlation functions, but the randomness associated with Monte Carlo summations and integrations remains in the simulation. By varying the rate parameter r in our simulations, we can go all the way from continuous integration (large r) to perturbation theory (small r). The simulation corresponds to the perturbation theory given by Eqs. (46) and (47). However, the formal regularization parameter ǫ is replaced by the rate parameter r, where 1/r is an average time step and ir is added to the interaction. We have thus found a very simple recipe for how to proceed from perturbation theory to an integration scheme with average time step 1/r.

B. Procedure
We are now in a position to describe the stochastic simulation procedure for calculating correlation functions. During the simulation of a trajectory in Hilbert space, we always have exactly one of the base vectors (42) together with a corresponding complex coefficient. The free evolution between interactions is easy to evaluate because, during the entire simulation, we always keep eigenstates of the energy of the noninteracting system, as given in Eq. (42). Free evolution merely changes the coefficient but not the base vector.
Let us assume that we wish to calculate the time correlation function 0| Ae −iHt B |0 or its Fourier transform 0| AR ω B |0 occurring in Eqs. (44) or (50). The basic idea is to start from the vacuum state |0 (0) of the free theory, to initiate by a sufficiently large number of time steps, to apply the operator B, to perform the time steps for the evolution from B to A, to apply the operator A, to finalize by a sufficiently large number of time steps, and to project onto 0 (0) |. "Initializing" and "finalizing" should be recognized as the counterparts of the two parts of Eq. (47); in a sense, we need to find the vacuum state |0 of the interacting theory by "equilibration" when we start from the vacuum state of the free theory. For these equilibration steps, we should use Eq. (70) with ω = r = 0. Going to the limit r = 0, or a very large time step, is certainly appropriate for producing the vacuum state of the interacting theory. When we apply R ω for the evolution from B to A, the rate parameter r and the number of interactions n in Eq. (70) should be adjusted to the frequencies ω we would like to explore (nω ∼ r).
Whenever we apply H (1) according to Eq. (43), there is a number of summations and integrations to be carried out. For each sum or integral we choose only one single term; the full sums and integrals are obtained by averaging over many realizations. For example, we can choose the probability p rule (j) for using one of the collision rules j = 1, . . . 8 to be 1/8 and then multiply the estimate for the randomly chosen term j by 8 (more generally, by 1/p rule (j)). The other summations and integrations in Eq. (43) are performed with randomly selected particles, spins, polarizations or randomly selected momenta. Whenever we make a random choice, we need to introduce the reciprocal of the corresponding probability as a weight factor. In the final part of the simulation, one should only admit random selections that can still lead towards the free equilibrium state, with correspondingly larger probabilities and smaller weight factors.
A subtlety is introduced by the δ function in rule number seven in Eq. (43), which reduces the number of integrations by one. As the rule is associated with the annihilation of a positron and we typically choose initial and final states that do not contain any positrons, we can always use the δ function to fix the momentum of the positron when it is created. As long as a positron exists in the simulation, we need to leave it open whether the positron can be randomly selected (if it is eventually annihilated by rule 3) or it is fixed (if it is eventually annihilated by rule 7).
Note that, in principle, the proper averages are obtained for any choice of the probabilities for randomly selecting rules, particles, spins, polarizations or momenta. However, the efficiency of the simulation, or the size of the statistical error bars, depends crucially on the choice of these probabilities for the simulation algorithm. The more uniform the weight factors, the more efficient is the simulation. The art of developing good Monte Carlo simulations consists in finding probabilities that reflect the summations and integrations in Eq. (43) as faithfully as possible.

C. Truncations
In a practical implementation of the simulation procedure, a number of truncations are required. Most obviously, when the expression (70) for R ω is used, only a finite number of interactions can be kept (the number of interactions corresponds to the number of time steps). The smaller the frequency ω, or the larger the probed time scale, the more interactions need to be kept. The required computer time grows linearly with the number of interactions.
It is convenient to limit the maximum number of electrons, positrons and photons present in a simulation at any given time. Limiting the number of simultaneously present field quanta corresponds to limiting the order of perturbation theory. Note however that, in the course of the simulation, a much larger number of electrons, positrons and photons can occur (and disappear). In that sense, we deal with a renormalized perturbation theory with interactions on all scales.
In practice, we actually set minimum and maximum values for the photon momenta to avoid possible infrared and ultraviolet problems. This is done through the choice of momenta according to a probability density f (q). If one starts with a probability density for q, one needs to find an algorithm to sample according to that distribution. We hence prefer to start with an algorithm for generating q and then determine the corresponding probability density. As we anticipate logarithmic divergences at large and small wave vectors in QED, a natural choice for the length of the wave vector of a photon is given by where q min and q max are the lower and upper cutoffs and r is a uniform random number on the interval [0, 1]. If the orientation of the vector q is chosen randomly, this choice corresponds to an isotropic power law probability density for the three-dimensional wave vector given by As an alternative, one can try a uniform distribution of q, which implies the probability density The most efficient choice of the probability density f (q) is expected to depend on the quantity of interest. Fig. 2.

D. Example
In order to illustrate our simulation approach, we consider a very simple simulation in which the interaction is applied only three times. Even such a toy simulation contains interesting physics. According to Fig. 1, such a simulation should be sufficient to get the leading-order correction to the effective collision amplitude and hence to the magnetic moment of the electron. The Feynman diagram in Fig. 1 represents the contribution to the vertex function given in Eq. (58), that is, it involves momentum four-vectors and four-dimensional integrations. Our Hamiltonian time-evolution approach, on the other hand, is focused on time evolution, momentum three-vectors and three-dimensional integrations, as is clearly visible in all expressions for the various Hamiltonians. We need to identify the possible sequences of three interactions that are consistent with the Feynman diagram in Fig. 1 in our approach. If we amputate the external legs from this Feynman diagram and focus on the core triangle representing Γ µ(3) , we obtain the Feynman diagrams for the corresponding three-step simulation shown in Fig. 2. For our kinetic-theory-type simulation, these six diagrams are not equivalent because they correspond to the six possible permutations of three interaction events occurring in the course of time; they correspond to the six terms of Eq. (48), in the same order, from which we get the respective frequencies and signs coming with each term (possible minus signs are actually taken into account by the amputation procedure so that we can use positive signs for all terms contributing to the vertex function). The information about the full diagrams with external legs can be reconstructed by extending the fermion line by an incoming electron of momentum p and spin σ and an outgoing electron of momentum p ′ and spin σ ′ , and by attaching an incoming photon line with momentum q = p ′ − p and spacetime index µ to the kink of the fermion line. Therefore, the six diagrams in Fig. 2 specify the three-step simulations for the vertex function in a unique way.
The first two Feynman diagrams in Fig. 2 correspond to absorption of an external photon by an electron or positron, the next two diagrams correspond to pair production, and the last two diagrams correspond to pair annihilation. The simulation steps associated with the first diagram are shown in detail in Fig. 3. In the middle, we show the sequence of states of the form (42) for the trajectory in Fock space; for this diagram, the visited states consist of only one electron and of up to two photons. On the left-hand side, we show the collision amplitudes occurring for each of the three interactions according to the collision rules in Eq. (43). For the con-tribution associated with the kink in the fermion line, the polarization vector is missing because, in the definition (62) of the effective collision amplitude, we have multiplied Γ µ with the spinors but not with the polarization vectors. On the right-hand side of Fig. 3, we show the factors resulting from the free evolution according to R (0) with the frequencies obtained from Eq. (48) where, in view of Eq. (60), we choose ω = E p and ω ′ = E p ′ (Eq. (55) then implies that the incoming and outgoing electrons are on the mass shell). To amputate the external legs, only the propagating particles shown in Fig. 2 need to be taken into account for the energy calculations.
For example, in the first step of Fig. 3, an electron with momentum p and spin σ emits a photon with momentum q ′ and polarization α ′ , whereupon the electron momentum changes to p − q ′ and the spin toσ. From the fourth rule in Eq. (43) we hence read of the collision amplitude hσ σν p−q ′ p (n α ′ q ′ ) ν , which is the first factor on the left-hand side of Fig. 3. Of the electron and the two photons of the resulting state, only the electron (with energy E |p−q ′ | ) and one photon (with energy q ′ ) are propagated in the lower half of the first Feynman diagram in Fig. 2 because the external photon (with energy q) is amputated. For the action of the free evolution operator, this implies where the frequency ω has been identified from the first term of Eq. (48), which corresponds to the first Feynman diagram in Fig. 2. With ω = E p , we obtain the first factor on the right-hand side of Fig. 3. The remaining factors can be found by analogous arguments. A sequence of steps like the one shown in Fig. 3 can be associated with each of the six Feynman diagrams in Fig. 2 (a further example is given in Fig. 4). By collecting all the factors associated with each of the six diagrams and using Eq. (9), we obtain the following effective collision amplitude for a three-step simulation (the order of the terms corresponds to the order of the Feynman diagrams in Fig. 2 and hence also to the order of the terms in Eq. (48)): With the exception of photon emission by a positron, all possible collision rules have been used to obtain Eq. (76), so that our toy simulation is a serious test of the proposed simulation approach. For a formal proof of the suitability of the amputation rules in our simulation, one should compare Eqs. (59) and (76). After multiplying Eq. (59) with the proper factors to obtain the effective collision amplitude (62), we see the claimed equivalence after using the definitions (29) and (30) for the electron and positron projectors. This reformulation demonstrates nicely that perturbation theory focuses on propagators whereas our simulations focus on collision rules and collision amplitudes.
For our actual simulations of the various terms in Eq. (76), we choose the convenient momentum vectors and we hence have The corresponding polarization vectors n α q defined in Eqs. (3)-(6) for α = 0, 1, 2, 3 are given by respectively.
To identify the form factors F and G, we have listed the two building blocks occurring in the representation (60) for our choice (77) of the momenta p and p ′ in Tables  I and II. According to Eq. (79), the rows for the spatial components µ = 1 and µ = 3 correspond to the physically observable transverse photon polarizations. For symmetry reasons, we can focus on the two double-underlined components in these tables. p ′ γ µ u σ p for the momenta given in Eq. (77); the pairs of spin values σ, σ ′ are given in the first row, the values of the spacetime index µ in the first column; our data analysis is focused on the two double-underlined components  77); the pairs of spin values σ, σ ′ are given in the first row, the values of the spacetime index µ in the first column; our data analysis is focused on the two double-underlined components 2m(Ep +m) 2m(Ep +m) By extrapolating to x = 0, we then obtain the g-factor for the electron according to Eq. (63).
To gain a better intuition why the deviation of the g-factor from 2 is given by G(0), one should look at Table III. The observation " Table I = Table II − Table III" is known as the Gordon decomposition. By means of this decomposition, one can rewrite the representation (60) of the vertex functions in terms of the momentum-and p for the momenta given in Eq. (77); the pairs of spin values σ, σ ′ are given in the first row, the values of the spacetime index µ in the first column; our data analysis is focused on the two double-underlined components spin-dependent terms in Tables II and III, so that the occurrence of the magnetic moment becomes plausible. For our toy simulation of H σ ′ σµ p ′ p , the summations over spin values and spatial components in Eq. (76) can actually be performed in a deterministic manner rather than by Monte Carlo methods. Only the integrations over the photon momentum q ′ in Eq. (76) need to be done by Monte Carlo simulations. We have tried the sampling procedures (71) and (73) for the photon momentum; it turned out that, for our simulation of the form factor G, the error bars for the procedure (71) are about two times larger than for the procedure (73), so that we adopted the latter one. To reproduce the symmetry properties in Tables I-III rigorously, we have actually symmetrized the integrands in Eq. (76) in the two-component of q ′ . All quantities have been made dimensionless by setting m = 1. The full simulation code consists of about a hundred commands.
Our simulation results for the form factor of the electron for y = 2 are shown in Figure 5. According to Eq. (80), the result for G(x) should be independent of our choice of y; an intermediate value of y seems to be most natural for the simulations (we have checked that y = 5 works equally well). For each data point, 40 million photon momenta have been sampled for the integrations. The required computer time on a single Intel Xeon 2.6 GHz processor is about two minutes per data point. Curves for q min = 0 and three different values of q max are shown Figure 5. For small x, each of the curves is linear in x 2 and can easily be extrapolated to x = 0. The extrapolated results are collected in Table IV. Larger q max leads to larger error bars. One should hence keep the cutoff as small as possible and then extrapolate to infinite cutoff, as shown in Figure 6. Our final extrapolated simulation result for infinite cutoff corresponds to in perfect agreement with Eq. (65). Of course, we could easily achieve higher precision by investing more than two minutes of CPU time per data point. However, the main point of our toy simulation is to demonstrate that a short and simple simulation program can reproduce the anomalous magnetic moment of the electron with minor computational effort. In view of our unsophisticated Monte Carlo sampling procedure, the efficiency of the simulations is truly amazing.
With the help of symbolic computation (series expansions and integrations handled by Mathematica), the four  Table IV to infinite cutoff; the dashed line represents a linear fit to the data points. integrals in Eq. (76) in the limit of small x and y have actually been calculated in closed form by Martin Kröger [18], thus reproducing Schwinger's famous result for g given in Eq. (65). The resulting relative contributions to g − 2 of the four integrals in Eq. (76) are given by (47 + 7π + 34 ln 2)/30 ≈ 3.085 for absorption of a photon by an electron, (−21 − 7π + 34 ln 2)/30 ≈ −0.647 for absorption of a photon by a positron, and equal values of (1 − 17 ln 2)/15 ≈ −0.719 for pair production and pair annihilation. Simulations show that, for increasing y, the contributions from absorption of a photon by a positron, pair production and pair annihilation become less relevant (while the total result, apart from the magnitude of the error bars, is independent of y, as pointed out before). For large y, say y ≥ 100, the result for g is entirely dominated by the absorption of a photon by an electron, that is, by the first Feynman diagram in Fig. 2 or, in other words, by the sequence of collisions shown in Fig. 3.

VI. SUMMARY AND DISCUSSION
We have rewritten QED as a kinetic-theory-like evolution equation with complex collision amplitudes in the Fock space of electrons, positrons and photons. By using longitudinal and temporal in addition to transverse photons, the collision rules for photon emission/absorption and electron-positron pair production/annihilation become particularly simple. This specific procedure of introducing unphysical degrees of freedom corresponds to choosing a gauge. The Fourier transforms of timeordered correlation functions can then be evaluated by stochastic simulations of trajectories in Fock space consisting of sequences of such collision processes. Summations and integrations introduced by the interaction Hamiltonian are evaluated by Monte Carlo techniques. In the limit of the average time step of the simulations going to infinity, we recover perturbation theory. Even for fundamental theoretical developments, the kinetictheory-like formulation of QED with four types of photons may be advantageous.
An illustrative toy simulation used here to obtain the leading-order contribution to the anomalous magnetic moment of the electron is very close to perturbation theory, which is known as an extremely successful tool in QED. However, the purpose of the proposed stochastic particle simulations is not to obtain high-precision results or to compete with perturbation theory. The goal rather is to provide a flexible tool that can go beyond perturbation theory (by choosing a finite average time step).
In more typical simulations with many time finite steps one can, for example, study propagators for different values of the ultraviolet cutoff q max and track how the (bare) mass and charge need to be changed with q max to get consistent results for the frequency dependence of the propagator. With such a procedure, one can extract both the physical properties of the propagators and the renormalization-group flow of the model parameters (mass and charge). Additional physical properties, such as form factors, can also be obtained by such simulations.
Even in QED there are questions of fundamental importance that cannot be conclusively answered by perturbation theory. Most important is the question whether the bare charge diverges at a finite ultraviolet cutoff, at the so-called Landau pole. Such a pole would intrinsically indicate a limitation of QED at high energies and would imply a vanishing coupling constant in the absence of an ultraviolet cutoff (that is, triviality of pure QED). This question has been addressed by extensive lattice simulations of QED (see, for example, [19][20][21][22][23][24][25]), which reveal a second-order chiral phase transition with a nonperturbative critical point in the strong-coupling regime (which does not correspond to a continuum ver-sion of QED). The lattice simulations imply that QED is a valid theory only for small renormalized charges [23]. Full consistency of QED can only be obtained in the bigger settings of electroweak interactions or the standard model. In view of these severe limitations of pure QED, it will be interesting to verify the consistency of the insights obtained from field-and particle-based simulations in the low-energy regime. In particular, one should compare the renormalization-group flow of the bare mass and charge with the ultraviolet cutoff. A major advantage of our particle-based simulations is that significantly larger variations of the cutoff can be achieved than with fourdimensional lattice simulations.
A regularization mechanism through an ultraviolet cutoff in momentum space is significantly different from a lattice formulation. Consistent simulation results in the low-energy domain of QED would hence be encouraging; although they cannot put QED on solid ground as a stand-alone quantum field theory, they would explain the enormous success of perturbative QED. For our stochastic simulation technique, still another regularization mechanism would be very natural: a rapid dissipative smoothing of the fields on short length scales [26]. Such a dynamic coarse-graining approach to quantum field theory actually provided the original motivation for kinetic-theory-like stochastic simulation technique developed in the present work.
A generalization of the proposed simulation techniques from QED to nonabelian gauge theories would clearly be desirable. For quantum chromodynamics (QCD), a naive generalization of the ideas of Gupta [4] and Bleuler [5] would suggest to introduce longitudinal and temporal gluons leading to particularly simple gluon propagators. However, additional ghost particles related to the Fadeev-Popov method for treating nonabelian gauge theories in the path-integral approach are required. The generalization of the ideas of Gupta and Bleuler to obtain a manifestly covariant canonical quantization of Yang-Mills theories has actually been developed in a series of papers by Kugo and Ojima [27][28][29][30]. As a key ingredient to their approach, Kugo and Ojima use the conserved charges generating BRST transformations [31,32] to characterize the physical states. In addition to the ghost particles, three-and four-gluon collisions make the kinetic theory for QCD more complicated than for QED. For electroweak interactions, an additional complication arises: the Higgs particle needs to be included into the kinetic-theory-like description. Although the details are considerably more complicated, the ideas of the present paper can be generalized to simulate the Yang-Mills theories for QCD and electroweak interactions and hence all parts of the standard model.