Quark production in heavy ion collisions I. Formalism and boost invariant fermionic light-cone mode functions

We revisit the problem of quark production in high energy heavy ion collisions, at leading order in $\alpha_s$ in the color glass condensate framework. In this first paper, we setup the formalism and express the quark spectrum in terms of a set of mode functions of the Dirac equation. We determine analytically their initial value in the Fock-Schwinger gauge on a proper time surface $Q_s\tau_0\ll 1$, in a basis which makes manifest the boost invariance properties of this problem. We also describe a statistical algorithm to perform the sampling of the mode functions.


Introduction
In heavy ion collisions at ultrarelativistic energies, such as those performed at the RHIC or the LHC, the bulk of particle production originates from soft gluons that carry a small fraction of the projectile longitudinal momentum [1,2,3,4,5,6]. Because of the infrared singularity present in the emission probability of soft massless gluons, the occupation number of these gluons increases as an inverse power of the longitudinal momentum fraction x, according to the BFKL evolution equation [7,8]. When it reaches values of order 1/α s , non-linear processes such as recombinations become important and tame the growth of the occupation number -a phenomenon known as gluon saturation [1].
By virtue of this large occupation number, the dynamics of these soft gluons is essentially classical, but non-perturbative because highly non-linear [9,10]. The Color Glass Condensate (CGC) effective theory [11,12,13] provides an organization principle for this regime of strong interactions, and a calculational framework for computing observables relevant to hadronic or nuclear collisions involving such densely occupied projectiles.
In the study of heavy ion collisions, the CGC has been applied to calculate the gluon yield at early times [14,15]. At leading order (tree level), this amounts to solving the classical Yang-Mills equations with light-cone currents representing the fast color charges of the two projectiles. At next-to-leading order (one loop) [16], one can extract the terms that contain logarithms of the collision energy and show that they can be absorbed into the renormalization group evolution -according to the JIMWLK equation [17,18] -of the probability distribution of the above currents.
The CGC framework can also be used in order to study the production of quarks in heavy ion collisions. In this framework, the light-cone color currents couple only to gluons (because gluons are the dominant constituents of high energy hadrons or nuclei), and quarks are produced indirectly from the gluons by the process gluons → qq. Thus, the quark spectrum is one order higher in α s than the gluon spectrum. Equivalently, one may say that the quark spectrum is a 1-loop quantity while the gluon spectrum is a tree level quantity. It is well known that the single inclusive quark spectrum can be expressed in terms of a basis of solutions of the Dirac equation, with a color background field that corresponds to the LO gluons (i.e. the classical solution of the Yang-Mills equations). The choice of this basis of solutions is not unique. When they are chosen in such a way that they coincide with the free spinors v s (k)e ik·x or u s (k)e −ik·x in the remote past, they are often called mode functions in the literature, and we will also adopt this terminology in the rest of this paper.
This approach can be used to study the production of fermions or scalar particles in any situation where the production is due to a classical background field [19,20,21,22,23]. In the context of heavy ion collisions, this formulation 1 has been used first in studies of electron production in nuclear collisions [29,30] (although this is a pure QED process, its treatment in the "equivalent photon" approximation is very similar to the CGC), and later in a computation of quark production in heavy ion collisions [31,32]. This earlier work was limited in a number of ways: (i) the basis of mode functions that was used was expressed in terms of proper time and the Cartesian coordinates x, y, z, making the boost invariance of the problem highly non-obvious, (ii) the sum over the modes was restricted to a subset of all the possible modes, and (iii) the resulting quark spectrum may be contaminated by spurious lattice doublers at high momentum.
The goal of this work is to revisit this study in order to overcome all these limitations. In this first paper, we first obtain a new basis for the Dirac mode functions, that naturally depend on the proper time τ , on the rapidity η and on the transverse position x ⊥ . These mode functions are indexed by the transverse momentum k ⊥ and a wave number ν which is the Fourier conjugate to η, making them very convenient for a lattice implementation where the grid covers a fixed range in η. In order to improve the sampling of the mode functions, we use spinors that are random linear superpositions of all the possible mode functions. A proper choice of the distribution of the random weights ensures that the exact result is recovered in the limit of infinite statistics. With finite statistics, this procedure provides a straightforward way to estimate the statistical errors. Numerical results based on a lattice implementation of this framework will be presented in a forthcoming paper.
The contents of the paper is the following : In the section 2, we briefly remind the reader of the Color Glass Condensate and of the expression of the quark spectrum in this framework. We also show in this section how to choose a basis of mode function that makes boost invariance manifest. In the section 3, we present a statistical method to sample the modes, and derive the formula for the corresponding statistical errors. The initial value of the mode functions on the forward light-cone (i.e. just after the collision of the two nuclei) is derived in the section 4. The section 5 is devoted to concluding remarks. A few appendices collect more technical material. The derivation of the expression of the quark inclusive spectrum in terms of Dirac mode functions is recalled in the appendix A, and an alternate derivation following more closely standard Feynman perturbation theory is presented in the appendix B. The appendix C discusses a technicality in the derivation of the initial value of the mode functions, and the appendix D is devoted to the study of a conserved inner product between the mode functions. We make an extensive use of this inner product in order to properly normalize the mode functions, and as a consistency check at various stages of the calculation. In the section E, we use the QED version of the mode functions derived in the section 4 in order to recover the electron production amplitude in the collision of two electrical charges.
2 Quark yield in the CGC framework 2

.1 Color Glass Condensate
The Color Glass Condensate framework is an effective theory that can be used to study the early stages of heavy ion collisions, summarized by the following Lagrangian density, written here for one family 2 of quarks of mass m. The color charge content of the incoming nuclei is described by the two currents J µ 1,2 , whose supports are restricted to the light-cones, J µ 1 ∝ δ(x − ) and J µ 2 ∝ δ(x + ), in a collision at very high energy. These currents fluctuate event-by-event, with a Gaussian probability distribution in the McLerran-Venugopalan (MV) model 3 [9, 10] that we use in this paper. If one is interested in the production of quarks in a given collision, one could draw randomly one configuration of J µ 1,2 , and not perform an average over these currents 4 . In the gluon saturation regime, the currents J 1,2 are inversely proportional to the gauge coupling, For this reason, gluonic observables at leading order are expressible in terms of a classical color field that obeys the Yang-Mills equation with the source J 1 +J 2 , One should in principle also impose the covariant conservation of the current [D µ , J µ 1 +J µ 2 ] = 0. This constraint becomes trivial in the Fock-Schwinger gauge, x + A − + x − A + = 0, since it ensures that the gauge potential vanishes on the support of the current. We adopt this gauge in the following. Furthermore, for inclusive observables, one can prove that this equation of motion must be supplemented by a retarded boundary condition [33], such that the gauge field vanishes in the remote past, thereby making this classical solution unique. In the saturation regime, this classical gauge field is strong, of order A µ ∼ 1/g.

