Quantum radiation reaction spectrum of electrons in plane waves

In previous works we derived equations for the average momentum of high-energy electrons experiencing quantum radiation reaction (RR) in strong electromagnetic plane-wave background fields. In this paper we derive similar equations for the momentum spectrum. We formulate the equations in terms of the cumulative function and study the relation between the equations for the spectrum and the equations for the moments, analyze the structure of the low-energy expansions, and finally explain how our formulation is essentially in terms of a ``Green's function'' which allows us to study the dynamics without choosing a specific initial wave packet (or particle-bunch distribution).

are always calculating some inclusive probability.Here we have the inclusive probability that an electron ends up with a certain momentum and spin given some initial momentum and spin, summed over the probabilities that this happens after emitting 0, 1, 2... photons.The probability that n photons are emitted is itself given by an infinite sum over loops.Clearly, one cannot calculate multi-photon emissions or higher-order loops exactly, even after simplifying by assuming the background field is a plane wave (for which the solution to the Dirac equation is particularly simple).In [32][33][34] we showed how higher-orders-in-α processes can be approximated by incoherent products of "strong-field-QED Mueller matrices".In [22,23] we showed how to sum the individual probabilities and how to actually evaluate the sum.The results are recursive and integrodifferential matrix equations.To summarize some of the differences between ours and other approaches, we do not introduce rates, our equations are formulated in terms of a single electron rather than an electron bunch described by a classical particle distribution, and we are not restricted to the locally-constant-field (LCF) regime.
As a motivation for considering plane waves, note that an electron with sufficiently high energy will effectively see a more general background field as if it were a plane wave, which follows from a short Lorentz boost argument.There are of course exceptions to this.For example, if, on the trajectory of the electron, the only nonzero component of the background field is an electric field parallel to the electron momentum, see e.g.[35,36].But in general one can expect a plane wave to be a good approximation of a general field.
In the literature it is common to simplify further by treating the plane wave in a LCF approximation.This additional step requires that 1 a 0 = E/ω is sufficiently large, where E is the maximum field strength and ω a characteristic frequency scale.In Sec.II and IV we derive general equations that are valid beyond the LCF regime, see in particular e.g.(38), (49), (121) and (129).The general Mueller matrices can also be expressed compactly in terms of known special functions if, instead of 1 We use units with the electron mass me = 1 in addition to c = ℏ = 1.

arXiv:2310.09863v1 [hep-ph] 15 Oct 2023
assuming large a 0 , one assumes that the field has circular polarization and is sufficiently long, which allows one to use a locally monochromatic field approximation, see [34].For different ways of treating the field as locally monochromatic, see [37,38].
Another important parameter is χ = −(F µν p ν ) 2 , where p µ is the electron momentum.Quantum effects can be neglected for sufficiently small χ.Here we are interested in values of χ that are large enough for significant quantum effects in RR, but small enough so that we can neglect pair production.To lowest order the probability of trident pair production scales as [39,40] P(e − → e − e − e + ) ∝ exp(−16/[3χ]).In [22,23] we studied the χ ≪ 1 expansions of the average momentum and the spin transition probability, i.e. the first two moments of the spectrum, and showed that they are asymptotic and can be resummed with the Borel-Padé method.In Sec.III we show how to obtain χ ≪ 1 expansions of the spectrum.We find that the expansion parameter is √ χ, so there is significant room for quantum effects in RR while pair production is still negligible.

