Resummation of quantum radiation reaction in plane waves

We propose a new approach to obtain the momentum expectation value of an electron in a high-intensity laser, including multiple photon emissions and loops. We find a recursive formula that allows us to obtain the $\mathcal{O}(\alpha^n)$ term from $\mathcal{O}(\alpha^{n-1})$, which can also be expressed as an integro-differential equation. In the classical limit we obtain the solution to the Landau-Lifshitz equation to all orders. We show how spin-dependent quantum radiation reaction can be obtained by resumming both the energy expansion as well as the $\alpha$ expansion.

An electron in an electromagnetic field emits photons and the recoil it experiences is called radiation reaction (RR) [1,2]. Classically, the standard equation is the Abraham-Lorentz-Dirac (LAD) equation, but since it leads to unphysical solutions it is common to replace it with the Landau-Lifshitz (LL) equation. The first RR experiments with high-intensity lasers were performed recently [3,4] (see also [5]), and in upcoming experiments one may expect significant quantum effects [6][7][8][9][10]. In this paper we propose a new method for obtaining quantum RR in high-intensity lasers. It is based on incoherent products of O(α) terms (cf. [11]), with Mueller matrices taking spin into account (cf. [12,13]).
Let p µ and P µ be the initial and final momentum of the electron, k µ the wave vector of the plane wave, which is described by a gauge potential a ⊥ (σ) where σ = kx = ω(t + z) is the (rescaled) lightfront time and a ± = 0 1 . We will focus on kP = ∞ n=0 kP (n) , where kP (n) is O(α n ) and the dependence on the initial and final spin can be expressed as kP (n) = (1/2)N 0 ·M (n) ·N f , where N = {1, n} is the Stokes vector for spin along the unit vector n and M (n) are "strong-field-QED Mueller matrices". Averaging and summing over the initial and final spin gives {1, 0} · M · {1, 0}. We know from [14][15][16] how to obtain O(α n ) probabilities to leading order for long pulses or large a 0 = E/ω by multiplying O(α) Mueller matrices. In this case we need M C for (nonlinear) Compton scattering and M L for the cross term between the O(α 0 ) and O(α) parts of the amplitude for e − → e − . From this we now find the following recursive formula where b 0 = kp, q = kl/kp and l is the photon momentum, and M (0) = b 0 1 (M (0) = 1 for the probability). The lower integration limit has been introduced so that the products of Mueller matrices are lightfront-time ordered. The final result is obtained by setting σ = −∞. The This integro-differential matrix equation gives quantum RR to all orders in α with spin taken into account. Note that even if we are only interested in unpolarized initial and final states, we still need to solve a matrix equation before we can project with {1, 0} · M · {1, 0}. Note also that (2) holds even if a 0 is not large, provided the field is sufficiently long and the full version of M C,L is used. Eq. (2) looks similar to kinetic RR equations [9] (see also [17]), so one may expect that it can be solved with similar methods. However, in this paper we will instead use (1) to obtain the first M (n) terms before performing a resummation. We will first show that the classical limit of (1) agrees with the solution to LL. Classically there is no spin dependence, so we neglect the spin components and focus on M = {1, 0} · M · {1, 0}. Using M L = −M C and expanding to first order in the photon momentum, we find which is solved by where J is obtained by matching with M (1) , which we have already calculated in [18,19], Resumming this geometric series we finally find In A we show that the P ⊥,+ components can be calculated in a similar way. We find and P + satisfies the mass-shell condition P + = (P 2 ⊥ + 1)/(4P − ). This resummation 2 agrees exactly with the exact solution to LL [20] (see also [9]).
This does not mean that LAD does not agree with the classical limit of QED, because our approximation neglects terms that have subdominant scaling with respect to the pulse length and/or a 0 . From Eq. (4.39) in [19] we see that the difference between LL and LAD at O(α 2 ) is a term proportional to a 2 , i.e. a term without an integral. This term vanishes at asymptotic σ, but is also subdominant at finite times because it has no pulse-length scaling. We can therefore not rule out that the classical limit of all QED contributions may agree with LAD rather than LL.
In the limit of a very long pulse we have so P − becomes small, P ⊥ stays O(1) and P + becomes large, and (since ∆ ∝ p − ) all components of P µ become independent of the initial momentum.
The ratios N (n+1) /(χN (n) ) of neighboring coefficients in the α expansion as a function of χ, for an unpolarized initial state, N (n) = {1, 0} · N (n) . The ratios have been divided by χ in order to show the convergence to the classical limit Having shown that (1) gives the expected classical limit, we now turn to quantum RR. For simplicity we focus on P − , and we sum over the final spin, i.e. N f = {1, 0} is fixed and we have an overall factor of 2, so we use N (n) := M (n) · {1, 0}. In this first application we consider a constant field, so the σ integrals give dσ n ...dσ 1 → (∆φ) n /n!, and it is convenient to rescale N (n) so that where T = αa 0 ∆φ is our effective expansion parameter, which can be T > 1 for large a 0 ∆φ. In general N has four elements, but here we consider only initial and final states that are either unpolarized or polarized (anti-)parallel to the magnetic field, and then only two elements contribute. Omitting the irrelevant elements we have N 0 = {1, 0} and N 0 = {1, ±1} for unpolarized and polarized states. We have where N (0) (χ) = χ{1, 0}, 1/n comes from (∆φ) n /n!, χ = a 0 b 0 is the quantum nonlinearity parameter (gives the electric field in the rest frame of the electron), and where ξ = (r/χ) 2/3 with r = (1/s 1 ) − 1, κ = (1/s 1 ) + s 1 and s 1 = 1 − q.  In order to compute N (n) we have used two completely different methods. In the first we compute N (1) (χ) between χ = 0 and some 3 χ max and make an interpolation function of it, which is then used in (10) to compute an interpolation function of N (2) (χ), and so on.
In the second method we expand N (1) (χ) in a power series in χ, which is then used to obtain a power series of N (2) (χ) and so on, see B. This gives Fig. 1, the coefficients have alternating sign and grow factorially a m ∼ (−1) m m!. The χ expansion is therefore asymptotic with zero radius of convergence and needs to be resummed. There is no unique resummation method (unless, of course, one is able to find an exact expression). Different resummations are obtained by matching the se-  [23][24][25] (see [26] for further applications in strong-field QED). Such resummations require fewer terms, but usually some extra information e.g. about the opposite limit (large χ in our case). However, in the present case, we can without problems obtain a large number of terms quickly. We can hence use the general Borel-Padé resummation method [27][28][29][30][31][32][33][34][35][36], which only requires the a  Using either of these two approaches we obtain a set of functions, N (n) (χ), for 0 < χ < χ max . The result for an unpolarized initial state is shown in Fig. 2. In the low χ limit we find convergence towards the classical prediction.
In general one would also expect the α expansion to be asymptotic. However, an approximation does not have be asymptotic, see e.g. [31,[37][38][39][40] for the weak and strong field approximations of the QED effective action. (See [41,42] for other recent α resummations.) In classical RR, LAD leads to an asymptotic series [43], while LL has a finite radius of convergence. Since the coefficients we have calculated suggest a finite radius of convergence, we propose to resum the α expansion with a Padé approximant where A i and B i are obtained by matching a 0 kP re = χ + I+J n=0 T n N 0 · N (n) + O(T I+J+1 ). Since we expect a finite limit for T 1, we take I = J. This makes it possible for kP re → 0 as T → ∞, which is what we expect (cf. (8)). If we do not impose this limit, then we can take |A I /(χB I ) + 1| as an upper-limit estimate on the relative error (the error goes to zero as T → 0). Alternatively, we can demand kP re = O(1/T ) for T 1, which fixes A I = −χB I . LL predicts that the leading order in T 1 is independent of the initial momentum (8). If we assume that that holds in general, then the O(1/T ) term must be the same as the classical prediction, which implies Fig. 3 shows that the resummation (14) converges quickly. On the scale of this plot, the I = 4 and I = 5 resummations are virtually indistinguishable, which are obtained from the first eight and ten N (n) terms. For short to moderately large T , we see that the classical prediction overestimates the effect of RR, as expected [1]. However, for larger T the classical and quantum predictions seem to converge. This is what one would expect if the leading order at T 1 is independent on the initial momentum. This motivates us to use the modified Padé approximant based on (15). With the two extra conditions in (15) we indeed find an even faster convergence, with a decent approximation already with I = 2, i.e. using only the O(α) and O(α 2 ) terms.
These resummations can be compared with the solution to the integro-differential equation corresponding to (10), i.e. (16) where N = ∞ n=0 T n N (n) and N(T = 0) = χ{1, 0}. We have solved (16) numerically with the Euler method and found good agreement with the resummations above. However, it takes much more time to solve (16) because we need T max /dT steps, which is typically more than the ten (or fewer) steps we needed in the resummation approach. At , which gives a condition for c since the right-hand side of (16) is not automatically O(1/T 2 ). As mentioned, we expect c to be independent of χ, and now we can confirm that this is consistent with (16).
In order to obtain N (n) from N (n−1) we need to calculate both components of N (n−1) , even if we at the end are only interested in unpolarized particles. Hence, in calculating the results above we have obtained the necessary information to study polarized initial state as a byproduct. We begin with the leading order for χ 1.
. This is an example of a first-order difference equation which can be solved with a general method, see e.g. Eq. (2.2.7) in [33]. We find where H n is the n-th harmonic number. This compact formula allows us to find an exact resummation, This all-order result has the same radius of convergence as the solution to LL. However, for ∆ P we see that the convergence-limiting singularity is a combination of a pole and a branch point. We also see that ∆ P too becomes independent of χ to leading order in T 1, although this time the next-to-leading order is only logarithmically suppressed. Another difference is that (19) is non-monotonic with a maximum |∆ P | at T ∼ 1/χ.
The full quantum result can be obtained as described above. Although (19) might suggest using a resummation involving logarithms, we still find a fast convergence with a Padé approximants as in (14), except that P − > 0 implies that |∆ P | must be smaller than P ↑ + P ↓ and for that to hold at large T we need J > I. In practice, different choices of I and J can give good approximations, and a "wrong" choice just means that we need to include more terms. Fig. 4 shows that the classical prediction overestimates the peak of |∆ P | by about a factor of two for χ = 0.2, but is close to the quantum result for large T .
In conclusion, we have derived a recursion formula for calculating the expectation value of an electron in a plane wave background field, where O(α n ) is obtained by multiplying the spin dependent O(α n−1 ) result with a Mueller matrix for nonlinear Compton scattering and the corresponding loop and then integrating over the photon momentum. In the classical limit we find the solution to LL to all orders. We have shown that quantum RR can be obtained either by constructing interpolation functions of each order, or by expanding each O(α n ) in χ and then resumming these asymptotic series with e.g. the Borel-Padé method. We have also found that the α expansion obtained in this way can be be resummed with Padé approximants. The fast convergence of these approximants is encouraging for the generalization to nonconstant fields. Our approach takes the spin into account using Mueller matrices, so we have also studied spin dependent RR. In the main text we showed that the classical limit of P − agrees with the solution to LL. Now we will extend this to the transverse components. Since this introduces more variables, we use for convenience a more compact notation. Let p −â C and p −â L be the contribution to P (1) coming from Compton scattering and the loop, respectively, where p = p − l is the momentum of the electron after emitting a photon with momentum l, and a = a(σ) whereσ is some upper limit for σ. In general we would need to use the Mueller matrices for P ⊥ (1) . However, to compare with classical physics we only have to consider M = {1, 0} · M · {1, 0}. We have which corresponds to Eq. (9) in [18]. To obtain higher orders we let the lower integration limit for the σ integral be σ rather than −∞. We can obtain this from [18,19] We could drop the π term since it does not contribute to the leading order for long pulse length or high intensity, and we can only expect the Mueller matrix/incoherent product approach to give the leading order. However, we will keep it since at least in this case doing so allows us to obtain all terms in the solution to LL. To obtain O(α 2 ) we should prepend a Compton scattering or a loop step before the first-order result, where the σ integral is restricted by σ > σ. Since π ⊥ = −a ⊥ is independent of the momentum, the first term, ∆π ⊥ (σ), depends only on the longitudinal momentum and with a linear scaling, and can therefore be treated in exactly the same way as the above calculation for kP . For the second term we have by expanding to linear order in the photon momentum kp (p −a) ⊥ −kp(p−a) ⊥ ≈ −kl(p−a) ⊥ +kp(−l) ⊥ , (A4) so we only need to calculate −l C which we have already obtained for P (1) , we just have to remember the σ ordering. In general we find where a sum over a = −, 1, 2 is implied. For the contribution to P ⊥ (n) from the part without π we find where I = −a 2 and I ⊥ = −a 2 (p − a) ⊥ . Thus, we again find a geometric series and by resumming this we find (7). The classical limit of P µ is on shell, i.e. P + can be obtained from (7). So, the following, direct calculation of P + can be seen as an additional check. We now drop the π term. We will use (A5) (where the sum is still over a = −, 1, 2). From whereππ = 1 + (1/2)(â − a) 2 is independent of p µ , we find ∂ P + (1) /∂p − = 0, so only a = 1, 2 contribute and This is quite different from P + (1) , but at third order we find from which we start to see a pattern. Thus, after P + we find P + (n) ≈ (n − 2)!∆ n 2p − dσ 1 ...

Appendix B: χ expansion
The χ expansion can be obtained by changing integration variables from q = χγ/(1 + χγ) to γ, which removes χ from the argument of the Airy functions, and then simply expanding the integrand. The resulting integrals are These gamma functions are responsible for the factorial growth of the coefficients in the χ expansion. Note that the terms in the χ expansion are completely analytical, but higher orders involve integers with a very large number of digits. Hence, for finding Padé approximants of the Borel transform and for performing the Laplace integral, one may have to work with very high precision.