Inclusive quark spectrum
In the CGC framework, the fermions do not couple directly to the currents J µ 1,2 , but only indirectly through the gauge field that appears in the covariant derivative in the Dirac operator i/ D − m. Therefore, the natural order of magnitude of the spinors is ψ ∼ 1, in accordance with the fact that the occupation number of fermions is bounded by unity. Thus, observables that contain quark fields are of higher order in the gauge coupling. For instance, the g 2 power counting for the quark spectrum at LO is the same as that of the gluon spectrum at NLO: both are 1-loop quantities, the only difference being the nature of the field running in the loop (quark versus gluon). In a fixed background color field, the quark spectrum at LO is given by the following formula 5 , whose derivation is recalled in the appendix A : 4 Note however that the theoretical basis for doing this is not very robust, and becomes inconsistent beyond leading order. Indeed, for fixed J µ 1,2 , loop corrections to observables contain unphysical logarithms of the longitudinal momentum cutoff that separate the gluon modes that are described by the sources J µ and those that are described as gauge fields. These logarithms can be absorbed into a redefinition of the probability distribution W [J] of these sources (that must now evolve according to the JIMWLK equation). But this procedure only works if one performs an average over the sources. 5 This formula is true to all orders in the currents J µ 1 and J µ 2 . If one expands it to lowest order in these currents (i.e. dNq/d 3 p ∝ O((J 1 J 2 ) 2 )), one recovers the standard result for the process gg → qq with off-shell incoming gluons, derived in the framework of k T -factorization in refs. [34,35] (see ref. [24] for this comparison).
pσb is a free positive energy spinor of momentum p, spin σ and color b (since quarks live in the fundamental representation of the gauge group SU (N c ), this color index runs from 1 to N c ). In the absence of background field, these spinors are given by 6 , However, it may also happen that the gauge fields at x 0 → +∞ evolve into a nonzero pure gauge configuration. In this case, the above spinor should be replaced by a color rotated one: where SU (N c ) is the SU (N c ) matrix defining the pure gauge background. In contrast, ψ − ksa is a spinor that has evolved over the background color field, starting at x 0 = −∞ from a negative energy free spinor of momentum k and spin s : Note that the subscripts a, b refer to the initial color of the quarks. The color they carry at the point x is not written explicitly, and is encoded in the N c (color) × 4 (Dirac) components of the spinors. The inner product · · x 0 that appears under the integral in eq. (4) is defined by (In the product ψ † χ, all the unwritten color and Dirac indices are contracted.) The properties of this inner product are studied in detail in the appendix D.
In this appendix, we also use its conservation as a consistency check of the results that will be derived in the section 4. Note that the formula for the quark spectrum requires that one takes the limit of infinite time. As we shall discuss later in this section, this is also a requirement for the quark spectrum to be gauge invariant. Eq. (4) is the expression for the fully inclusive spectrum that we are going to use in the rest of this paper. The virtue of this formula is that it reduces the calculation of a one-loop graph in a background field to solving a (linear) partial differential equation with retarded boundary conditions. Even if this can be done analytically only for very simple backgrounds, this problem can in principle be tackled numerically for completely general backgrounds. 6 In this equation and in the rest of this paper, we write explicitly only the indices that characterize the initial value of the spinor. A more complete notation would read : where α is the Dirac index and b is the color of the quark at the point x.
Note that in eq. (4), the spectrum is summed over all the possible final states and over the spin of the tagged quark. In order to obtain the polarized spectrum, one simply needs to remove the sum over the spin index σ. This formula also contains sums over the colors of the initial and final fermion. These sums should not be undone, as the spectrum of quarks with a given color has no gauge invariant meaning. k and s can be viewed as the momentum and spin of the antiquark that must be produced along with the quark to satisfy the conservation of the flavor quantum number, as reflected by the initial condition for the spinor ψ ks in the remote past.

Change of coordinates
Since collisions in the high energy limit are invariant under boosts along the longitudinal axis 7 , it is convenient to trade the longitudinal components of the momenta p z , k z in favor of the corresponding rapidities y p and y k . Likewise, the proper time τ and spatial rapidity η are more suitable than x 0 and z to map the space-time: Besides the obvious change in the measure dp z /ω p = dy, one must alter the definition of the inner product so that the integration is on a surface of constant τ (instead of a constant x 0 ), (See the appendix D.) When doing this, eq. (4) becomes The boost invariance of the problem implies that the inner product depends only on the difference of the rapidities y p − y k . After integration over y k , the resulting quark spectrum is independent of the rapidity y p .

Boost invariant spinors
The boost invariance can be made manifest at the level of the spinors ψ 0+ p ⊥ ypσb and ψ − k ⊥ y k sa themselves. Even when the background field is invariant under boosts in the z direction, these spinors depend separately on the momentum 7 Here, we are disregarding the small-x evolution of the color sources in the incoming nuclei. This effect would break the boost invariance of the problem due to gluon loop corrections, and make the quark spectrum depend on rapidity on scales ∆y ∼ α −1 s . rapidity y and on the spacetime rapidity η. This can be trivially seen on the vacuum spinors, whose rapidity dependence can be made explicit as follows where M p ≡ p 2 ⊥ + m 2 is the transverse mass. To turn the prefactor into a function of y p − η, it is convenient to define transformed spinors as follows, The factor √ τ has been introduced for later convenience. After this transformation, the new spinors ψ k ⊥ ysa are boost invariant, in the sense that they depend on the spatial rapidity η and on the momentum rapidity y only through the difference y − η (provided that the background field does not depend on η).
The boosted spinors introduced in eq. (13) also offer the advantage of obeying a simpler form of the Dirac equation where rapidity does not appear explicitly in the coefficients. In order to see this, first note that Then, multiply this operator on the left by exp(− η 2 γ 0 γ 3 ). A simple calculation gives From this observation, we conclude that the modified spinors ψ obey the following Dirac equation : One can see that the coefficients of this equation are independent of the rapidity η when the background field is boost invariant (so that there is no η dependence hidden in the covariant derivatives). In terms of the boost invariant spinors defined in eq. (13), the inner product on a constant proper time surface takes a particularly simple form,