II. DERIVATION
Lightfront coordinates are defined as and v ⊥ = {v 1 , v 2 }, so that the background field is given by a ⊥ (ϕ) and a ± = 0, where ϕ = kx = ωx + .We consider general pulse shape and polarization.We are interested in the dependence on the lightfront longitudinal momentum, kP .Rather than considering the spectrum directly, we consider instead a partially integrated spectrum or a cumulative distribution function, which we define as the probability that an electron which initially has longitudinal momentum kp emerges, after interaction with the background field, with momentum where 0 < x < 1.We also consider general spin transition.If we sum over the final spin and set x = 1 then we obtain a probability P = 1 as we should.x = 0 gives the probability of electrons that have lost no or almost no longitudinal momentum.Differentiating the final result with respect to x gives the spectrum.As our starting point is the Furry-picture expansion in α, we initially have the probability as a Taylor expansion in α, where each ) is a nontrivial function of the field strength 2 .We can write each term as a multiplication of the (4D) Stokes vectors, N 0 and N f , for the 2 Recall that we have absorbed a factor of e as eE → E.
initial and final spin and a (4 × 4) Mueller matrix, Note that we first calculate M n or so we do not need to choose any specific initial or final spin until the very end of the calculation, where we simply have to project the result for M with the Stokes vectors, The initial state is described by a wave packet as where is the usual Lorentz-invariant integration measure, written here in lightfront coordinates, and We consider first a wave packet which is sharply peaked, and then in Sec.IV we show that the results for a sharply peaked wave packet essentially give what can be thought of as a Green's function, from which one can afterwards obtain the results for an arbitrary, wide wave packet.A general spin can be written as where ρ and λ are two real constants.To zeroth order we hence find where the Stokes vectors are related to the angles in (10) as Here we have assumed that the wave packet is sufficiently narrow compared the values of x we consider.If we were to consider x too close to 0, then the step function would essentially be θ(kp ′ − kp) and we would find as only about one half of the peak of f would be integrated over.For example, if we were to consider a Gaussian wave packet proportional to then for sufficiently small λ we would find This is nonzero but exponentially suppressed for x < 0. A nonzero λ gives a regularized step function for the cumulative distribution, and a regularized delta function for the spectrum.However, until Sec.IV, we will assume that the wave packet is sufficiently narrow so that we can always approximate Thus, from (11) we have At O(α) we have the loop correction to scattering without emission and emission of one photon without loops, which can be written as M (1) (b 0 , x, −∞), where [33,34] M (1) where σ = (ϕ 2 +ϕ 1 )/2 is the average lightfront time, with ϕ 2 being the lightfront-time variable for the amplitude M and ϕ 1 for the complex conjugate M , q = kl/kp = kl/b 0 is the ratio of the longitudinal momentum of the photon and the initial electron, and M C and M L are the Mueller matrices for photon emission and the electron self-energy loop.The reason for introducing the lower integration limit for the lightfront-time integral, θ(σ ′ − σ), will be explained below.The step function θ(x − q) is due to the restriction in (2), which to this order and for this term can be rearranged into an upper cut-off for the photon momentum, q < x.There is no such step function for M L because that q integral corresponds to a photon that is emitted and reabsorbed and therefore does not change the final electron momentum.
For an arbitrary field, M C is given by Eq. ( 24), ( 25), ( 26), ( 27) and ( 29) in [33], and M L by Eq. ( 67) to (71) in [34], which we write here as 3 where θ = ϕ 2 −ϕ 1 , the integration contour goes above the pole at θ = 0, r = (1/s) − 1, s = 1 − q is the ratio of the longitudinal momentum of the electron after and before emitting the photon, b 0 = kp, and M 2  21 is an effective mass where For photon emission, the 4 × 4 matrix is given by where and where and 3 When comparing the overall factor of 2 in Eq. ( 24) in [33], note that here we get one extra factor of 2 when summing over the photon polarization, and one factor of 1/2 has been factored out, writing Note that no special notation has been used for outer products, so, for example, ( k X•V) j = kj (X•V).For the loop we have so some of the elements are identical to R C , while where Y = w 2 − w 1 , gives spin rotation.Note that where so if we sum over the final spin state and if we integrate over all momenta, i.e. x = 1, then P (1) = 0 for any initial spin state.We also have P (n>1) = 0, so P = P (1) = N 0 •e 0 = 1, which is what it has to be due to unitarity.
To O(α 2 ) we have where, in the last term, we have used where (in this term) kp ′ = kp − kl 1 − kl 2 = kp 1 − kl 2 , q 1 = kl 1 /kp, q 2 = kl 2 /kp 1 , and kp 1 /kp = 1 − q 1 .Note that q j is the longitudinal momentum of the photon emitted (or emitted and reabsorbed) at step n divided by the longitudinal momentum of the electron immediately before this step.Writing the step function as in (35) allows us to obtain M (2) from M (1) by prepending M L and M C as •M (1) Note that if we replace M (2) → M (1) and M (1) → M (0) = 1 then we have (18).Higher orders can be obtained recursively 4 from Summing over n ≥ 1 and differentiating with respect to σ gives where we have suppressed the lightfront time argument σ as it is now the same in all terms.We have an "initial" condition at σ → +∞ rather than σ → −∞, and then we integrate backwards in lightfront time, where the final result is given by M(σ → −∞).
An alternative form that might be illuminating is obtained by using b ′ := (1 − x)b 0 instead of x for the final momentum, and by changing integration variable from q to b = (1 − q)b 0 , As a check, for x = 1, which means no restriction on the final momentum, we have which is the same equation as in [22,23].For x = 0, which means observing only electrons that have lost no or very little longitudinal momentum, the photon-emission term becomes negligible due to θ(x − q), and we have The solution is given by a lightfront-time-ordered exponential, where T means anti-lightfront-time ordering.This agrees with the result in [34].As shown in [34], for a field with linear or circular polarization, one can write the result in an explicit form without a time-ordered exponential of a matrix.Thus, in general, the new equation (38) interpolates between the result in [34] for x = 0 and the one in [22,23] for x = 1.