Mode functions in the ν basis
Another useful transformation is to go from a basis where the spinors have a definite momentum rapidity y to a basis where they have a fixed Fourier conjugate ν to the space-time rapidity η, When the background field is boost invariant (i.e. independent of η), ψ k ⊥ ysa depends only on y − η and the spinor ψ k ⊥ νsa in the new basis has a trivial η dependence in exp(iνη): ν is a conserved quantum number and the η dependence of these spinors is not altered by their propagation over the background field. Moreover, the Dirac equation obeyed by these spinors is effectively 2 + 1 dimensional, since the η dependence can be factored out. When calculating the inner product of two such spinors (see eq. (141)), the integration over η trivially yields a delta function, where we denote by · · τ the "reduced" inner product that remains after one has factored out the delta function. In this basis, the quark spectrum is given by In this formula, it is tempting to ignore the limit τ → +∞ and to interpret the resulting expression as the quark spectrum at the finite proper time τ . One should however consider such a generalization with caution, since it is not possible to rigorously define asymptotic states at a finite time.

Gauge invariance
Under a local gauge transformation, a spinor ψ(x) is transformed as follows where Ω(x) is an SU (N c ) matrix. Since the inner product that enters in the quark spectrum given by eqs. (4) or (22) is local, it is gauge invariant provided of course that the spinors ψ − and ψ 0+ † are gauge rotated consistently. In eq. (6), we have already indicated that the spinors ψ 0+ † should be obtained from the free solutions in a null background (5) by an appropriate color rotation, if the background field at the time of quark measurement is a nonzero pure gauge. The square of the inner product that appears in eq. (4) can be written as For the time being, let assume that the background field is a pure gauge at the time x 0 where the quarks are being measured 8 . The factor Ω(y)Ω † (x) depends on this pure gauge, and can be obtained by a Wilson line between the points x and y, where γ is a path from x to y in the hyperplane of fixed time x 0 . Thus, in practice one would calculate When the background field is a pure gauge, this does not depend on the path γ chosen between x and y. If one extends the use of eq. (26) to a situation where the background field at the time x 0 is not a pure gauge, one still obtains a gauge invariant result, but there is now an ambiguity due to the choice of the path. Indeed, Wilson lines U γ and U γ evaluated on two different paths differ by a Wilson loop, defined over the closed loop γ ∪ γ −1 made of the path γ followed by the reverse of the path γ . This Wilson loop measures the chromo-magnetic flux across the closed loop, and is therefore equal to the identity only if the background field is a pure gauge. One expects this irreducible ambiguity to decrease with the quark momentum. On the one hand, in the spectrum of quarks of momentum p, the typical spatial separation x − y is of order 1/|p| (since x − y and p are Fourier conjugates in eq. (24)), and the closed loops that one would have to consider in the above argument have a typical area of order 1/|p| 2 . On the other hand, numerical studies of the gauge fields produced in the McLerran-Venugopalan model indicate that the expectation value of Wilson loops decreases exponentially with the area of the loop when it becomes larger than Q −2 s [36,37]. This suggests that this ambiguity should not affect much the quark spectrum for |p| Q s .

Sketch of a direct algorithm
Eq. (22) contains the essence of our procedure for calculating the quark spectrum: i. Draw randomly a pair of sources J µ 1,2 , and solve the classical Yang-Mills equation (3) with null retarded initial conditions.
ii. For a given transverse momentum k ⊥ , wavenumber ν, spin s and color a, initialize the spinor ψ k ⊥ νsa as a free negative energy spinor.
iii. Solve the reduced 2+1 dimensional Dirac equation (20) for the time evolution of this spinor over the color field found in the step i.
iv. At some sufficiently large time, project the spinor ψ k ⊥ νsa on a free positive energy spinor of transverse momentum p ⊥ , same wavenumber ν, spin σ and color b. This gives the reduced inner product ψ 0+ Repeat the steps ii to iv in order to sum over all k ⊥ , ν, s, a's.
vi. (optional) Repeat steps i to v in order to average over the configurations of the color sources J µ 1,2 .
If we discretize the transverse plane into a N x × N y grid and the rapidity axis into a grid with N η points, the computational cost for solving the Dirac equation for one mode function scales as N x N y N τ where N τ is the number of timesteps. Note that this cost does not depend on N η since the η dependence can be factored out for a boost invariant background. This must be repeated for the N x N y N η mode functions. Therefore, the total computational cost scales as (N x N y ) 2 N η N τ .

Statistical method
The algorithm described in the previous subsection is deterministic, but it suffers from an unfavorable scaling with the size of the transverse grid, since the computational cost scales as (N x N y ) 2 . Instead of doing in full the sum over all the modes, it is possible to do it by a Monte-Carlo sampling in which one uses random linear superpositions of the mode functions where the coefficients c k ⊥ νsa are random numbers with the following variance The justification of this approach is to note that the projector on the subspace of negative energy spinors can be rewritten as follows where P[c] is a normalized probability distribution (in practice a Gaussian distribution) with zero mean value and a variance given in eq. (29). In eq. (22), the integrals over k ⊥ , ν and the sums over s, a are then replaced by a statistical average over these random numbers, (31) (L η is the total length of the η interval represented in the lattice implementation.) Since the random sum in eq. (28) mixes 9 the various ν's, the evolution of ψ − c is governed by the 3+1 dimensional Dirac equation (16), and the computational cost of this approach scales as N x N y N η N τ N conf where N conf is the number of samples used in the statistical average. Compared to the direct deterministic method, a power of N x N y has been replaced by N conf , which is advantageous for large grids if N conf N x N y .

Statistical errors
The statistical method summarized by the eqs. (30) and (31) is exact only in the case of a perfect sampling of the Gaussian distribution P[c]. In practice, this sampling is performed by generating a large but finite number N conf of configurations. In doing this, the left hand side of eq. (30) is replaced by where we have used discrete notations that correspond to the lattice implementation, and we use the following shorthands : (The integers j x,y,η label the momentum modes k ⊥ , ν in the lattice implementation, and s, a are the spin and color quantum numbers.) The coefficients that appear in the sum in eq. (32) result from N conf samplings of the Gaussian distribution : In this equation, c (n) J is the n-th random sample for the mode J . With N conf samples, the observables we are interested in (e.g. the inclusive quark spectrum given by eq. (31)) are generically approximated by where F denotes the quantum numbers of the final state, f a partial sum over these quantum numbers (in the case of eq. (31), this partial sum is over the wave number ν, the spin and color of the produced quark), and N a normalization factor. C N conf ( J , J ) is itself a random number, whose distribution can be determined in the large N conf limit by a method similar to the derivation of the central limit theorem. The mean value and variance of this approximation can be obtained from those of C N conf ( J , J ). Let us first recall that This leads easily to The formula for the variance is exact if the distribution of c J is Gaussian. Like in the central limit theorem, the variance of C N conf decreases as 1/N conf . From the first of eqs. (37), we obtain the mean value of O N conf which is indeed the exact value of the observable. Using the second of eqs. (37), we get We see that the standard deviation of O N conf decreases as 1/ √ N conf , with a coefficient that has a non-trivial covariance. It can itself be estimated by the statistical method as follows : i. Define two random linear superpositions of the negative energy mode functions : with uncorrelated random weights c iii. Compute the following quantity : iv. Repeat the steps i-iii in order to average over the random numbers c Since this is just an error estimate, a small number of samples is sufficient. In practice, one may divide the N conf samples already calculated in two subsets, and use these subsets to evaluate the error.
The standard deviation of O N conf is the square root of the result of this computation, divided by √ N conf . Note that the summand in eq. (41) is a complex number, which can lead to phase cancellations when summing over the final quantum numbers f . These cancellations are more effective for more inclusive observables, thanks to a more extended sum on f .