A. Moments
The spectrum is given by from which we can obtain the moments as The lower integration limit is not x = 0 because then we would miss the essentially delta-function-like peak at x = 0.The width of the wave packet determines the width of this peak.For a sharply peaked wave packet the probability for x < 0 is very small but it is nonzero.We can avoid the x < 0 region by making a partial integration, where in the second step we have neglected the x < 0 part since now the integrand has no delta-function-like peak and is negligibly small for x < 0 (because we have assumed a sharply peaked wave packet).We can integrate (38) over x before solving it, i.e. we can express entirely in terms of M. This is trivial for the first two terms, while for third we change variable from x to x ′ = (x − q)/(1 − q).We find We thus have the same equation for all n.n only appears in the "initial" condition Note that (49) allows us to obtain each moment separately.In other words, if we want, say, the second moment then we do not need to consider any other moments and we do not need to obtain the spectrum.For n = 0 we have M(n = 0) = M(x = 1) and so ( 49) is the same equation as (41).n = 1 gives the expectation value ⟨kP ⟩ which we studied in [22].(49) together with (50) are what one should expect by a generalization of the derivation for n = 1 in [22], but this indirect derivation of ( 49) from ( 38) serves as a check of (47).
Apart from checking that these two sets of equations are consistent, this also allows us to check the numerical results obtained from them.After having obtained a numerical solution of M(x) (for x > 0) we can check the result by integrating it as in (47) for the first couple of moments and comparing with the moments obtained by instead using (49), which is much faster to solve since M only has two integration variables, σ and b 0 , while M also has x.
Thus, all the moments obtained from M(b 0 , x) by solving the new equation (38) and using (47) agree with the moments obtained using the approach in [22].Two different distributions can in principle have the same moments 5 .However, in the above derivation we did not need to assume that n is an integer.If we let n be a continuous variable then we essentially have the Mellin transform of the spectrum.The inverse would be given by an integral over n in the complex plane.Or we could replace (1 − x) n b n 0 in (46) with some arbitrary function f ([1 − x]b 0 ).If f (0) is nonzero then we can simply deal with that constant separately, so we can assume without loss of generality that f (0) = 0. Now the resulting M is determined by the same equation as the moments, i.e. ( 49), but with "initial" condition Thus, the result is again what we would find if we instead started with the approach in [22].