Relation to "low-cost fermions"
In real-time lattice simulations of fermions, the so-called low-cost fermion method [38] has been used in several works, e.g. [39,40,41,42,43,44]. Let us briefly compare this method with our approach. In our statistical method, we compute the following quantity using the stochastic field (28): where J comprises all quantum numbers including momentum, and O is a matrix that depends on the observable we wish to compute. In the case of the spectrum (31), In the low-cost fermion method, instead of using one stochastic field (28), one employs two kinds of stochastic fields called "male" and "female" fields: where c J and d J are independent random numbers that have the same variance as eq. (36). Combining these two fields, one can compute (1l is the unit matrix in the spin, color and position of the spinors at the time t), we can relate the quantities evaluated in our method (42) and in the low-cost fermion method (44) by (for 2 spin states and 2 colors.) Therefore, the two methods provide the same result 10 up to a trivial additive term, that can be interpreted as a constant vacuum contribution. However, our method has two advantages over the low-cost fermion method. Firstly, it is numerically less costly than the low-cost fermion method, simply because it uses only one kind of stochastic field. Secondly, the statistical errors are smaller for the evaluation of the spectrum. In our method (42), the spectrum is directly obtained from the statistical ensemble, without any subtraction. On the other hand, in the low-cost fermion method (44), one gets directly access to 1 2 − f (f being the fermion occupation number), and the vacuum 1/2 must be subtracted later. Because this vacuum 1/2 also contains statistical errors due to the Monte-Carlo sampling, the low-cost fermion method suffers from comparatively larger statistical errors, especially when the value of the occupation number is small compared to 1/2.

Fermionic mode functions on the light-cone 4.1 Background gauge field
In order to use the classical-statistical method in the calculation of the quark spectrum, we should solve the Dirac equation with a background color field, starting with a free spinor at x 0 = −∞. However, it is not straightforward to do this numerically, due to the singular nature of the gauge field on the lightcones x ± = 0. The field strength on these lines is proportional to a δ(x ± ), which cannot be handled easily in a numerical program. Similarly to the case of gluon production, one should first solve the Dirac equation analytically up to a surface τ = τ 0 Q −1 s , and perform the numerical resolution only in the forward light-cone for τ ≥ τ 0 .
In this section, we generalize to the case of spinors the derivation that was performed in ref. [45]. Since the Dirac equation is linear, its solution can be written as the sum of a left-moving and a right-moving partial waves, as illustrated in the left diagram in the figure 1. These two partial waves are totally independent. We take advantage of this fact to choose the gauge for the background field differently for each of them, in order to simplify the resolution of the Dirac equation. Of course, before adding up the two partial waves in order to construct the full solution, we must perform a gauge rotation that brings them to a common gauge.
x + x - Figure 1: Left : decomposition of the solution of the Dirac equation into leftmoving and right-moving partial waves. Right : structure of the gauge field produced by the two nuclei in the A − = 0 gauge, that we use for solving the Dirac equation for the right-moving partial wave. x ± denote light-cone coordinates, x ± ≡ (t ± z)/ √ 2.
For the right-moving partial wave, it is convenient to work in the light-cone gauge A − = 0, in the same spirit as what was done in ref. [45] for the gluons. After having determined the right-moving spinor at the proper time τ 0 inside the forward light-cone, we will rotate it to the Fock-Schwinger gauge. The calculation for the left-moving spinor (to be performed in the light-cone gauge A + = 0) will not be performed here, since the result can be guessed from the other partial wave by symmetry. The structure of the background field in the A − = 0 is shown in the right diagram of the figure 1. It is made of the following two elements : 1. the nucleus moving in the +z direction produces a field A + 1 , proportional to a δ(x − ) and independent of x + .
2. in the region x + > 0, x − < 0, the nucleus moving in the −z direction produces a field A i 1 , which has the form of a transverse pure gauge, 3. in the strip corresponding to the shock-wave of the nucleus moving in the +z direction, this field A i 2 receives a color precession induced by the first nucleus. This color rotated field reads where Note that eq. (48) is equivalent to Therefore, by starting the evolution at x 0 = −∞, the right-moving partial wave encounters first the field A i 2 and then the field A + 1 . In order to express the quark spectrum, we need the negative energy spinors, ψ − ksa (k is the 3-momentum of the incoming quark, s its spin and a its color). Before they encounter the background color field, they read simply Since the background color field has only a finite jump on the half-line x + = 0, x − < 0, the spinors are continuous across this line, and the above formulas remain valid up to x + = 0 + (just above this line).

Evolution in the region
The next step is to solve the Dirac equation in the region x + > 0, x − < 0. Since the background field is a pure gauge in this region, the covariant derivative can be written as Therefore, the new spinor defined byψ − ksa ≡ U 2 ψ − ksa obeys the free Dirac equation : The solution of this equation can be expressed in terms of the initial value of ψ − ksa on the surface x + = 0 + by the following Green's formula, where S 0 R (x, y) is the bare retarded propagator for a quark, Note that, even if we have not restricted the integration range for the variable y − in eq. (54), the fact that the support of the retarded propagator S 0 R is limited to x 0 > y 0 , (x − y) 2 ≥ 0 imposes that y − ≤ x − . Therefore, for a point x located in the region x + > 0, x − < 0 (shaded in green in the figure 1), only points with y − < 0 on the initial surface can contribute. This is the reason we can solve independently the equation for the left-and right-moving partial waves.
By introducing the Fourier representation (55) of the retarded propagator in the Green's formula (54), one can perform most of the integrations analytically except the final integration over transverse momentum. This leads to where M p ≡ p 2 ⊥ + m 2 is the transverse mass and where we have introduced the projector P + ≡ γ − γ + /2. We have also introduced the Fourier transform of the Wilson line U 2 One can check that eq. (56) falls back to the free spinor v s (k) exp(ik · x) if the background field is turned off (U 2 = 1).

Evolution across the line
The next step is to propagate the spinor (56) across the field A + 1 of the nucleus moving in the +z direction. The region of space-time supporting this field is infinitesimal (since A + 1 ∼ δ(x − )), but the infinite strength of the gauge field in this region nevertheless produces a finite change of the spinors. In this region, there is also an O(1) transverse component of the gauge field, given in eq. (48). The Dirac equation, can be separated into a part that depends on the background field A + 1 and a constraint independent of A + 1 by multiplying 11 it by the projectors P + or P − ≡ γ + γ − /2 : The first equation, independent of the background field, is a constraint that relates the two projections of the spinors at every x − (note that this equation does not contain the ∂ + derivative). The second equation determines the dynamical evolution (in the variable x − ) of the P − projection under the influence of the background field A + 1 . Inserting the first equation into the second gives a second order equation that drives the evolution of the P − projection, In this equation, the ∂ + derivative and the field A + 1 are large (inversely proportional to the thickness of the shock-wave that supports the A + 1 field), while all the other terms do not have this large factor. Physically, keeping only the term in ∂ + − igA + 1 leads to the eikonal approximation, where the fermion would propagate on a straight line along the x − axis, while the terms −D 2 ⊥ + m 2 lead to some transverse diffusion with respect to this axis. In the limit where the thickness of the shock-wave goes to zero, we can neglect it : (Note that this approximation is only valid to cross the shock-wave, and should not be used to evolve at a finite distance from the shock-wave.) The solution of this equation is very simple, where the x − -dependent Wilson line U 1 was defined in eq. (49). The spinor at x − = 0 that appears in the right hand side is given by eq. (56). Then, the constraint equation can be solved to give Note that the constraint defines the P + projection of the spinor only up to an arbitrary function of x − and x ⊥ . This "integration constant" can be determined by requesting that we recover the P + projection of a free spinor when all the Wilson lines are set to the identity.

Transformation to Fock-Schwinger gauge
In order to gauge transform these spinors to the Fock-Schwinger gauge, it is sufficient 12 to multiply them by U † 1 , When this transformation is applied to eqs. (62) and (63), the Wilson line U 1 appears in two types of combinations : 12 In the case of the background gluon field, an additional transformation was necessary, see the eq. (18) of ref. [45]. But since the color rotation Ω that characterizes this transformation is equal to the identity for proper times τ Q −1 s (see the eq. (19) in ref. [45]), it has no effect on the spinors.
The second of these structures can be simplified if we recall that D i is the covariant derivative built with the field of eq. (50). Therefore, where A i 1 is defined in the same way as A i 2 (see the eq. (47)), Note that the field A i 1 + A i 2 that appears in this equation is nothing but the transverse component of the gauge potential in the Fock-Schwinger gauge at τ = 0 + , hence the notation D i FS for the resulting covariant derivative. Therefore, the two projections of the right-moving negative energy spinors in the Fock-Schwinger gauge read (at x − = 0 + , just above the shock-wave) The eqs. (68) for the right-moving spinors must be completed by a set of similar equations for the left-moving spinors. These can be obtained from the above formulas by the following substitutions At this point, we have the components of the negative energy mode functions on the light-cone τ = 0 + (i.e. just after the collision), which provides all the necessary initial data for studying their evolution after the collision. As explained in the appendix C, these formulas can also be used as initial conditions at some initial time τ 0 > 0, provided that τ 0 a ⊥ where a ⊥ is the transverse lattice spacing used in the numerical resolution.

Mode functions in the ν basis
So far, we have derived the mode functions ψ ksa in terms of the Cartesian 3momentum k. However, as explained in the section 2.3,the boost invariance of a high energy collision is more manifest if one uses the modified spinors defined in eq. (13) and if one further goes to a basis where the spinors are labeled by the quantum number ν (Fourier conjugate to η) instead of y. It is easy to obtain these new mode functions by the following transformation : For the right-moving partial waves, the integral that enters in the transformation of P ψ − is of the form while for the left moving partial waves, one needs These integrals can be expressed in terms of the Γ function, Let us now recapitulate our results for the fermionic mode functions on the light-cone after this transformation, after summing the right-moving and leftmoving partial waves, and including both the P + and P − projections In this formula, we have kept only the terms that are non-vanishing in the limit τ → 0 + . Therefore, its use should be restricted to very early times.

Summary and outlook
In this paper, we have reconsidered the problem of quark production in heavy ion collisions in the color glass condensate framework. Our approach is closely following that of refs. [31,32], which expresses the leading order inclusive quark spectrum in terms of a set of mode functions of the Dirac equation, but overcomes a number of limitations of this earlier work.
Firstly, we have defined a basis of fermionic mode functions that are more appropriate for the boost invariant expanding geometry of a high energy collision. In particular, since these mode functions are indexed by the Fourier conjugate ν to rapidity, they are especially suitable for a lattice implementation in which one discretizes the rapidity axis. We have calculated analytically the value of these mode functions just after the collision, at a proper time Q s τ 1 and in the Fock-Schwinger gauge, in terms of the Wilson lines that represent the classical color background field of the two colliding nuclei. Thanks to these analytical initial values, one will not have to deal with crossing the light-cones in the numerical resolution of the Dirac equation.
Secondly, we have exposed a statistical method for sampling the set of modes over which one must sum in the calculation of the quark spectrum. This method ensures that no mode is left out, while considerably reducing the computing time compared to a complete sum over all the modes. This approach also provides a more robust way of estimating the error one makes in the sum over the modes functions.
In a forthcoming paper, we will apply the formalism that we have setup in the present paper to a study of quark production in two situations. Firstly, we will present a test of the method in the case of a background field for which one can solve analytically the Dirac equation for the mode functions (a constant SU (2) chromo-electrical field). Then, we will present results on quark production in heavy ion collisions, in the case where the background color field is given by the MV model. In order to mitigate the problems caused by the lattice fermion doubler modes, we will also explain how our framework must be modified in order to include a Wilson term in the fermionic action.

A Quark spectrum in the Schwinger-Keldysh formalism
In this appendix, we recall the derivation of the formulas (4) for the inclusive quark spectrum, starting from the standard LSZ reduction formulas. For simplicity, we consider only one flavor of quarks. Since in the CGC framework, the external sources are only coupled to gluons and quarks can only produced in quark-antiquark pairs (the net flavor number is always zero).