B. Locally constant field approximation
For sufficiently large a 0 we can use a LCF approximation.Here the θ integrals in (19) can be performed in terms of Airy functions, Ai, Ai ′ and Ai and the Scorer function Gi, (52) From [34] we have and where ξ = (r/χ(σ) the local directions of the electric and magnetic are defined initially as 3D vectors as but made 4D trivially as and similar for B(σ), k = e 3 , 1 4 = 4 j=1 e j e j and The e 0 B and Be 0 terms lead to induced polarization along the magnetic field direction, as a Sokolov-Ternov effect, while the k Ê − Êk term gives spin rotation, including the effect of the anomalous magnetic moment (cf.[42,43]).For a rotating field, we have in general B(σ 1 ) • Ê(σ 2 ) ̸ = 0, so all 4 components couple, and an induced polarization along the local direction of the magnetic field will later be along the electric field direction, which undergoes spin rotation.However, for a linearly polarized field, the e 3 and Ê components decouple from the e 0 and B components.When we for such fields consider initial or final Stokes vectors with e 3 • N 0,f = Ê3 • N 0,f = 0, then we drop the irrelevant components and use 2D Stokes vectors and 2 × 2 Mueller matrices.

III. CONSTANT FIELD
For a constant field we have so we have an effective expansion parameter We separate T n from M (n) , changing notation slightly so that We also separate a factor of αa 0 from M L,C , and then the recursive formula simplifies (62) We can sum this into an integrodifferential equation where the "time" variable is T instead of σ used in (38), with "initial" condition M(T = 0) = M (0) .In contrast to (38), where the physical result is only obtained by, in the end, setting σ → −∞, all values of T in (63) give physical results, and if we are interested in a field with e.g.T = 10 then we obtain the results for T < 10 as a byproduct because we integrate the equation starting with the initial condition at T = 0.
The results for χ = 0.01 are shown in Fig. 1.For these results we have used a step size of ∆T = 0.1.This is smaller than what it might first seem because the natural O(1) "time" variable is actually χT (see (70)).Fig. 1 shows the results for every tenth time step (or every twentieth if one counts the midpoints).For χ we have used an evenly distributed grid from χ = 0 to χ = χ max = 0.01 with ∆χ = χ max /100.As can be seen in the plot for M 11 in Fig. 1, ∂M/∂x diverges at x → 0 for small T .Recall that this is because we have assumed a sharply peaked wave packet and is why we work with the cumulative distribution.While M is finite and all the integrals converge, we want to avoid having to use a large number of points in x in order to obtain a good interpolation near x = 0, because that would make the code much slower.For small T we therefore use an evenly distributed grid in y = x 1/3 instead of x, which is better because ∂M/∂y does not diverge at y = 0 and is easier to interpolate using fewer points.We can illustrate this by looking at the part of the integral in (63) with q > x, Since this is independent of M, we have made a separate interpolation of it so that we do not have to calculate it over and over again when solving (63).that, as a function of x, ∂J/∂x diverges as x → 0, but ∂J/∂y is finite and therefore easier to interpolate.We have used ∆y = 0.01.As T increases beyond a certain finite point, M becomes exponentially small near x = 0, and then we no longer have the problem of a divergent ∂M/∂x.When this happens we switch to using an even grid in x instead of y, with ∆x = 0.01.Thus at each point in T we make a cubic-polynomial interpolation of the list of M(χ, x) evaluated on a 100 × 100 size grid.
The results for S for χ = 0.1 are shown in Fig. 3.The spectrum is initially a delta-function-like peak (because we hav neglected the width of the wave packet), then becomes a relatively low and wide Gaussian peak, but as T increases further the peak eventually starts to become higher and narrower as it is squeezed towards the upper limit x = 1.This evolution of {1, 0} • S • {1, 0} is similar to the time evolution of distribution functions found in Fig. 1 in [1] and Fig. 2 in [44], even though both methods and setups are different, e.g.we consider a single electron with a wave packet that is initially sharply peaked, while [1,44] considered electron beam distributions and neglected spin.as in (47) for every tenth time step.The corresponding ∂xM is shown in Fig. 3.

A. Moments
The moments 6 ⟨(a 0 kP ) m ⟩ can be obtained by solving with "initial" condition The derivation of this equation is a trivial generalization of the m = 0 and m = 1 cases in [22,23].Note that with this approach we obtain the moments without calculating the spectrum.We thus have two independent ways of calculating the moments: either by solving (65) or by first obtaining the spectrum via (63) and then the moments by integrating as in (47).
We have solved (65) numerically using the midpoint method for the variable T , as described in [23].At each step in T we make an interpolation function for the χ dependence, with 0 < χ < χ max , where χ max is the maximum value of χ we want to consider.We have used an adaptive grid for χ, i.e. adding points until a certain precision and accuracy is reached, which results in more points where the function has a larger curvature and less points where it is flat.The interpolation between these points are done using Mathematica's built-in function "Interpolation" with cubic polynomials. 6We have included a factor of a 0 here so that a 0 only appears in χ or T .
M in (65) already gives the resummation of all orders in α.We can also obtain the moments by first calculating each order in α separately by solving the recursive formula starting with and then resum the α expansion at the end.Also in this case we can make a numerical interpolation function for the χ dependence for each step.After having obtained terms up to e.g.n = 10 or n = 20, we can resum the α expansion using Padé approximants as described in [22,23].The convergence of this resummation method is often fast.Alternatively, we can expand each order in α as a power series in χ.Then the χ expansion of Mn is obtained by inserting the χ expansion of Mn−1 into (67) and expanding the result in χ.After having obtained χ expansions for the first e.g.n = 10 orders in the α expansion, we first resum the χ expansions for each order in α separately using standard Borel-Padé resummation, and finally we resum the α expansion using Padé approximants (this last step is the same with or without using the χ expansion).A detailed explanation for this doubleresummation approach can be found in [22,23].It turns out to be quite fast.
As an example, we consider χ = 0.1.The zeroth and first moments look very similar to the results already presented in [22,23].Calculating also the second moment allows us to obtain the standard deviation, In Fig. 4 we plot S obtained by using either the double resummation approach or by integrating the cumulative function as in (47).We find perfect agreement.Comparing with the double-resummation approach serves as a quick way to check the precision of M(x), which takes much, much longer time to obtain.The shape of S is quite similar to plots in Fig. 4 in [44], Fig. 3 in [45] or Fig. 9 in [6] for the energy spread of an electron beam/bunch, i.e. first a rapid increase and then a slower decrease after a maximum.But, again, it should be noted that we are considering somewhat different quantities.

B. Low energy expansion of moments
We will now calculate the first quantum correction in a low energy expansion.χ ≪ 1 is the expansion parameter, but we need to consider as O(1) to keep a nontrivial dependence on T 7 .The zeroth order is then given by (71) This is just the longitudinal-momentum component of the solution to LL [47] to the power of m.The standard deviation (69) vanishes at O(χ 0 ), as it has to in the classical limit.To obtain the first nonzero correction we can solve the recursive equation (67) approximately using the ansatz where δ M only depends on m but not on χ.With the problem has been reduced to an algebraic recursive equation for δ M(n) as function of n.There are standard methods to solve such recursive equations [48], but this is most conveniently done using Mathematica's "RSolve".The result for δ M(n) is not particularly illuminating.But one thing to note though is that it is an expression valid for all orders in α, which we can obtain here because we only consider the first two orders in the χ expansion.Contrast this with the general case described above where we would only be able to calculate the first e.g. 10 or 20 orders and then use e.g.Padé resummation based on those finite number of terms.Now for the low-energy approximation, we have access to all coefficients in the α expansion and we can resum this using Mathematica, for example.For the average momentum, i.e. m = 1, we find We recognize the {0, 1}•δ M•{1, 0} element from Eq. ( 13) in [22], which corresponds to the difference in the final momentum due to initial spin being either parallel or antiparallel to the magnetic field and after summing over the final spins.For m = 2 we find (75) 7 There can however be experimental signals in the spectrum of the emitted photons even if u ≪ 1 [46].
If we sum over the final spin then we find a standard deviation that to leading order does not depend on the initial spin, Thus, the width of the spectrum goes to zero at both u ≪ 1 and u ≫ 1, which is also what we see in Fig. 3.
If we expect the spectrum to have a more or less symmetric peak around the average momentum, then it is not surprising that S → 0 asymptotically, because the average momentum decreases as 1/(1 + u) and the lightfront longitudinal momentum kP > 0, so the lower half of the peak will be squeezed between 1/(1+u) and 0, and hence the peak must become narrower.However, even if we normalize the width by dividing S by 1/(1 + u) the result still goes to zero.In [22] we showed that the average momentum converges to the classical LL momentum asymptotically, and now we can also see that the standard deviation decreases asymptotically.We also note that the width scales as S ∼ √ χ.So, if, say, χ = 0.01 then we would expect S ∼ 0.1, which corresponds to 10% of the physical interval for this scaled momentum variable (i.e.0 to 1).Thus, even if χ is so small that there is little difference between the average momentum and its classical limit (the solution to LL), the width can still be a significant fraction.While it would be possible to extend the above approach to include both higher powers of χ as well as higher moments, it becomes inconvenient to solve the recursive equations for the expansions in u and then resumming them, so we have instead used the following approach for higher orders.We change variable in (65) from T to u in (70).We separate out the overall factor of χ m in (72) as The integrodifferential equation in terms of W is given by with the same "initial" condition for all moments.Note that, while (65) is local in T , (78) is not local in u.The reason for making this change of variables is to have a natural χ expansion right from the start and to avoid expanding in u.This expansion can now be written as To find the first couple of w's, we insert (80) into (78) and match the two sides order by order in χ.To do so we need a suitable change of variable for the q integral.We have used γ as defined by After this change of variable, we can expand the integrand before performing the integral.This results in the following type of integrals, where the second line shows the explicit numbers for n = 0, 1, 2 . . ., These are the same integrals that we used in [22] to calculate the χ expansion of each of the orders in α separately.Now we instead work with (80), which is already resummed in α (recall u ∝ α).
To zeroth order we find which implies in agreement with ( 72) and what one should expect from the solution to LL. w m,k (u) starts contributing at O(χ k ) in the expansion of (78).At this order, w m,0 (u) to w m,k−1 (u) also contribute.Moving all the w m,k (u) terms to the left-hand side gives a differential equation on the following form We find that w m,k (u) can be expressed in terms of linear combinations of where r and s are integers.Derivatives of L r,s appear when we replace u → (1 − q)u, change variable to γ and expand in χ, but so this class of functions is closed under differentiation.We can integrate (87) using and so this class of functions is also closed under integration.For m = 0 we find which agrees with Eq. ( 10) in [23], and for m = 1 and m = 2 we recover (74) and (75).At the next order in χ we find, for example, As we go to even higher orders we get more and more terms, and the expressions become too large to write down here.We will use them, though, in the next section to determine integration constants in the derivation of an analytical approximation of the spectrum.

C. Low energy approximation of the spectrum
To obtain an analytical approximation for the spectrum, we start by differentiating (63) with respect to x to obtain an equation directly expressed in terms of S, The last term makes this an inhomogeneous equation.It is given by the time-ordered exponential in (43).However, for low energies it is exponentially suppressed and can be neglected.
We know from (76) that the width of the peak in the spectrum is O( √ χ).In the limit χ → 0 this would be an increasingly narrow peak, centered at a point which we denote x = x cl (u), where we again use u in (70).In the classical limit, we have ⟨kP ⟩/kp → 1/(1 + u) (the solution to LL), so from (2) we have If we were to plot S(x) for very small χ, it would be natural to not plot it over the entire interval 0 < x < 1, because then we would just see a very sharp peak.Instead we would plot it over where c 1 and c 2 are roughly ∼ 5 − 10 or so.Beyond this interval S(x) would anyway be negligible small.Similarly, in order to obtain an analytic approximation, we also want to figure out how to choose parameters such that they can be considered O(χ 0 ).Here we keep u as "time" variable, but switch to as a momentum parameter.λ(u) is related to the standard deviation and will be determined below.As an educated guess, and after some trial and error, we take as an ansatz The idea now is to insert this expansion into (94) to determine the functions λ, E and ρ.We again change integration variable from q to γ as in (81), which allows us to expand the integrand in a series in χ before performing the integral.This expansion is done with u and X considered as O(1) parameters.In particular, for the M C • S term, the replacements χ → (1 − q)χ and x → (x − q)/(1 − q) lead to With and we see that the expansion of ∂S/∂T and hence also the right-hand side of (94) starts at O(χ 0 ).Rearranging (94) as RHS−LHS= 0, we find at O(χ 0 ) With initial condition x cl (0) = 0, (101) implies that x cl is given by (95).
At O(χ 1/2 ) we would in general get derivatives on ρ 0 , but as part of the ansatz we take ρ 0 = 1.We then find This equation should hold for any values of u and X, and λ (E) should only depend on u (X).If we first set X = 0, then we obtain a differential equation for λ, which we solve with initial condition λ(0) = 0, which is motivated by the fact that the standard deviation should go to zero at u = 0, as we found in (76).We find Inserting ( 103) into (102) gives and equation that only involves X, (104) The solution to this equation is At first it might seem like we have two integration constants to determine, E(0) and E ′ (0).However, from ( 96) and ( 103) we see that the constants in the exponent actually cancel.This is not surprising since we could have included any O(1) constant in the definition of X.We can therefore set E ′ (0) = −E(0) so that the exponent is simply e −X 2 .The remaining constant is determined by which to leading order implies Thus, expanding (94) to next-to-leading order (i.e.O(χ 1/2 ) has allowed us to determined S to leading order, and we find that the width of the peak is determined by and the shape of the peak by However, since the expansion parameter in (97) is only √ χ rather than e.g.χ, we expect the next couple orders to be important to obtain a precise approximation even for χ ∼ 0.01, so we continue.
At O(χ) we find a differential equation for ρ 1 (u, X).This equation contains terms that are linear in ρ 1 (u, X) and derivatives of ρ 1 (u, X) with respect to u and X, and terms without ρ 1 (u, X).After dividing away the overall factor of e −X 2 , the terms without ρ 1 (u, X) can be written f 1 (u)X + f 3 (u)X 3 .The solution is therefore on the form The equation for ρ 1 thus separates into one part proportional to X and another part proportional to X 3 .Both parts have to be separately zero.Solving these equations gives and where c i are some constants.We cannot determine these constants by considering some initial conditions at u = 0, because our approximation breaks down when u becomes too small.We can see this from Fig. 1, where the x derivative diverges at x = 0 for T below some finite point.A diverging S is of course not approximated by (97).To determine the constants c i we will instead compare with the moments calculated with the method in the previous section.
We will use the same method also for higher orders, and we have found that ρ n (u, X) is in general a polyno-mial in X.To compare with the moments we consider with w defined in (80).By matching the first and the last lines for each m and for each order in χ we obtain equations that we can use to determine the constants in (111) and (112), and other constants at higher orders.
By evaluating (113) with m = 1 and m = 3 and selecting the term proportional to χ 1/2 we are able to determine the constants in (111) and ( 112 In Fig. 5 we compare the approximation in (97) with a numerical solution of (63).We have chosen χ = 0.01, which at first might seem like a quite small value.However, the approximation is a series in √ χ = 0.1, so we should not expect higher orders to be negligible.In Fig. 5 we see that, while the leading order gives a decent approximation, adding higher orders does indeed lead to a noticeable difference.Even after adding the next-toleading order there is still a noticeable difference from the numerical result.For {1, 0} • S • {1, 0}, it is only after adding the first three terms that we obtain a result that is more or less indistinguishable from the numerical solution.For {0, 1}•S•{1, 0} we need the first four terms.
We have mostly considered the cumulative function as a computationally convenient tool for in the end finding the spectrum, but to obtain this low-energy approximation of the spectrum we worked directly with S. If we want M it is now straightforward to integrate (97).We have found that the ρ functions are in general polynomials in X. Regardless of what sort of polynomial they are, we can always write them as sums of Hermite polynomials, H m (X).We can therefore write (97) as This is useful because integrating over X becomes trivial using Rodrigues' formula, and we find IV. FINITE WAVE PACKETS So far we have focused on sharply peaked wave packets.Now we consider wave packets with finite width.In principle we could consider a wave packet with two different functions, f 1 (p) and f 2 (p), for two different spin states, but we will for simplicity consider just a single function f .When we considered a wave packet sharply peaked at b 0 , we found it natural to factor out b 0 as in (2).But now when we have an integral over the initial momentum, we replace (1 − x)b 0 → b ′ , so that we can describe the final momentum without referring to the initial momentum.To avoid confusion of the integration variable kp ′ with b ′ we also replace as kp ′ → kp out .To generalize we first note that each contribution to the cumulative function involves momentum integrals on the form where in the last line the sum is over all real photons emitted in this particular term, and the factor of 1/k + is just due to our normalization of the amplitude M .In this paper we only consider observables that do not depend on the transverse momenta.After integrating over the transverse momenta of the final-state particles, the probabilities in plane waves no longer depend on the transverse momentum of the initial particle.We write where ρ now describes the initial longitudinal momentum distribution.We thus find a simple relation between the Mueller matrix M(f ; b ′ ) for a wide wave packet and M(b 0 , b ′ ) for a sharply peaked wave packet9 , where M(b 0 , b ′ ) is determined by (40).Note that even if we were only interested in a sharply peaked wave packet for only one value of b 0 = B, to solve (40) we anyway need to consider M(b 0 , b ′ ) for 0 < b 0 < B. Thus, once we have solved (40) for a sharply peaked wave packet, there is very little extra work to perform the integral in (121) for various types of wave packets.In other words, all the really nontrivial stuff is included in M(b 0 , b ′ ), which is obtained without having to choose spin states or wave packets.We can therefore think of M(b 0 , b ′ ) as a Green's function.
Partly in order to compare with kinetic equations, we will rewrite this in terms of the spectrum (this differ by a factor of b 0 compared to the previous definition (44)) We begin by re-expanding as We obtain higher orders from ( 18) and ( 34) by replacing and writing the q j integrals in terms of kl j .For example, in the last term in (34) we have where where the first argument is for the momentum of the electron before emitting the photon.In (34) and (36) we noted that M (2) can be obtained by prepending M L and M C to M (1) .From (124) we see that S (2) can instead be obtained by appending M L and M C to S (1) .We can do this if we also replace Higher orders can be obtained in the same way.Thus, we find where Summing over n and differentiating with respect to σ gives with initial condition If one is only interested in a definite initial Stokes vector, then one can project N 0 • (129) before integrating and solve for N 0 • S with initial condition N 0 • S(−∞) = ρN 0 .On the other hand, (129) • N f does not give an equation for S • N f .For (40) we can instead project, before integration, with N f but not with N 0 .In contrast to the equations for cumulative function M and the moments, see e.g. ( 40), which we integrate backwards in lightfront time from σ = +∞ to −∞, (129) is integrated forwards from σ = −∞ to +∞.We can understand this difference as follows.For S(f ; b ′ ) the earliest step is special due to the appearance of the initial wave packet, while the subsequent steps are obtained recursively by appending M L and M C , which gives an equation that should be integrated forwards in σ.For the moments M(n, b 0 ) it is instead the last step that is special, because there we have an additional factor of kP n in the integral for the final electron, while the preceding steps are obtained recursively by prepending M L and M C and there is no nontrivial distribution in the first step, which gives equations that should be integrated backwards in σ.

V. CONCLUSIONS
We have derived recursive and integrodifferential matrix equations for the momentum spectrum for electrons in a plane wave.This is formulated in terms of what can be thought of as a Green's function, M(b 0 , b ′ ), which gives a complete description of RR and spin transition.After M(b 0 , b ′ ) has been calculated one obtains the result for particular initial and final electron states by projecting with the initial and final Stokes vectors for the spin, N 0 • M • N f , and integrating over the initial longitudinal momentum, b 0 , weighted by the absolute square of the wave packet, as in (121).
Due to the singularity at b ′ = b 0 during early ligthfront times, we have found it convenient to work with the cumulative distribution function, M(b 0 , b ′ ), rather than the spectrum itself.The spectrum is then obtained at the end of the calculation by differentiation, ∂ b ′ M(b 0 , b ′ ).
By integrating the resulting spectrum we find agreement with moments calculated using a generalization of the approach in [22,23], where we considered the zeroth and first order moments.With the second moment we can calculate the standard deviation, and we find that it scales as √ χ and is therefore relatively large even for small χ.
We have derived a low energy expansion of the spectrum for a constant field and found that it takes the form of a Gaussian multiplied by a power series.To zeroth order this tells us how the momentum is distributed around the solution to LL [47].We find that the expansion parameter is proportional to √ χ, so even for χ = 0.01 we need to go beyond the leading order and sum the first 3 or 4 orders in the χ ≪ 1 expansion in order to obtain a precise approximation.
For the numerical results of this paper, we have chosen a constant field.We do not actually expect it to be more difficult to solve (38) in the LCF approximation for a nonconstant field compared to (63) for a constant field.Indeed, both equations are on the form ∂M/∂t = F, where t = σ or t = T , and whether F only depends on t via M, as in the constant field case, or also has an explicit dependence on t, as in the nonconstant case, we could still use the same method (e.g. the midpoint method) when integrating over t.In fact, we did so in [51] for the all-order probability of trident.The issue is instead that for the spectrum we have an additional parameter, x, compared to the equations in [22,23,51], so it takes much longer time to compute an interpolation function of M(χ, x) at each point in "time".While this is an issue for both (38) and (63), the "time" parameter in (63) is a physically relevant parameter, so when we integrate from T = 0 to e.g.T = 100, all the intermediate time steps give physical results.In contrast, for (38) we integrate from σ = +∞ to σ = −∞, but it is only the final result at σ = −∞ that gives something physical.Moreover, in order to check the numerical results by comparing with analytical approximations, we have to generalize the lowenergy approximations to nonconstant fields.We expect the approximations for a constant field to serve as a useful guide for generalizing to nonconstant fields.Thus, we leave a numerical study of different nonconstant fields and generalization of the approximation of the momentum spectrum for future studies.
Moreover, our general methods are not restricted to the LCF regime, but work as long as the field is sufficiently strong or long.So, one could use for example a locally monochromatic field approximation, which we also leave for future studies.