A.1 Single quark pair amplitude
Firstly, let us consider the amplitude for producing a single quark-antiquark pair, In order to keep the notations concise, we are not writing explicitly the color and spin indices of the quark and antiquark. The operator b † out (p) (resp. d † out (q)) creates a quark of momentum p (resp. an antiquark of momentum q).
Using the decomposition of the free field operator ψ out (x) as a superposition of free modes, we have (In these formulas, p 0 = p 2 + m 2 and q 0 = q 2 + m 2 .) Note that the time t at which these formulas are evaluated do not change the result. Using these formulas, standard manipulations lead to the LSZ reduction formulas for the single quark pair production amplitude, Note that the 2-point correlation function that appears in this formula is a Feynman (i.e. time-ordered) propagator. Its evaluation to all orders in the background gluon field cannot be performed in practice. Indeed, even though it obeys the Dirac equation, its determination is made extremely complicated by the fact that it must satisfy mixed boundary conditions, both at x 0 = −∞ and at x 0 = +∞.

A.2 Inclusive quark spectrum
There is no practical way to calculate the single quark amplitude considered in the previous subsection, in the presence of a strong background gluon field, as is the case in the high energy heavy ion collisions. This is not a big limitation however, because this quantity is also not very phenomenologically useful in such a context. Indeed, for light quarks (quark flavors for which m Q s ), quark production is not a rare process and more than one quark pair are produced in a collision. Therefore, the probability P 1 of producing exactly one quark pair (given by the square of M 1 ) is not very interesting. Much more useful would be the complete probability distribution, P 1 , P 2 , P 3 , etc. Unfortunately, calculating them is as complicated as calculating P 1 . But it turns out that the moments of the probability distribution are much easier to compute. Besides the trivial one ( ∞ n=0 P n = 1), the simplest of these moments is the first moment, n nP n , that gives the mean number of produced pairs. A little more information can be gathered by considering the same quantity in differential form, which is precisely the quark spectrum. In terms of transition amplitudes, this quantity reads where we have used the shorthand for the invariant phase-space of final state quarks and antiquarks. In this formula, the differential probability for producing n + 1 quark-antiquark pairs is weighted by the number of quarks (n + 1), and then integrated over the phasespace of all the antiquarks and of n of the quarks. This quantity is normalized in such a way that it gives the mean number of produced quarks after integration over d 3 p, Most of eq. (78) is in fact the projector on the subspace of states with net flavor number −1, dΦ pi dΦ qi × [p 1 · · · p n ] Q [qq 1 · · · q n ] Q out [p 1 · · · p n ] Q [qq 1 · · · q n ] Q out , (81) that has a trivial action on the state b out (p) 0 in (since this state has also flavor number −1), and eq. (78) can then be reduced to the much more compact form The interpretation of this formula is quite transparent, since it amounts to evaluating the expectation value of the final quark number operator, for a system prepared in the pure state 0 in . Using eqs. (76), this can be rewritten in terms of the quark field operator as follows, The main differences between this expression and the similar expression for the amplitude M 1 are the following : 1. The vacuum state is the in-vacuum state on both sides.

2.
The two spinors are not time ordered.
These differences in fact lead to considerable simplifications in the evaluation of this quantity in the presence of a strong background gluon field. The first step is to note that the 2-point correlator that appears in the integrand of eq. (83) is the component S −+ (x, y) of the fermion 2-point function in the Schwinger-Keldysh formalism. In the presence of a background field A µ , its tree level expression (to all orders in the background field) can be obtained by noticing that it obeys the following equations Next, one should recall the expression of the vacuum propagator S vacuum −+ (x, y), It is then easy to construct a semi-explicit expression for the dressed propagator S −+ (x, y) in terms of a basis of solutions of the Dirac equation: When this expression is inserted into eq. (83), one must evaluate the following expression It is useful to note that Adding this identity to the integrand of eq. (87), we obtain which can be rewritten as a 3-dimensional integral since it is the integral of a total derivative. The boundary at spatial infinity can be dropped if there are no background fields there. The boundary at x 0 → −∞ does not contribute because it leads to the vanishing overlap u † (p)v(p). The only remaining contribution comes from x 0 → ∞, (90) This formula leads immediately to eq. (4).

B Quark spectrum from Feynman amplitudes
In the previous appendix, we presented a derivation of the quark spectrum based on fairly standard many-body manipulations: we first related this spectrum to the expectation value of the quark number operator, and then we evaluated the latter using the Schwinger-Keldysh formalism. However, a more elementary derivation is possible, where the many-body aspects of the problem are treated "by hand". In the present appendix, we present such an alternate derivation, starting from the Feynman amplitudes for producing 1,2,3,... quark-antiquark pairs, and combining them in the appropriate way to obtain the single quark spectrum. This method is a bit more involved since it requires to account for all the final state particle permutations, but it has the advantage of making more tangible the combinatorics that happens under the hood in the derivation of the appendix A.

B.1 Pair production amplitudes
The starting point is the amplitude M 1 (p, q) for producing one quark-antiquark pair, already introduced in eq. (75). This amplitude is made of a time ordered 2-point function connecting the quark of momentum p and the antiquark of momentum q, times a disconnected sum of vacuum-vacuum graphs. The latter is crucial in the presence of a background field, since the sum of the vacuumvacuum graphs is not a pure phase (unlike when the background is the vacuum). In practice, we do not need to calculate this factor, since it is the same in all amplitudes and can therefore be determined at the very end by the request that the sum of all probabilities to produce 0,1,2,3,... quarks is equal to one. For now, we will simply write where M c 1 is the connected part of the pair production amplitude (only this factor carries a dependence on the momenta of the produced quark and antiquark).
In the rest of this appendix, we limit the discussion to the lowest order for the factor M c 1 . At this order, it is simply made of a Feynman propagator connecting the produced quark and antiquark, dressed by the background field: For the sake of simplicity, we will represent this dressed propagator as follows: A more explicit expression of the connected part of the single pair production amplitude is where T F is the dressed Feynman propagator, amputated of its final free propagators. In terms of the dressed (G F ) and free (G 0 F ) Feynman propagators, it is given by The convention for the momenta in eq. (94) is that the 4-momentum −q enters on one side of the propagator and the 4-momentum p exits on the other side. At the same level of approximation, the amplitude for producing n pairs is obtained as follows: M n (p 1 · · · p n , q 1 · · · q n ) = 0 out 0 in where S n is the symmetry group of the set [1, n]. The sum over all the permutations σ ∈ S n is necessary in order to account of all the possible ways to connect the quarks with antiquarks. (σ) is the signature of the permutation σ, resulting from the signs collected when permuting fermion fields.

B.2 Final state combinatorics
It is possible to encapsulate all the information about the distribution of the produced quarks in the following generating functional, where we have used the following shorthand for 1-particle phase-space integrals: The + superscript on the integration symbol indicates that we keep only the positive on-shell energy. Likewise, a − superscript will indicate that the negative on-shell energy is retained: From its definition, it is easy to see that the single quark spectrum is obtained as Note also that unitary (the sum that all probabilities should be one) implies that F[z ≡ 1] = 1. Inserting eq. (96) into eq. (97), we obtain When all the spin and Dirac indices are summed over, the product in the second line forms closed quark loops, from which the spinors can be eliminated by using u(p)u(p) = / p + m and v(q)v(q) = / q − m. It is convenient to change q → −q for all the antiquarks, so that we have: where we have defined 13 In the diagrammatic representation used for this quantity, the dotted line represents the final state. Right of this line is an amplitude and left of this line is a complex conjugated amplitude. In order to go from eq. (101) to eq. (102), we have used the symmetry of the n-quark and n-antiquark phase-spaces and the permutation ρ is defined as ρ ≡ σ −1 σ . Every permutation ρ ∈ S n has a unique decomposition in a product of cycles 14 . In eq. (102), each of these cycles will produce a closed quark loop, that depend on z through the quantity L[z] (linear in z) defined in eq. (103). The degree in z of such a quark loop is the order of the cycle (the number of iterations before the cycle returns to the starting point). For a cycle of order r, we will denote the value of the corresponding quark loop tr ((L[z]) r ), where the trace symbol compactly encapsulates the integrals over all the momenta along the loop, as well as the contractions over Dirac and color indices. The important point is that in eq. (102) the product of L's depends only on the orders of the cycles into which the permutation ρ can be decomposed: if ρ = c 1 c 2 · · · c l where the c j 's are cycles of orders r 1 , r 2 · · · r l , then we have In order to perform the sum over the permutations ρ, it is sufficient to know the number of ρ's that admit a decomposition into a 1 cycles of order 1, a 2 cycles of order 2, ..., a n cycles of order n (with the constraint a 1 + 2a 2 + · · · + na n = n), n! a 1 ! · · · a n ! 1 1 a1 · · · n an , and its signature (ρ) = (−1) n n j=1 (−1) aj .
Combining these results, we can rewrite 15 The second and third lines are equivalent because the constraint j≥1 ja j = n prevents a j 's with j > n from being nonzero. Going from the third to the fourth line merely corresponds to a different way of slicing the sum over all a j 's. The only missing ingredient is the prefactor 0 out 0 in 2 , that can be trivially determined in order to satisfy unitarity. The final expression for the generating functional is therefore Taking a functional derivative leads to the following expression for the quark spectrum: At this stage, we have a representation of the quark spectrum in terms of the object L[z], which is itself quadratic in the Feynman propagator of a quark in a background field. However, this expression contains terms of arbitrary order in this propagator, because of the prefactor (1 − L[1]) −1 . Note that in a weak background field, the quantity L [1] is small and the first term, quadratic in the Feynman propagator, dominates the sum. This is the limit where the pair production probability is small and where at most one pair is produced in a collision. Therefore, the quark spectrum is almost equal to the differential probability of producing a single quark, given by the first diagram in the above series. The second, third, etc... diagrams encode Fermi-Dirac correlations that are only important when more than one quark are likely to be produced (which is the case in heavy ion collisions for quarks whose mass is comparable to the gluon saturation momentum or smaller).

B.3 Expression in terms of retarded amplitudes
It turns out that a much simpler expression, quadratic in the propagator, can be obtained if we rewrite eq. (109) in terms of the retarded propagator of the quark. One can define a retarded analogue T R of T F , built with free retarded propagators instead of free Feynman propagators. The free Feynman and retarded propagators are related by If we denote by V one insertion of the background field, the following two equations define T F and T R recursively from which one first obtains: In terms of these compact notations, L [1] can be written as where ρ + (p) ≡ 2π(/ p + m)θ(+p 0 )δ(p 2 − m 2 ), and manipulations similar to above lead to an expression that depends only on the retarded T R : Note that since L[z] is a linear functional of z(p), the derivative δL[z]/δz(p) is similar to L[1], but with the intermediate momentum p fixed instead of being integrated over. We will denote this as follows: In the same fashion, one also obtains the following identities: Therefore, The outcome of these algebraic manipulations can be pictorially summarized as where the blue lines in the right hand side are retarded propagators. The connection with the previous appendix and the formula (4) resides in the following identities: The eqs. (68) (and their counterparts for the left-moving partial waves) give the value of the fermionic mode functions immediately after the collision, one the semi-axis x − = 0 + , x + > 0 for the right-moving partial waves and on the semi-axis x + = 0 + , x − > 0 for the left-moving partial waves. Therefore these formulas provide for each mode function the complete initial data on the lightcone τ = 0 + . However, the transformation to the quantum number ν (Fourier conjugate of the spatial rapidity η) introduces a non-analyticity in the time dependence, in the form of factors τ ±iν . For this reason, the numerical resolution of the Dirac equation should be started at a strictly positive proper time τ 0 > 0. Therefore, one should evolve the mode functions from the time τ = 0 + to the time τ 0 (see (120) This Green's formula in principle contains the Dirac propagator dressed by the background color field, which is not known analytically inside the forward light-cone. However, if one is interested only in very small propagation times Q s τ 0 1 one can use the bare Dirac propagator instead, the effects of the background field on the evolution of the spinors start becoming important only at Q s τ 0 1. But even with a bare propagator, the evaluation of eq. (120) is cumbersome. It turns out that this is not necessary, if the time τ 0 is chosen small enough.
For elementary plane waves, the basic formula to justify this is the following : where k ± = (M k / √ 2) e ±y and x ± = (τ 0 / √ 2) e ±η . This formula can be checked by an explicit calculation of the integrals on both sides, which can be expressed in terms of Hankel functions for the right hand side and in terms of Gamma functions in the left hand side, and by doing the Taylor expansion to first order of the Hankel functions. The interpretation of this formula is the following: • the left hand side is the sum of the left-and right-moving partial waves, evaluated as if the time was τ = 0 + (i.e. neglecting the evolution from τ to τ 0 ), • the right hand side contains the plane wave evolved to the time τ 0 (i.e. the result of using the Green's formula from τ to τ 0 ).
In other words, eq. (121) shows that it is legitimate to neglect the time evolution of the spinors in the forward light-cone, provided that the time obeys M k τ 0 1. This exercise shows that if the time τ 0 obeys this condition, it is sufficient to add up the two partial waves on the light-cone, and to apply the transformation y → ν to their sum. Let us end this appendix by an important remark regarding the condition M k τ 0 1: it must be satisfied for all the transverse masses M k = k 2 ⊥ + m 2 that can exist in the problem. In the lattice discretization of the Dirac equation, this implies that one must have τ 0 a ⊥ where a ⊥ is the transverse lattice spacing.

D Conserved inner product D.1 Definition and main properties
Let us consider a locally space-like surface 16 Σ. For every point y ∈ Σ, we can define an orthogonal vector such that n µ dy µ = 0 for any displacement dy µ around y ∈ Σ n 0 > 0 Given two spinors ψ and χ, one can define the following inner product on Σ, where d 3 S y is the 3-dimensional measure 17 on the surface Σ. One sees immediately that this inner product is Hermitean, The main property of the inner product defined in eq. (123) is that it is independent of the surface Σ if both ψ and χ are solutions of the same Dirac equation 18 , (This is true for any real valued background potential A µ .) From now on, we can thus drop the subscript Σ in our notation for this inner product. Note also that this inner product is gauge invariant, since it involves the product of a spinor and the Hermitean conjugate of another spinor at the same space-time point.
Let us give more explicit expressions of the inner product (123) for two important types of surface Σ. On a constant x 0 surface, it takes the following form, On a surface of constant proper time, the integration measure is d 3 S y = τ dηd 2 y ⊥ and the normal unit vector n µ has the following components Therefore, the inner product reads D.2 Inner product at x 0 = −∞ This conserved inner product can be used as a consistency check for the various analytic formulas that we have obtained for the mode functions in the section 4. The first step is to evaluate the inner product on the surface y 0 = −∞, where the mode functions ψ ± ksa are not yet modified by the background field. We get D.3 Inner product at τ = 0 + in LC gauge and y basis Now, let us use the eqs. (68) and their counterparts for the left-moving partial wave in order to check that the spinors (in light-cone gauge, and in the y basis) evolved to the forward light-cone are consistent with eqs. (129). First of all, from the definition of eq. (123), we immediately see that In words, on the right branch of the light-cone we need only the P − projection of the spinors, and their P + projection on the left branch of the light-cone.
The other projections do not contribute to the inner product evaluated on the light-cone 19 .
Adding the contributions of the two branches of the light-cone, we find the following expression for the inner product, where y, y are the momentum rapidities corresponding to the 3-momenta k, k . Using Gordon's identities, one sees that the imaginary part of the right hand side vanishes. Thanks to
Using the fact that exp(−yγ 0 γ 3 /2) acts on spinors as a boost of rapidity −y in the z direction, we have v † s (k ⊥ , y) e y P − +e −y P + v s (k ⊥ , y) = v † s (k ⊥ , y)e − y 2 γ 0 γ 3 e − y 2 γ 0 γ 3 v s (k ⊥ , y) Therefore, we obtain 20 which is consistent with the conservation of the inner product. A similar verification can be performed for the positive energy mode functions ψ + ksa .

D.4 Inner product at τ = 0 + in FS gauge and ν basis
We can perform the same check for the mode functions in the Fock-Schwinger gauge and ν basis given by eq. (74). Firstly, let us apply the transformation ψ k ⊥ νsa ≡ dy e iνy ψ k ⊥ ysa (138) to eqs. (129) in order to know what to expect for the inner product in the ν basis. This trivial calculation tells us that we should obtain Let us now calculate this inner product directly from eq. (74). The inner product on a surface of constant τ is given by eq. (128). It is convenient to absorb the factor τ exp(−ηγ 0 γ 3 ) by defining new spinors that are the original ones boosted to the comoving frame at the rapidity η, times a factor √ τ , ψ(τ, η, y ⊥ ) ≡ √ τ e − η 2 γ 0 γ 3 ψ(τ, η, y ⊥ ) .

E Abelian case: electron spectrum in QED
The initial value of the mode functions just above the forward light-cone, given in eq. (74), can also be used in the Abelian case. The analogous problem in QED would be that of the production of electrons in a high-energy collision of two large electrical charges Z 1 and Z 2 . When they collide, the electromagnetic field of the two charges can produce electron-positron pairs. The spectrum of the produced electrons can be calculated by a formalism which is very similar to the Color Glass Condensate considered in this paper. In this description, the two colliding charges are replaced by the electrical currents they carry along their trajectories. These currents are the source of electromagnetic fields, that can be obtained by solving the Maxwell's equations with sources : It is then this electromagnetic field that produces the charged fermions (in QED, this approach is known as the equivalent photon approximation.) Since QED is an Abelian gauge theory, it is much simpler than QCD in certain respects. Firstly, the Wilson lines U 1 and U 2 are simply complex valued phases, that can be commuted at will. Secondly, each of the colliding charge produces transverse E and B which are localized in a shockwave transverse to its trajectory. In the Fock-Schwinger gauge, the corresponding gauge potential is a pure gauge in the half-space located after the shockwave. But since Maxwell's equations are linear, the solution for the 2-charge problem is the sum of the solutions for individual nuclei, i.e. a sum of two pure gauge fields, which is itself a pure gauge. In QED, the fields are thus trivial inside the forward light-cone, unlike the QCD case.
The formula (74), that gives the mode functions just above the forward light-cone, is therefore sufficient to obtain in closed form the amplitude. One can evaluate it by computing the inner product of eq. (10) at an infinitesimal time τ → 0 + , since the evolution at τ > 0 is trivial. The only subtlety when doing so is that, since there is a non-zero pure gauge field in the forward lightcone, one should use a gauge rotated free positive energy spinor instead of ψ 0+ pσ , This free spinor must be first transformed as in eq. (70), In order to perform the integral over the momentum rapidity y, we need first to extract the y dependence hidden in the spinor u σ (p), u σ (p ⊥ , y) = e y 2 γ 0 γ 3 u σ (p ⊥ , 0) = e y 2 P + + e − y 2 P − u σ (p ⊥ , 0) .
Therefore, we have × dy e iνy e −iMpτ cosh(y) e y 2 P + + e − y 2 P − u σ (p ⊥ , 0) , where M p = p 2 ⊥ + m 2 denotes the transverse mass (in this equation, we have redefined the integration variable y − η → y). The result of the integration over y can be expressed in terms of Hankel functions, thanks to thanks to which we finally obtain This formula is identical to the eq. (52) in ref. [29]. Note that in order to recover this known result, it was crucial to gauge rotate the free spinor used in the projection, because the gauge field inside the forward light-cone is a non-zero pure gauge in the Fock-Schwinger gauge.