text:
# Coulomb corrections to 𝑒⁺𝑒⁻ production in ultra-relativistic nuclear collisions

## 1 Introduction

In the past three years, the problem of $`e^+e^{}`$ pair production induced by the collision of two nuclei has attracted a lot of interest, due partly to a series of papers showing that the Dirac equation of an electron can be solved exactly in the electromagnetic background field created by the two nuclei, in the limit where the two nuclei are ultra-relativistic . This solution was then used to derive an expression for the pair production cross-section. An unexpected consequence of this formula was that it lead to a pair production cross-section equal to its Born value (given by the diagram in figure 1). In other words, this result seems to indicate that all the corrections due to the multiple exchange of photons (called “Coulomb corrections” in the following) cancel exactly in the total cross-section.

This observation then prompted another series of papers on the same topic, basically saying that this property cannot be true. One of them  re-derived with new methods the classical result of Davies, Bethe and Maximon  about the fragmentation of a photon into a pair, induced by a electromagnetic background field. Then, this result was applied to pair production in nuclear collisions . In this perturbative calculation, the Coulomb corrections do not cancel and the resummed cross-section is about $`25\%`$ smaller than the Born cross-section at RHIC energy. Others claimed the discrepancy to be due to a breakdown of crossing symmetry . Another paper  attributed the contradiction to a mishandling of ill-defined integrals, and claimed that the result of Ivanov et al. can be recovered from the solution of the Dirac equation if a proper regularization is used. But those papers apparently focused on technical details, without first asking if an agreement between  and  should be expected at all. In other words, are those two approaches calculating the same physical quantity ? Given the fact that one obtains a retarded amplitude by solving the Dirac equation, while the perturbative approach calculates a time-ordered amplitude, such an assumption is highly non-trivial, and in fact very unlikely.

In this paper, we re-examine from first principles the problem of pair production by two colliding nuclei. Our main focus is on showing that the two approaches used to study pair production in the peripheral collision of two ultra-relativistic nuclei are in fact calculating different physical quantities. In particular, we give the relationship between several physical quantities and the appropriate correlators. The exclusive cross-section of single pair production is related to a time-ordered (or Feynman) correlator, and is the object usually calculated in perturbation theory. On the other hand, the inclusive cross-section of pair production can be related to retarded amplitudes, which are the ones accessible by solving the Dirac equation with retarded boundary conditions. Therefore, the main discrepancy stems from the fact that the two approaches are calculating different quantities, despite giving them the same name; and the two approaches do not have to give the same answer.

We then pursue this analysis by studying the retarded and Feynman propagators in the background field of two ultra-relativistic nuclei, and show that these two types of propagators have dramatically different perturbative expansions: the retarded propagator can be obtained in closed form in the ultra-relativistic limit, as opposed to the Feynman propagator.

We also resolve a puzzle related to the unitarity of the pair production probabilities calculated in perturbation theory, which seems to give probabilities larger than 1 . We show that this calculation must include as a corrective factor the vacuum-to-vacuum transition probability, and that this factor precisely restores unitarity.

The structure of this paper is as follows: after setting up some notations and details about the model in section 2, we recall in section 3 some relations between observable quantities and various types of correlators. We also point out some uncommon aspects of the Feynman rules that are due to the fact that the background field is non-static.

We study in section 4 the propagators in the background field of only one nucleus. This case can be solved exactly in the ultra-relativistic limit, and serves as a starting point for the more complicated case of two nuclei.

In section 5, we begin with a formal solution to the problem of two nuclei, which enables one to write both the retarded and Feynman propagators in terms of the scattering matrices associated to the individual nuclei. Then, we explain why the ultra-relativistic limit simplifies this formula in the case of the retarded propagator, but does not help for the Feynman one.

In section 6, we focus on the problem of the unitarity of the perturbative predictions, and show how it is resolved. We also show that the multiplicity distribution is not exactly Poissonian.

Section 7 discusses two strategies one can use to get an approximate answer for the Feynman propagator, and hence for the pair production amplitude.

Finally, section 8 is devoted to concluding remarks.

## 2 Ultra-relativistic limit. Notations

In this paper, we consider two nuclei of electric charge $`Z_1e`$ and $`Z_2e`$ respectively, colliding with an impact parameter $`b`$. We study the problem in the frame where the two nuclei have opposite velocities. For the sake of definiteness, we decide that the nucleus $`Z_1`$ is moving in the positive $`z`$ direction, while the other nucleus is moving in the negative $`z`$ direction. We denote by $`A_1^\mu (t,x)`$ and $`A_2^\mu (t,x)`$ the electromagnetic vector potentials created by those nuclei. Because the superposition theorem holds in QED111This property is not true in QCD. The problem of $`q\overline{q}`$ pair production is therefore complicated by the necessity of first finding the color vector potential due to the two nuclei., the total vector potential is $`A^\mu (t,x)=A_1^\mu (t,x)+A_2^\mu (t,x)`$.

In the ultra-relativistic limit, we can assume that the two nuclei do not recoil, and simply treat them as a classical background field. The trajectories of the two nuclei are respectively $`z=\pm t`$ ($`c=1`$), $`x_{}=\pm b/2`$, and there is a gauge in which $`A_1^\mu `$ and $`A_2^\mu `$ have the following form :

$`A_1^\mu (t,x)=Z_1e\delta (v_+x)v_+^\mu \mathrm{ln}\left({\displaystyle \frac{(x_{}b/2)^2}{b^2}}\right),`$
$`A_2^\mu (t,x)=Z_2e\delta (v_{}x)v_{}^\mu \mathrm{ln}\left({\displaystyle \frac{(x_{}+b/2)^2}{b^2}}\right),`$ (1)
where $`v_\pm ^\mu (1,0,0,\pm 1)/\sqrt{2}`$. This choice of gauge makes the following fact obvious: the electromagnetic field of an ultra-relativistic nucleus is confined in its transverse plane due to Lorentz contraction. In the ultra-relativistic approximation, we neglect any modification of the trajectories of the two nuclei. Therefore, one can see this problem as a field theory for fermions in a classical electromagnetic background field. An essential property of the field generated by two colliding nuclei is that there is no frame in which it is time-independent (static). This property has important consequences on the field theory describing fermions in this background; one of them is that pairs can be created222This is to be contrasted with the case of a single nucleus. In the frame of that nucleus, the background field is static, and it is well know that a single nucleus moving on a straight line does not produce pairs (except via the non-perturbative Schwinger mechanism of vacuum instability , which we ignore here)..

However, one should realize that the limit of infinite momentum for the nuclei is only an approximation of the real physical situation. Indeed, the potentials of Eq. (1) cannot be expected to describe the motion of a lepton comoving with one of the nuclei. This restriction also appears in the work of . In particular,  shows that the gauge transformation that leads to Eq. (1) is well-defined only if the lepton is not comoving with one of the nuclei. There is another way to see the problem with Eq. (1): this background potential is short-ranged in the $`z`$ and $`t`$ directions, while the potential before the gauge transformation was long-ranged. The absence of interaction at asymptotic times allows the construction of plane-wave asymptotic states, and makes the usual formulation of reduction theory applicable. But it is also clear that those interactions cannot be removed by a gauge transformation for a comoving lepton, the states of which are described by distorted waves. In conclusion, those potentials are valid only if the typical interaction time between the lepton and the nuclei goes like $`1/\gamma `$, i.e. for leptons produced in the mid-rapidity region.

In addition, the description of the background field by a delta function is probably not accurate at large transverse distances if the Lorentz factor $`\gamma `$ is finite. Indeed, at large distances, the lines of force of the electric field should have some curvature and depart from the transverse plane. Therefore, Eq. (1) is expected to break down also at large impact parameter $`b`$ or if the transverse separation between the lepton and a nucleus is large. This problem is expected to show up via infrared divergences in the transverse momentum integrals. In , it was argued that those divergences are regularized by an effective cutoff of order $`\omega /\gamma `$. Although this argument cannot tell precisely what this cutoff should be, it is obvious from the derivation of the ultra-relativistic limit of the potentials that it is completely determined by the kinematics and is therefore independent of $`Z_{1,2}\alpha `$ (the atomic numbers and coupling constants appear only in the numerator of the Coulomb potential). Therefore, the uncertainty in the cutoff can affect only the overall normalization of integrated cross-sections, but not its dependence on atomic numbers.

In view of those caveats, one should take with a grain of salt the results obtained in this model for completely integrated (over the pair phase-space and over impact parameter) cross-sections, because there are regions in such an integral where the above approximation may not be appropriate.

One should realize that due to the large electric charge $`Ze`$ of a nucleus like gold, the addition of an extra photon connecting the electron line to a nucleus brings a factor $`Z\alpha `$ ($`\alpha e^2/4\pi `$), which may be of order $`1`$. Therefore, we want to include as much as possible of these “perturbative” corrections. On the contrary, photons coupling an electron line to itself (or to another electron line) yield only a factor of $`\alpha `$, which is indeed a small correction. Therefore, we also neglect the interactions of the electrons with dynamical photons, and keep only the classical photon background. The Lagrangian density for the electron field $`\psi (x)`$ is therefore

$$\overline{\psi }(x)(i/_xe/A(x)m)\psi (x)$$
(2)
in this model.

In the following, we make extensive use of the light-cone coordinates. For any 4-vector $`x^\mu `$, we define:

$$x^\pm \frac{x^0\pm x^3}{\sqrt{2}},$$
(3)
and denote by $`x_{}`$ the transverse part of the 3-vector $`x`$. With these notations, the invariant norm of $`x^\mu `$ is $`x^2=2x^+x^{}x_{}^{}{}_{}{}^{2}`$, and the scalar product of $`k^\mu `$ and $`x^\mu `$ is $`kx=k^+x^{}+k^{}x^+k_{}x_{}`$. The invariant measure $`d^4x`$ becomes $`d^4x=dx^+dx^{}d^2x_{}`$. Note also that $`x^\pm =xv_{}`$ with the $`v_\pm `$ defined above.

## 3 Reduction formulae and Feynman rules

In this section, we relate observable quantities to correlation functions of the fermionic field operator. In fact, these considerations do not depend on the nature of the background field, and rely only on the Lagrangian of Eq. (2). In particular, the formulae of this section are independent of the ultra-relativistic approximation.

### 3.1 Amplitude to produce one pair

#### 3.1.1 Reduction formula

The amplitude to produce one $`e^+e^{}`$ pair is

$$M_1(p,q)e^+(p)e^{}(q){}_{\mathrm{out}}{}^{}|0_{\mathrm{in}}=0_{\mathrm{out}}|d_{\mathrm{out}}(p)b_{\mathrm{out}}(q)|0_{\mathrm{in}},$$
(4)
where $`d_{\mathrm{out}}^{}(p)`$ (resp. $`b_{\mathrm{out}}^{}(p)`$ ) is the operator that creates a positron (resp. an electron) of 3-momentum $`p`$ in the final state. At this stage, it is very important to carefully distinguish the in- and out- states and operators333For instance, $`b_{\mathrm{out}}`$ does not annihilate the in-vacuum: $`b_{\mathrm{out}}|0_{\mathrm{in}}0`$, while $`b_{\mathrm{out}}|0_{\mathrm{out}}=0`$. It is precisely this fact that makes the pair production amplitude non-zero in a non-static background field..

Making use of the following relations between the annihilation/creation operators444For the sake of brevity, we do not write explicitly the spin indices on spinors and creation/annihilation operators. In formulae for cross-sections, it is implicitly assumed that we sum over the spin of the final state particles. and the field itself (see , page 61):

$`b_{\mathrm{out}}(q)={\displaystyle d^3x\overline{u}(q)\gamma ^0\psi _{\mathrm{out}}(t,x)e^{iqx}},`$
$`d_{\mathrm{out}}(p)={\displaystyle d^3x\overline{\psi }_{\mathrm{out}}(t,x)\gamma ^0v(p)e^{ipx}},`$ (5)
where $`p_0`$ and $`q_0`$ are, respectively, $`\sqrt{(p^2+m^2)}`$ and $`\sqrt{(q^2+m^2)}`$ (the time $`x^0=t`$ that shows up in these formulae is irrelevant, and disappears in the course of the calculation), one can show by standard manipulations (, pages 205-207) that the pair production amplitude defined above can be related to a 2-point time-ordered correlator by the following reduction formula:

$`0_{\mathrm{out}}|d_{\mathrm{out}}(p)b_{\mathrm{out}}(q)|0_{\mathrm{in}}=\left[{\displaystyle \frac{i}{\sqrt{𝒵_2}}}\right]^2{\displaystyle d^4xd^4y}`$
$`\times e^{iqx}\overline{u}(q)(i\underset{x}{\overset{}{/}}m)0_{\mathrm{out}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}(i\underset{y}{\overset{}{/}}+m)v(p)e^{ipy},`$ (6)
where the arrows indicate on which side the derivatives act (preventing their action on the exponentials), and where $`𝒵_2`$ is the wave function renormalization factor for an electron555With the Lagrangian of Eq. (2), this factor is equal to $`1`$.. The important point to note here is that the average value of the time-ordered product is taken between the in- and out- vacua, which are different states. To expand a little on this, let us add that we necessarily have $`|0_{\mathrm{out}}|0_{\mathrm{in}}|^2<1`$ in background fields that can produce pairs, because of unitarity (see section 6 for the role played by $`0_{\mathrm{out}}|0_{\mathrm{in}}`$ in issues related to unitarity).

#### 3.1.2 Perturbative expansion of $`0_{\mathrm{out}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$

We need to calculate the correlator $`0_{\mathrm{out}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$. From that, Eq. (6) tells us how to obtain the pair production amplitude: amputate the external legs of the correlator, take its Fourier transform, and insert the result between the appropriate spinors.

The Feynman rules to calculate perturbatively this correlator can be obtained by switching to the interaction picture. The Heisenberg field can be expressed in terms of the field in the interaction representation, via the following relation:

$$\psi (x)U(t__I,x^0)\psi __I(x)U(x^0,t__I),$$
(7)
where $`t__I`$ is the time at which the Heisenberg and interaction pictures coincide (ultimately, we will take $`t__I`$ to $`\mathrm{}`$), and where the evolution operator $`U`$ is related to the interaction part $`_{\mathrm{int}}e\overline{\psi }(x)/A(x)\psi (x)`$ of the Lagrangian by

$$U(t_2,t_1)=\mathrm{P}_{12}\mathrm{exp}i_{t_1}^{t_2}d^4x_{\mathrm{int}}(\psi __I(x)),$$
(8)
where $`\mathrm{P}_{12}`$ is an ordering operator along the path connecting $`t_1`$ to $`t_2`$ (ordinary time-ordering if $`t_1<t_2`$, and reverse time-ordering if $`t_1>t_2`$). Using this transformation, and taking $`t__I\mathrm{}`$ (in this limit, $`\psi __I\psi _{\mathrm{in}}`$), we have

$`0_{\mathrm{out}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}=<0_{\mathrm{out}}|U(\mathrm{},+\mathrm{})\mathrm{T}\overline{\psi }_{\mathrm{in}}(y)\psi _{\mathrm{in}}(x)`$
$`\times \mathrm{exp}i{\displaystyle _{\mathrm{}}^+\mathrm{}}d^4x_{\mathrm{int}}(\psi _{\mathrm{in}}(x))|0_{\mathrm{in}}>.`$ (9)
Noticing that $`0_{\mathrm{out}}|U(\mathrm{},+\mathrm{})=0_{\mathrm{in}}|`$, we have simply:

$$0_{\mathrm{out}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}=<0_{\mathrm{in}}\left|\mathrm{T}\overline{\psi }_{\mathrm{in}}(y)\psi _{\mathrm{in}}(x)\mathrm{exp}i_{\mathrm{}}^+\mathrm{}d^4x_{\mathrm{int}}(\psi _{\mathrm{in}}(x))\right|0_{\mathrm{in}}>.$$
(10)
Everything being now expressed in terms of in-fields and in-states, the right hand side of the last equation can be evaluated by the standard perturbative expansion of the exponential.

This perturbative expansion has, however, a peculiarity due to the fact that the background field is non-static: the vacuum-vacuum diagrams (i.e. diagrams without any external legs) do not cancel, and their sum is not a phase. On the contrary, in a conventional field theory, $`0_{\mathrm{out}}|0_{\mathrm{in}}`$ is just an irrelevant phase, and one takes advantage of this fact to divide the r.h.s. of the previous equation by $`0_{\mathrm{in}}|U(+\mathrm{},\mathrm{})|0_{\mathrm{in}}`$. Then, one can show that this denominator cancels (, pages 266-267) order by order the vacuum-vacuum diagrams666This step is not mandatory though. One can live with the vacuum-vacuum diagrams and notice that they add up to a pure phase, so that they always drop out in the calculation of cross sections, even if they are present in the amplitude.. This trick cannot be used here due to the fact that $`0_{\mathrm{out}}|0_{\mathrm{in}}`$ is not a phase in the present problem . In fact, this is deeply rooted in the property that the background can produce particles.

We can keep this complication aside for a while by just writing:

$`0_{\mathrm{out}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}=0_{\mathrm{in}}|U(+\mathrm{},\mathrm{})|0_{\mathrm{in}}`$
$`\times {\displaystyle \frac{0_{\mathrm{in}}\left|\mathrm{T}\overline{\psi }_{\mathrm{in}}(y)\psi _{\mathrm{in}}(x)\mathrm{exp}i_{\mathrm{}}^+\mathrm{}d^4x_{\mathrm{int}}(\psi _{\mathrm{in}}(x))\right|0_{\mathrm{in}}}{0_{\mathrm{in}}|U(+\mathrm{},\mathrm{})|0_{\mathrm{in}}}},`$ (11)
so that the fraction on the right hand side has a perturbative expansion where the vacuum-vacuum diagrams do cancel. In fact, this fraction is nothing but the Feynman propagator

$$G_F(x,y)\frac{0_{\mathrm{out}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}}{0_{\mathrm{out}}|0_{\mathrm{in}}}$$
(12)
of an electron in the electromagnetic background field. Its perturbative expansion is the usual one: it is obtained by inserting the external potential $`ie/A(x)`$ on chains of free Feynman propagators $`G__F^0(x,y)0_{\mathrm{in}}|\mathrm{T}\overline{\psi }_{\mathrm{in}}(y)\psi _{\mathrm{in}}(x)|0_{\mathrm{in}}`$. In particular, the perturbative expansion for the full propagator can be generated by the following Lippmann-Schwinger equation:

$$G(x,y)=G^0(x,y)ied^4zG^0(x,z)/A(z)G(z,y).$$
(13)
Note that this equation is equally valid for the Feynman and for the retarded propagator (provided the free propagator $`G^0`$ is chosen accordingly). This remark will become important later when we also need to study the retarded propagator.

It turns out to be more convenient to work in Fourier space. If we define777We do not use a distinct symbol for Fourier transforms, as the context always enables to tell the difference.

$$G(x,y)\frac{d^4p}{(2\pi )^4}\frac{d^4q}{(2\pi )^4}e^{iqx}e^{ipy}G(q,p),$$
(14)
then the Lippmann-Schwinger equation becomes

$$G(q,p)=(2\pi )^4\delta (pq)G^0(p)ieG^0(q)\frac{d^4k}{(2\pi )^4}/A(k)G(q+k,p),$$
(15)
with $`G^0(p)i/(/pm)`$ (and an unspecified $`iϵ`$ prescription depending on whether we are studying the retarded or the Feynman propagator) and

$$A^\mu (k)d^4xe^{ikx}A^\mu (x).$$
(16)
If one introduces the interacting part $`𝒯_F`$ of the Feynman propagator by the relation

$$G_F(q,p)=(2\pi )^4\delta (pq)G_F^0(p)+G_F^0(q)𝒯_F(q,p)G_F^0(p),$$
(17)
then Eqs. (6) and (11) lead to the following expression for the probability to produce exactly one pair in a collision at impact parameter $`b`$:

$$P_1=|0_{\mathrm{out}}|0_{\mathrm{in}}|^2\frac{d^3q}{(2\pi )^32\omega _q}\frac{d^3p}{(2\pi )^32\omega _p}\left|\overline{u}(q)𝒯_F(q,p)v(p)\right|^2.$$
(18)
Note that this formula has been derived independently of the details of the background field, and is therefore completely general.

### 3.2 Average number of produced pairs

#### 3.2.1 Expression as a correlator

The approach based on solving the Dirac equation was motivated by the papers , where a formula giving the average number of pairs produced per collision is derived in terms of retarded amplitudes only. We present here a justification of this formula in the field-theoretical framework we are following in this paper.

Let us start from the expression of the average number of pairs $`\overline{n}`$ as the sum

$$\overline{n}=\underset{n=1}{\overset{+\mathrm{}}{}}nP_n=\underset{n=0}{\overset{+\mathrm{}}{}}(n+1)P_{n+1},$$
(19)
where we denote by $`P_n`$ the probability to produce exactly $`n`$ pairs in a collision at impact parameter $`b`$. Making explicit what this probability is, we find

$`\overline{n}`$ $`={\displaystyle \underset{n=0}{\overset{+\mathrm{}}{}}}(n+1){\displaystyle \frac{1}{(n+1)!^2}}{\displaystyle \underset{i=1}{\overset{n+1}{}}\frac{d^3p_i}{(2\pi )^32\omega _{p_i}}\frac{d^3q_i}{(2\pi )^32\omega _{q_i}}}`$ (20)
$`\times \left|<0_{\mathrm{out}}\left|d_{\mathrm{out}}(p_1)b_{\mathrm{out}}(q_1)\mathrm{}d_{\mathrm{out}}(p_{n+1})b_{\mathrm{out}}(q_{n+1})\right|0_{\mathrm{in}}>\right|^2,`$
where we denote $`\omega _p\sqrt{(p^2+m^2)}`$. Singling out one of the electrons (say the one with momentum $`q_{n+1}`$, which we call simply $`q`$) and expanding the squared modulus, we can write:

$`\overline{n}`$ $`={\displaystyle \frac{d^3q}{(2\pi )^32\omega _q}<0_{\mathrm{in}}|b_{\mathrm{out}}^{}(q)}`$ (21)
$`\times {\displaystyle \underset{n=0}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}{\displaystyle \frac{1}{(n+1)!}}{\displaystyle }{\displaystyle \underset{i=1}{\overset{n+1}{}}}{\displaystyle \frac{d^3p_i}{(2\pi )^32\omega _{p_i}}}{\displaystyle \underset{j=1}{\overset{n}{}}}{\displaystyle \frac{d^3q_j}{(2\pi )^32\omega _{q_j}}}`$
$`\left|d_{\mathrm{out}}^{}(p_1)b_{\mathrm{out}}^{}(q_1)\mathrm{}d_{\mathrm{out}}^{}(p_n)b_{\mathrm{out}}^{}(q_n)d_{\mathrm{out}}^{}(p_{n+1})\right|0_{\mathrm{out}}>`$
$`<0_{\mathrm{out}}\left|d_{\mathrm{out}}(p_1)b_{\mathrm{out}}(q_1)\mathrm{}d_{\mathrm{out}}(p_n)b_{\mathrm{out}}(q_n)d_{\mathrm{out}}(p_{n+1})\right|`$
$`\times b_{\mathrm{out}}(q)|0_{\mathrm{in}}>.`$
Noticing now that the three intermediate lines are the identity operator on the subspace of states with electric charge $`+e`$, we find:

$$\overline{n}=\frac{d^3q}{(2\pi )^32\omega _q}<0_{\mathrm{in}}\left|b_{\mathrm{out}}^{}(q)b_{\mathrm{out}}(q)\right|0_{\mathrm{in}}>.$$
(22)
This formula simply tells that in order to count the number of pairs produced if the initial state is the vacuum, it is sufficient to count the number of electrons in the final state. Had we decided to single out a positron annihilation operator in Eq. (21), we would have obtained instead:

$$\overline{n}=\frac{d^3p}{(2\pi )^32\omega _p}<0_{\mathrm{in}}\left|d_{\mathrm{out}}^{}(p)d_{\mathrm{out}}(p)\right|0_{\mathrm{in}}>.$$
(23)
One can note that if the background potential has time-reversal symmetry, we have:

$$\overline{n}=\frac{d^3q}{(2\pi )^32\omega _q}<0_{\mathrm{out}}\left|b_{\mathrm{in}}^{}(q)b_{\mathrm{in}}(q)\right|0_{\mathrm{out}}>=\frac{d^3p}{(2\pi )^32\omega _p}<0_{\mathrm{out}}\left|d_{\mathrm{in}}^{}(p)d_{\mathrm{in}}(p)\right|0_{\mathrm{out}}>.$$
(24)
There are also reduction formulae for those correlators, that give for instance

$`\overline{n}`$ $`={\displaystyle \frac{d^3q}{(2\pi )^32\omega _q}\left[\frac{i}{\sqrt{𝒵_2}}\right]^2d^4xd^4y}`$ (25)
$`\times e^{iqx}\overline{u}(q)(i\underset{x}{\overset{}{/}}m)0_{\mathrm{in}}|\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}(i\underset{y}{\overset{}{/}}m)u(q)e^{iqy}.`$
Therefore, in order to calculate the average number of pairs, we need now the ordinary product of two fields, averaged with the initial vacuum $`|0_{\mathrm{in}}`$ on both sides.

#### 3.2.2 Perturbative expansion of $`0_{\mathrm{in}}|\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$

In order to switch to the interaction representation, it is easier to start with $`0_{\mathrm{in}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$ for which one can write directly

$`0_{\mathrm{in}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}=<0_{\mathrm{in}}|U(\mathrm{},+\mathrm{})\mathrm{T}\overline{\psi }_{\mathrm{in}}(y)\psi _{\mathrm{in}}(x)`$
$`\times \mathrm{exp}i{\displaystyle _{\mathrm{}}^+\mathrm{}}d^4x_{\mathrm{int}}(\psi _{\mathrm{in}}(x))|0_{\mathrm{in}}>.`$ (26)
This time, one cannot get rid of the $`U(\mathrm{},+\mathrm{})`$ in the right hand side, but there is a standard trick to incorporate it in the perturbative expansion . For that, one has to introduce a contour $`𝒞`$ going from $`\mathrm{}`$ to $`+\mathrm{}`$ just above the real-time axis, and then back from $`+\mathrm{}`$ to $`\mathrm{}`$ below the real axis:

One can then check that

$$0_{\mathrm{in}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}=<0_{\mathrm{in}}\left|\mathrm{P}\overline{\psi }_{\mathrm{in}}(y)\psi _{\mathrm{in}}(x)\mathrm{exp}i_𝒞d^4x_{\mathrm{int}}(\psi _{\mathrm{in}}(x))\right|0_{\mathrm{in}}>.$$
(27)
In the previous formula, $`\mathrm{P}`$ stands for an ordering of operators along the path $`𝒞`$, identical to the usual time-ordering $`\mathrm{T}`$ on the upper branch of the contour (and with the convention that points on the lower branch have a “larger” time than points on the upper branch). We therefore have formally similar Feynman rules for the perturbative expansion of $`0_{\mathrm{in}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$, except that the time integrations at each insertion of the external potential are performed on the contour $`𝒞`$.

It is customary to write this formalism in matrix form  (in this section, objects denoted by a boldface letter are $`2\times 2`$ matrices), by splitting the free propagator in four components according to where $`x^0`$ and $`y^0`$ lie on $`𝒞`$ (upper or lower branch):

$$G^0(x,y)\left(\begin{array}{cc}G_{++}^0(x,y)& G_+^0(x,y)\\ G_+^0(x,y)& G_{}^0(x,y)\end{array}\right),$$
(28)
with

$`G_{++}^0(x,y)0_{\mathrm{in}}|\mathrm{T}\overline{\psi }_{\mathrm{in}}(y)\psi _{\mathrm{in}}(x)|0_{\mathrm{in}},`$
$`G_{}^0(x,y)0_{\mathrm{in}}|\overline{\mathrm{T}}\overline{\psi }_{\mathrm{in}}(y)\psi _{\mathrm{in}}(x)|0_{\mathrm{in}},`$
$`G_+^0(x,y)0_{\mathrm{in}}|\overline{\psi }_{\mathrm{in}}(y)\psi _{\mathrm{in}}(x)|0_{\mathrm{in}},`$
$`G_+^0(x,y)0_{\mathrm{in}}|\psi _{\mathrm{in}}(x)\overline{\psi }_{\mathrm{in}}(y)|0_{\mathrm{in}},`$ (29)
where $`\overline{\mathrm{T}}`$ is the reverse time-ordering operator. The correlator $`0_{\mathrm{in}}|\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$ we need for the average number of produced electrons is the $`+`$ component of the exact matrix propagator.

One needs also to give a matrix structure to the external potential $`/A(x)`$:

$$/A(x)\tau _3/A(x),$$
(30)
where $`\tau _3\mathrm{Diag}(1,1)`$ is the third Pauli matrix888The $``$ sign for the $``$ component of $`/`$$`A`$ comes from the fact that the integral over the lower branch of the time contour goes from $`+\mathrm{}`$ back to $`\mathrm{}`$.. Then, at each vertex, the rule is to integrate only over the ordinary time axis, and to multiply the matrices corresponding to the propagators and potential insertions (in the order they appear in the Feynman diagram). For instance, the Lippmann-Schwinger equation for the exact matrix propagator in Fourier space is

$$G(q,p)=(2\pi )^4\delta (pq)G^0(p)ieG^0(q)\frac{d^4k}{(2\pi )^4}/A(k)G(q+k,p).$$
(31)
Note that in Fourier space, the expression of the free matrix propagator is999The expression of those propagators can be found in  (Eqs. (3.93)). To apply them here, the following substitutions must be made: $`\sigma =0`$, $`n(p_0)=0`$. Also,  is using different notations for the indices labelling the two branches of the contour: $`1+`$, $`2`$. :

$`G_{++}^0(p)=i{\displaystyle \frac{/p+m}{p^2m^2+iϵ}}=[G_{}^0(p)]^{},`$
$`G_+^0(p)=2\pi \theta (p_0)(/p+m)\delta (p^2m^2),`$
$`G_+^0(p)=2\pi \theta (+p_0)(/p+m)\delta (p^2m^2).`$ (32)
As it stands, the expansion of Eq. (31) has a very intricate matrix structure. However, it can be simplified by applying a “rotation”  on the previously defined matrices. Let us define

$`G_U^0(p)U(p)G^0(p)U^^T(p),`$
$`G_U(q,p)U(q)G(q,p)U^^T(p),`$
$`/A_U(q,q+k)U^^T{}_{}{}^{1}(q)/A(k)U^1(q+k),`$ (33)
where $`U`$ is an invertible matrix. The Lippmann-Schwinger equation hardly changes in the rotated formalism,

$$G_U(q,p)=(2\pi )^4\delta (pq)G_U^0(p)ieG_U^0(q)\frac{d^4k}{(2\pi )^4}/A_U(q,q+k)G_U(q+k,p),$$
(34)
but there are some choices of $`U`$ that simplify significantly the free matrix propagator. A convenient choice is

$$U(p)\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right),$$
(35)
which leads to the free propagator

$$G_U^0(p)=\left(\begin{array}{cc}0& G_A^0(p)\\ G_R^0(p)& G_S^0(p)\end{array}\right),$$
(36)
with the definitions (the index “$`S`$” stands for “on-shell”):

$`G_R^0(p)i{\displaystyle \frac{/p+m}{p^2m^2+ip_0ϵ}},G_A^0(p)i{\displaystyle \frac{/p+m}{p^2m^2ip_0ϵ}}`$
$`G_S^0(p)2\pi (/p+m)\delta (p^2m^2).`$ (37)
In this transformation, the external potential becomes

$$/A_U(q,q+k)=/A(k)\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right).$$
(38)
Therefore, we have a simplification in the propagator which has now a vanishing component, the non-zero components being the free retarded and advanced propagators and the very simple on-shell piece $`G_S^0`$. That this transformation helps to resum the perturbative expansion is readily seen by computing the building block that will be iterated in the expansion:

$$/A_U(q,q+k)G_U^0(q+k)=/A(k)\left[\left(\begin{array}{cc}G_R^0(q+k)& 0\\ 0& G_A^0(q+k)\end{array}\right)+G_S^0(q+k)\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)\right].$$
(39)
Indeed, the fact that this object is the sum of a diagonal matrix and a nilpotent matrix leads trivially to the following formula for a term with $`n`$ insertions of $`/`$$`A`$:

$`G_U^0(q_0)\left[/A_U(q_0,q_1)G_U^0(q_1)\right]\mathrm{}\left[/A_U(q_{n1},q_n)G_U^0(q_n)\right]=`$
$`=\left(\begin{array}{c}0G_A^0(q_0)/A(q_1q_0)G_A^0(q_1)\mathrm{}/A(q_nq_{n1})G_A^0(q_n)\\ G_R^0(q_0)/A(q_1q_0)G_R^0(q_1)\mathrm{}/A(q_nq_{n1})G_R^0(q_n)0\end{array}\right)`$
$`+{\displaystyle \underset{i=0}{\overset{n}{}}}G_R^0(q_0)/A(q_1q_0)\mathrm{}G_R^0(q_{i1})/A(q_iq_{i1})G_S^0(q_i)`$
$`\times /A(q_{i+1}q_i)G_A^0(q_{i+1})/A(q_{i+2}q_{i+1})\mathrm{}G_A^0(q_n)\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right).`$ (40)
At this stage, it is straightforward to sum over the number $`n`$ of insertions of the external potential, to obtain the exact propagator in the new basis (set $`qq_0`$ and $`pq_n`$ in the previous formula):

$$G_U(q,p)=\left(\begin{array}{cc}0& G_A(q,p)\\ G_R(q,p)& G_S(q,p)\end{array}\right),$$
(41)
where $`G_{R,A}(q,p)`$ are the exact retarded and advanced propagators, and where

$`G_S(q,p)`$ $`{\displaystyle \frac{d^4k}{(2\pi )^4}2\pi \delta (k^2m^2)}`$ (42)
$`\times G_R(q,k)G_R^0{}_{}{}^{1}(k)(/k+m)G_A^0{}_{}{}^{1}(k)G_A(k,p).`$
At this stage, we have an exact expression for the Fourier transform $`G_+(q,p)`$ of $`0_{\mathrm{in}}|\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$, which reads

$$G_+(q,p)=\frac{1}{2}\left[G_A(q,p)G_R(q,p)+G_S(q,p)\right].$$
(43)
Using Eq. (42), we see that this correlator can be written in closed form in terms of $`G_R`$ and $`G_A`$101010It is now clear why Eq. (31) had a very complicated matrix structure: trying to solve directly this matrix equation amounts to write the exact retarded and advanced propagators in terms of bare time-ordered propagators, which is not easy in practice. This explains why the rotation leading to the retarded/advanced basis simplified a lot this calculation.. According to Eq. (25), $`\overline{n}`$ is obtained by amputating the external legs of $`G_+(q,q)`$. To that effect, it is convenient to introduce the interacting parts $`𝒯_R(q,p)`$ and $`𝒯_A(q,p)`$ of the exact retarded and advanced propagators, defined by

$`G_R(q,p)=(2\pi )^4\delta (pq)G_R^0(p)+G_R^0(q)𝒯_R(q,p)G_R^0(p),`$
$`G_A(q,p)=(2\pi )^4\delta (pq)G_A^0(p)+G_A^0(q)𝒯_A(q,p)G_A^0(p),`$ (44)
which leads to

$`\overline{n}={\displaystyle \frac{1}{2}}{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^22\omega _q}}\overline{u}(q)[𝒯_A(q,q)𝒯_R(q,q)`$
$`+{\displaystyle }{\displaystyle \frac{d^4k}{(2\pi )^4}}2\pi \delta (k^2m^2)𝒯_R(q,k)(/k+m)𝒯_A(k,q)]u(q).`$ (45)
By using the Lippmann-Schwinger equations for the retarded and advanced propagators, one can check that111111This relation can be seen as a form of the optical theorem.

$`𝒯_R(q,p)𝒯_A(q,p)`$ $`={\displaystyle \frac{d^4k}{(2\pi )^4}2\pi [\theta (k_0)\theta (k_0)]\delta (k^2m^2)}`$ (46)
$`\times 𝒯_R(q,k)(/k+m)𝒯_A(k,p).`$
From there, we can further simplify $`\overline{n}`$, and obtain:

$`\overline{n}`$ $`={\displaystyle \frac{d^3q}{(2\pi )^22\omega _q}\frac{d^4k}{(2\pi )^4}2\pi \theta (k_0)\delta (k^2m^2)}`$ (47)
$`\times \overline{u}(q)𝒯_R(q,k)(/k+m)𝒯_A(k,q)u(q).`$
Changing $`kk`$ in the above formula, and integrating over $`k_0`$ thanks to the $`\delta (k^2m^2)`$ permits to simplify even more the expression of $`\overline{n}`$. Noticing that

$`𝒯_A(k,q)=\left[𝒯_R(q,k)\right]^{},`$ (48)
$`/km={\displaystyle \underset{\mathrm{spin}}{}}v(k)\overline{v}(k),`$ (49)
we finally obtain the following compact expression for the average number of pairs produced per collision at impact parameter $`b`$:

$$\overline{n}=\frac{d^3q}{(2\pi )^32\omega _q}\frac{d^3p}{(2\pi )^32\omega _p}\left|\overline{u}(q)𝒯_R(q,p)v(p)\right|^2.$$
(50)
Therefore, we have obtained for the average number of produced pairs a formula similar to Eq. (18), with however two major differences: the factor $`|0_{\mathrm{out}}|0_{\mathrm{in}}|^2`$ is not present here, and the formula for $`\overline{n}`$ involves the retarded propagator instead of the Feynman propagator.

### 3.3 Discussion of earlier literature

Equations (18) and (50) show clearly the main difference between the approach based on the Dirac equation and the perturbative approach concerning the “pair production cross-section”. Indeed, what is called “cross-section” in :

$$\sigma _Td^2b\overline{n},$$
(51)
is the inclusive cross-section121212Integrating over the entire range of impact parameters might include contributions from a region where our model is not expected to be valid. Cutoffs are implicitly understood in this integral, in order to make it finite. of pair production (obtained by counting all the pairs produced); while according to the nature of the diagrams considered in , the object they called “cross-section” is in fact131313Up to the factor $`|0_{\mathrm{out}}|0_{\mathrm{in}}|^2`$, which they seem to have overlooked.

$$\sigma _1d^2bP_1,$$
(52)
i.e. the exclusive cross-section to produce exactly one pair (measured by counting only those events that produce exactly one pair). It is now obvious that one should not expect an agreement between the two approaches, since they are calculating different observables. In fact, it is obvious that $`\sigma _1<\sigma _T`$ if there are collisions that can produce more than one pair. We are now in a position to discuss in detail the existing literature concerning the discrepancy between the two approaches.

#### Papers

It seems that those papers overlooked the results of . The latter papers, as well as the above derivation, only connect squares of retarded amplitudes to average numbers of particles (or other moments of the multiplicity distribution), but cannot give the pair production amplitude in terms of retarded solutions of the Dirac equation. The present derivation shows that  are not correct when they call $`\overline{u}(q)𝒯_R(q,p)v(p)`$ the pair production amplitude. The only use of this object is via its square in Eq. (50), where it leads to $`\overline{n}`$.

#### Papers

Eichmann and collaborators suggested in  that the discrepancy between the two approaches was a consequence of the fact that crossing symmetry is not valid in the ultra-relativistic limit. Here is how their argument goes: in , they claim that solving the Dirac equation in the background field of the two nuclei leads to the exact scattering amplitude of a lepton by the nuclei. From there, one would have to use crossing symmetry (i.e. change $`pp`$ for the incoming electron) in order to obtain the pair production amplitude, a procedure which they later claimed to be incorrect because crossing symmetry does not work when the nuclei are ultra-relativistic .

However, the premises of this explanation are incorrect, because solving the Dirac equation with retarded boundary conditions does not give the exact scattering amplitude in a relativistic theory if the background potential is time-dependent. Indeed, in a relativistic field theory the free retarded and Feynman propagators differ by their $`iϵ`$ prescription:

$$G_R^0(p)=i\frac{/p+m}{p^2m^2+ip_0ϵ},G_F^0(p)=i\frac{/p+m}{p^2m^2+iϵ}.$$
(53)
The two propagators will lead to equivalent results only if the background potential cannot change the sign of the energy $`p_0`$ (the energy of the incoming lepton is of course positive), i.e. only if the background potential is static. But if the background potential has some time dependence, it can change the sign of the energy carried by the propagator, and the retarded and Feynman propagators will lead to different results. Basic requirements of any field theory (Lorentz covariance, unitarity and causality) imply (see , pages 197-220) that scattering amplitudes are given by the time-ordered propagator, and not by the retarded one. Physically, the Feynman propagator takes into account the fact that an external field may create or annihilate pairs of leptons , an effect which affects lepton scattering. This effect is not taken into account if one uses the retarded propagator141414See also the section 5 of the present paper, which highlights the differences between the retarded and Feynman propagators.. Of course, if the external field is static, or if the problem is treated non-relativistically151515 The reason why scattering in non-relativistic quantum theory can be studied by solving the wave equation (even if the potential is non-static) is due to the fact that the Feynman and retarded free propagators are the same in a non-relativistic theory:
$$G_R^0{}_{}{}^{\mathrm{n}.\mathrm{r}.}(p)=G_F^0{}_{}{}^{\mathrm{n}.\mathrm{r}.}(p)=\frac{i}{p_0p^2/2m+iϵ},$$
(54)  because the propagator associated to Schrödinger’s equation has a single pole, at a positive energy., pair creation/annihilation is not possible, and both retarded and Feynman propagators give the same scattering amplitude.

The authors of  also proposed an alternate proof in section III for their expression of the scattering amplitude, this time by a direct calculation of Feynman diagrams. However, the diagrams displayed in the figure 3 of reference  do not vanish in general if calculated with Feynman propagators. Although it happens that the diagram with the configuration $`ABA`$ (see reference  for the notations) of external potentials vanishes, the diagram $`ABBA`$ for instance does not.

In conclusion, solving exactly the Dirac equation does not give the exact lepton scattering amplitude, which makes the discussion of crossing symmetry irrelevant for the present problem. Indeed, solving the Dirac equation (or, equivalently, calculating the retarded propagator), gives directly $`\overline{n}`$ thanks to Eq. (50). Of course, since one needs $`𝒯_R(q,p)`$ and not $`𝒯_R(q,p)`$, one should be careful not to assume that the incoming energy is positive when calculating $`𝒯_R`$. In the next two sections of the present paper, we discuss in more detail the differences between the retarded and Feynman propagators in the background field of the two nuclei. In particular, a closed form expression for $`G_R(q,p)`$ is derived (which is valid for energies of any sign, and therefore can be used in Eq. (50)). On the other hand, we show that it is not possible to obtain such a simple expression for the Feynman propagator.

#### Paper

In their paper , Lee and Milstein proposed a different explanation for the difference between the two approaches. Starting from a formula equivalent to our Eq. (51), and manipulating possibly ill-defined integrals with special care, they recover the formula obtained in the references  from Feynman diagrams at lowest order in $`Z_1\alpha `$. In our language, their identity reads:

$$_{b,p,q}\left|\overline{u}(q)𝒯_F(q,p)v(p)\right|_{Z_1\alpha 1}^2=_{b,p,q}\left|\overline{u}(q)𝒯_R(q,p)v(p)\right|_{Z_1\alpha 1}^2$$
(55)
However, one should refrain from trying to give a general interpretation to this result. Indeed, this identification works only at lowest order in $`Z_1`$. However, the integrand in the Eq. (51) for $`\sigma _T`$ is known to all orders in both $`Z_1`$ and $`Z_2`$, because the Dirac equation can be solved exactly with retarded boundary conditions. On the contrary, it is impossible to calculate all the corrections in $`Z_1`$ in the perturbative calculation of . We refer the reader to the sections 5 and 7.1 of the present paper, where we show that the retarded and time-ordered amplitudes are qualitatively different when considered to all orders in $`Z_1`$ and $`Z_2`$. In particular, the retarded amplitude has only a finite number of eikonal phases, while the time-ordered one can contain arbitrarily many of them. Therefore, it should be clear that the quantities calculated in  and in  cannot be the same in general.

The reason why Eq. (55) works at lowest order in $`Z_1`$ is the following: if one considers the time-ordered amplitude only at lowest order in $`Z_1`$, it simplifies dramatically into an expression that has at most two eikonal phases. If one focuses on the eikonal phases in $`\overline{u}(q)𝒯_F(q,p)v(p)`$ (see later the Eq. (106) for $`𝒯_F`$, with $`p^\pm ,q^\pm >0`$), one can show that (see , appendix E):

$`\overline{u}(q)𝒯_F(q,p)v(p)|_{Z_1\alpha 1}`$
$`={\displaystyle \mathrm{}\left[\left[e^{ie\mathrm{\Lambda }(x_{})}1\right]+\left[e^{+ie\mathrm{\Lambda }(y_{})}1\right]+\left[e^{ie\mathrm{\Lambda }(x_{})}1\right]\left[e^{+ie\mathrm{\Lambda }(y_{})}1\right]\right]}`$
$`={\displaystyle \mathrm{}e^{+ie\mathrm{\Lambda }(y_{})}\left[e^{ie\mathrm{\Lambda }(x_{})}e^{ie\mathrm{\Lambda }(y_{})}\right]},`$ (56)
where the dots represent factors that are not needed for the argument. Doing the same thing for $`\overline{u}(q)𝒯_R(q,p)v(p)`$ (see Eq. (75), truncated at lowest order in $`Z_1`$), one would get:

$$\overline{u}(q)𝒯_R(q,p)v(p)|_{Z_1\alpha 1}=\mathrm{}\left[e^{ie\mathrm{\Lambda }(x_{})}e^{ie\mathrm{\Lambda }(y_{})}\right],$$
(57)
where the dots represent exactly the same factors as in Eq. (56). One can see that the phase factors in the retarded amplitude differ from those in the time-ordered amplitude only by a global phase. It happens that this phase drops out when one is taking the modulus squared of those quantities and integrating over the momenta of the leptons, which explains Eq. (55). Because Eq. (55) relies on these properties of eikonal phases, it cannot be true for higher orders in $`Z_1`$.

Nevertheless, the calculation of reference  is very interesting, because it indicates that the integration in Eq. (51) should be handled with great care. In particular, it seems that Coulomb corrections to $`\sigma _T`$ survive the integration over impact parameter $`b`$ in the model defined by Eqs. (1) and (2).

In view of the experimental result of  which observed an almost exact $`Z^2`$ scaling of the positron yield (or, equivalently, the number of produced pairs) apparently incompatible with large Coulomb corrections, one should mention a possible shortcoming of such a model. As explained in reference , the Coulomb corrections in the integral Eq. (51) come entirely from the point $`k_{}=0`$. Therefore, they survive in $`\sigma _T`$ only if there is no cutoff preventing zero momentum transfers (either intrinsic to a theoretical calculation performed with a finite $`\gamma `$, or coming from the experimental setup).

## 4 Propagators in the presence of one nucleus

Even if the seeming discrepancy of the two approaches is now explained, there is still a paradox remaining. The work of  indicates that $`\sigma _T`$ (or, equivalently, the retarded propagator) can be calculated exactly in the case of the collision of two ultra-relativistic nuclei. On the contrary, the perturbative expansion of  leaves little hope that $`\sigma _1`$ (or the Feynman propagator) can be calculated exactly. Why are these two propagators, calculated in the same background field, so different? This is the question we address in the next two sections, by showing that in the ultra-relativistic limit there is a simplification that allows to express in closed form the retarded propagator, but does not help to calculate the Feynman propagator.

Let us first consider the case of the background field of only one nucleus of charge $`Z`$. It happens that for this case, the perturbative expansion of the scattering matrix can be summed into a closed expression when the nucleus is ultra-relativistic. This result will later appear as a building block in the expression of the propagator in the case of two nuclei. For a nucleus moving at the speed of light in the positive $`z`$ direction, the Fourier transform of the potential has the generic form $`A^\mu (k)=\delta (k^{})v_+^\mu \mathrm{\Lambda }(k_{})`$. For a point-like Coulomb interaction, $`\mathrm{\Lambda }(k_{})Zek_{}^{}{}_{}{}^{2}`$, but the following discussion does not depend on a specific form for $`\mathrm{\Lambda }(k_{})`$. Therefore, one could possibly take into account effects like an electromagnetic form factor for the nuclei.

Instead of the propagator itself, it is simpler to deal with the scattering matrix $`T(q,p)`$ (we reserve the calligraphic letter $`𝒯`$ for scattering matrices in the presence of two nuclei) defined by

$$G(q,p)(2\pi )^4\delta (pq)G^0(p)+G^0(q)T(q,p)G^0(p).$$
(58)
The term of order $`n`$ in the perturbative expansion of this object is given by

$`T_n(q,p)=(ie)^n{\displaystyle \frac{d^4k_1}{(2\pi )^4}2\pi \delta (k_1^{})\mathrm{\Lambda }(k_1)\mathrm{}\frac{d^4k_n}{(2\pi )^4}2\pi \delta (k_n^{})\mathrm{\Lambda }(k_n)}`$
$`\times /v_+G^0(p+k_1)/v_+\mathrm{}G^0(p+k_1+\mathrm{}+k_{n1})/v_+(2\pi )^4\delta (p+k_1+\mathrm{}+k_nq).`$
At this stage, one can use the $`\delta (k_i^{})`$ to perform for free all the integrations over the $`k_i^{}`$ components. Using $`/v_+(/p+m)/v_+=2p^{}/v_+`$, we obtain

$`T_n(q,p)=(ie)^n2\pi \delta (p^{}q^{})/v_+`$
$`\times {\displaystyle }{\displaystyle \frac{d^2k_1}{(2\pi )^2}}\mathrm{\Lambda }(k_1)\mathrm{}{\displaystyle }{\displaystyle \frac{d^2k_n}{(2\pi )^2}}\mathrm{\Lambda }(k_n)(2\pi )^2\delta (p_{}+k_1+\mathrm{}+k_nq_{})`$
$`\times {\displaystyle }{\displaystyle \frac{dk_1^+}{2\pi }}\mathrm{}{\displaystyle }{\displaystyle \frac{dk_n^+}{2\pi }}2\pi \delta (p^++k_1^++\mathrm{}+k_n^+q^+)`$
$`\times {\displaystyle \frac{i}{p^++k_1^+\frac{\omega _1^2}{2p^{}}+iϵ_p}}\mathrm{}{\displaystyle \frac{i}{p^++k_1^++\mathrm{}+k_{n1}^+\frac{\omega _{n1}^2}{2p^{}}+iϵ_p}},`$ (60)
where $`ϵ_p`$ is $`ϵ>0`$ in the case of the retarded propagator161616For an advanced propagator, $`ϵ_p`$ would be $`ϵ<0`$., and $`\mathrm{sign}(p^{})ϵ`$ in the case of the Feynman propagator, and where we denote:

$$\omega _a^2m^2+(p_{}+k_1+\mathrm{}+k_a)^2.$$
(61)
The next step is to perform the integrals over the $`k_i^+`$. A convenient trick is to introduce new variables $`A_ik_i^++(\omega _{i1}^2\omega _i^2)/2p^{}`$ (and $`\omega _0^22p^+p^{}`$) and write this integral as

$`I_n^+{\displaystyle \frac{1}{n!}}{\displaystyle \underset{\mathrm{perms}.\mathrm{of}\mathrm{the}A_i}{}}{\displaystyle \frac{dA_1}{2\pi }\mathrm{}\frac{dA_n}{2\pi }2\pi \delta (A_1+\mathrm{}+A_n+\frac{\omega _n^2}{2p^{}}q^+)}`$
$`\times {\displaystyle \frac{i}{A_1+iϵ_p}}\mathrm{}{\displaystyle \frac{i}{A_1+\mathrm{}+A_{n1}+iϵ_p}}.`$ (62)
Note that the $`1/n!`$ is exactly compensated by the sum over permutations of the $`A_i`$ (because these are interchangeable integration variables); it has been inserted for later convenience. Then, the following combinatoric formula

$$\underset{\sigma 𝔖_n}{}\frac{i}{A_{\sigma (1)}}\mathrm{}\frac{i}{A_{\sigma (1)}+\mathrm{}+A_{\sigma (n1)}}=\frac{i}{A_1}\mathrm{}\frac{i}{A_n}\frac{A_1+\mathrm{}+A_n}{i},$$
(63)
where $`𝔖_n`$ is the permutation group of $`[1,\mathrm{},n]`$, makes the various integrations almost independent (they are now coupled only by the $`\delta `$ function). At this stage, one begins with the $`\delta `$ function to get rid of $`A_n`$, and then performs successively the integrals over $`A_{n1}\mathrm{}A_1`$ in the complex plane. The final answer for $`I_n^+`$ is extremely simple:

$$I_n^+=\frac{(\mathrm{sign}(ϵ_p))^{n1}}{n!}.$$
(64)
The transverse integral (second line of Eq. (60)) factorizes completely in the space of transverse coordinates, so that we obtain

$$I_n^{}=d^2x_{}[\mathrm{\Lambda }(x_{})]^ne^{i(q_{}p_{})x_{}},$$
(65)
where $`\mathrm{\Lambda }(x_{})`$ is the inverse Fourier transform of $`\mathrm{\Lambda }(k_{})`$.

Collecting all the pieces, we have

$$T_n(q,p)=2\pi \delta (p^{}q^{})/v_+\mathrm{sign}(ϵ_p)d^2x_{}\frac{[ie\mathrm{sign}(ϵ_p)\mathrm{\Lambda }(x_{})]^n}{n!}e^{i(q_{}p_{})x_{}},$$
(66)
and summing over $`n`$ from $`1`$ to $`+\mathrm{}`$ to get the full $`T(q,p)`$ is now trivial171717For the scattering matrix associated with a nucleus moving in the $`z`$ direction, replace $`p^{},q^{},v_+`$ by $`p^+,q^+,v_{}`$, and use the appropriate $`\mathrm{\Lambda }(x_{})`$.:

$$T(q,p)=2\pi \delta (p^{}q^{})/v_+\mathrm{sign}(ϵ_p)d^2x_{}\left[e^{ie\mathrm{sign}(ϵ_p)\mathrm{\Lambda }(x_{})}1\right]e^{i(q_{}p_{})x_{}}.$$
(67)
We can now be more specific, and write explicitly the scattering matrix for the retarded ($`ϵ_p>0`$) and for the Feynman ($`\mathrm{sign}(ϵ_p)=\mathrm{sign}(p^{})`$) propagators181818For the advanced prescription, the result is:
$$T_A(q,p)=2\pi \delta (p^{}q^{})/v_+d^2x_{}\left[e^{ie\mathrm{\Lambda }(x_{})}1\right]e^{i(q_{}p_{})x_{}}.$$
(68)  :

$$T_R(q,p)=2\pi \delta (p^{}q^{})/v_+d^2x_{}\left[e^{ie\mathrm{\Lambda }(x_{})}1\right]e^{i(q_{}p_{})x_{}}$$
(69)
$$T_F(q,p)=2\pi \delta (p^{}q^{})/v_+\mathrm{sign}(p^{})d^2x_{}\left[e^{ie\mathrm{sign}(p^{})\mathrm{\Lambda }(x_{})}1\right]e^{i(q_{}p_{})x_{}}.$$
(70)
We observe that the retarded and Feynman results differ only by a $`\mathrm{sign}(p^{})`$ appearing in two places. This sign will turn out to be essential when we go to the case of two nuclei, basically because the interaction with the second nuclei can change the sign of the $`p^{}`$ of the electron191919One can already notice here that it is the retarded version of the scattering which appears in the solution of the Dirac solution . This should not be a surprise, since solving the Dirac equation with initial boundary conditions involves naturally the retarded propagator.  noticed the $`\mathrm{sign}(p^{})`$, but attributed to some mistake its absence in the solution of the Dirac equation. The above considerations show that the solution of the Dirac equation is correct, since it is $`T_R`$ and not $`T_F`$ that should appear in the solution with this type of boundary condition..

Using both Eq. (58) and the above results for the scattering matrix $`T`$, the interpretation of the result is rather straightforward: the propagator is the sum of two terms; one of them is the free propagator $`G^0`$ (nothing happens to the electron), and the second term contains the scattering matrix sandwiched between two free propagators. In space-time, the support of $`T`$ is on the hyper-plane $`t=z`$, where the field of the nucleus is non-zero, and we can represent the term $`G^0TG^0`$ by the diagram of figure 2.

The formulae of Eqs. (69) and (70) also illustrate the discussion of section 3.3. Indeed, they show that the retarded and Feynman prescription lead to the same result in the field of a single nucleus (i.e. in a background field that can be made static by a change of frame) if $`p_0>0`$. This is perfectly consistent with the fact that the difference only shows up in a time-dependent background field.

## 5 Propagators in the presence of two nuclei

### 5.1 Formal derivation: Watson’s series

Now that we have an exact result in the case of one nucleus, we present a formal solution for the case of two nuclei that uses the previously obtained scattering matrices $`T(q,p)`$ as its building blocks. This approach allows us to derive some results regarding the retarded and Feynman propagators in presence of two nuclei, that can be checked directly from perturbation theory.

Let us assume that we have to solve some generic Lippmann-Schwinger equation $`G=G^0+G^0VG`$ (in this section, we use very compact notations for the sake of brevity; the previous formula is in fact an integral equation), and that the potential $`V`$ receives contributions from a number of different scattering centers:

$$V=\underset{\alpha }{}V_\alpha .$$

Let us assume also that the scattering matrices $`T_\alpha `$ for the individual scattering centers are known. They satisfy

$$T_\alpha =V_\alpha +V_\alpha G^0T_\alpha .$$

Then, the full propagator $`G`$ resulting from the action of all the scattering centers can be formally written in terms of the $`T_\alpha `$. Indeed, it is a pure matter of algebra to check that the following object

$$G=G^0+\underset{\alpha }{}G^0T_\alpha G_\alpha $$
(71)
is a solution of the full Lippmann-Schwinger equation (see , pages 750-752), provided the $`G_\alpha `$ satisfy

$$G_\alpha =G^0+\underset{\beta \alpha }{}G^0T_\beta G_\beta .$$
(72)
This formal solution amounts to a reorganization of the initial perturbative expansion, which resums infinite subsets of terms corresponding to the $`T_\alpha `$. Despite this achievement, the problem of finding $`G`$ in closed form is far from being solved, because the equations (72) are coupled integral equations that are very difficult to solve. The expansion in powers of the $`T_\alpha `$ that emerges naturally from them is known in the literature as Watson’s series.

### 5.2 The case of moving scattering centers

The set (72) of coupled equations has a very intuitive interpretation in the context of wave propagation. Indeed, the “Lippmann-Schwinger” equation $`\psi =\varphi +G^0V\psi `$ ($`\varphi `$ being the incoming wave, $`G^0`$ being a propagator for the free wave equation) is solved exactly in the same way by

$`\psi =\varphi +{\displaystyle \underset{\alpha }{}}G^0T_\alpha \psi _\alpha ,`$
$`\psi _\alpha =\varphi +{\displaystyle \underset{\beta \alpha }{}}G^0T_\beta \psi _\beta .`$ (73)
In this solution, one can interpret $`G^0T_\alpha \psi _\alpha `$ as the partial wave scattered by the scattering center $`\alpha `$, and therefore $`\psi _\alpha `$ is the wave seen by the scatterer $`\alpha `$. The second equation then tells that the wave seen by the center $`\alpha `$ is made of the initial wave, and of contributions coming from the waves scattered by all the other centers. This interpretation is illustrated in figure 3 in the case of two centers.

On this classical example, it is also intuitive that multiple reflections of the wave on the two centers cannot happen if the centers are moving away at a speed larger than the velocity of wave propagation.

However, the latter remark is correct only in the case of a non-relativistic wave equation. In the case of a relativistic wave equation, the limit where the two centers are moving at the speed of light is more intricate. In the next subsection, we show that in the case of the Dirac equation, a similar conclusion holds for the retarded propagator, but not for the Feynman propagator.

### 5.3 Retarded vs. Feynman propagators

#### 5.3.1 Retarded propagator

If $`G^0`$ is the free retarded propagator, a typical contribution to the full retarded propagator is represented in the figure 4 (left) when the two nuclei move at a speed $`v`$ smaller than the velocity of light. Only this type of contribution can contribute to the retarded propagator, because the free retarded propagator connecting each scattering matrix is vanishing outside of the forward light-cone. This property has also important consequences in the case where the two nuclei are flying away at the speed of light, because the separation of two points lying respectively on the hyper-planes $`z=\pm t`$ is space-like, except if those points have $`z`$’s of the same sign. This restriction forbids terms with more than two insertions of the scattering matrices (see figure 4, right).

In the perturbative expansion of the retarded propagator, this simplification of the ultra-relativistic limit appears in the following way: terms where interactions with the second nuclei are inserted between interactions with the first nuclei, like in

$$T_{1R}G_R^0T_{2R}G_R^0T_{1R},$$
(74)
are vanishing ($`T_{1R}`$ and $`T_{2R}`$ are the retarded scattering matrices associated to the nuclei $`1`$ and $`2`$ respectively, in the retarded prescription (given by Eq. (69))). Indeed, it is immediate to check that such a term would have all its poles on the same side of the real-energy axis. More generally, the only terms allowed do not alternate interactions with the two nuclei: they can at most contain one “packet” of interactions with one nucleus, followed by a “packet” of interactions with the other nucleus. Using the above symbolic notations (integrations over the momenta exchanged with the nuclei are implicit), the full retarded propagator in presence of two ultra-relativistic nuclei reads:

$`G_R=`$ $`G_R^0`$ (75)
$`+G_R^0T_{1R}G_R^0+G_R^0T_{2R}G_R^0`$
$`+G_R^0T_{1R}G_R^0T_{2R}G_R^0+G_R^0T_{2R}G_R^0T_{1R}G_R^0.`$
These are the first three orders in the expansion of Eq. (72). All the following terms vanish because they involve factors like Eq. (74). This result is precisely the object that appeared in the solution of the Dirac equation , confirming the fact that this approach in fact derived the retarded propagator. That was to be expected given the boundary conditions used to solve the Dirac equation.

#### 5.3.2 Feynman propagator

Unlike the free retarded propagator, the free Feynman propagator can connect any pair of points in space-time. Therefore, the ultra-relativistic limit does not forbid any contribution to the Feynman propagator. Physically, the additional terms correspond to the creation of additional $`e^+e^{}`$ pairs, which are annihilated later so that only one electron is present when $`t+\mathrm{}`$. This is illustrated in figure 5.

In perturbation theory, the main difference with the retarded case is that the analogue of Eq. (74):

$$T_{1F}G_F^0T_{2F}G_F^0T_{1F},$$
(76)
is not vanishing. This can be traced back into the fact that the free Feynman propagator has poles on both sides of the real-energy axis, and this in turn is related to the fact that the free Feynman propagator can connect points with any time-ordering. As a consequence, the Watson’s series for the Feynman propagator is infinite even in the ultra-relativistic limit:

$`G_F=`$ $`G_F^0`$ (77)
$`+G_F^0T_{1F}G_F^0+G_F^0T_{2F}G_F^0`$
$`+G_F^0T_{1F}G_F^0T_{2F}G_F^0+G_F^0T_{2F}G_F^0T_{1F}G_F^0`$
$`+G_F^0T_{1F}G_F^0T_{2F}G_F^0T_{1F}G_F^0+G_F^0T_{2F}G_F^0T_{1F}G_F^0T_{2F}G_F^0`$
$`+\mathrm{}.`$
One can note that the difference between the retarded and Feynman propagators is a feature specific to the relativistic nature of the Dirac equation. Indeed, for a non-relativistic wave equation like Schrödinger’s equation, the retarded and Feynman $`iϵ`$ prescriptions give the same propagator because the propagator has a unique pole, which has a positive energy (see footnote 15). This remark is in agreement with the fact that the surviving terms involve the creation/annihilation of additional pairs, a purely relativistic effect. This interpretation in terms of pair creation/annihilation also highlights why the retarded and Feynman propagators are equally simple in the case of one nucleus, whereas they are not in the case of two nuclei: this is due to the fact that pairs cannot be produced by a single nucleus.

## 6 Unitarity

### 6.1 Calculation of $`|0_{\mathrm{out}}|0_{\mathrm{in}}|^2`$

So far, we have said nothing about the factor $`0_{\mathrm{out}}|0_{\mathrm{in}}`$ which appeared in the expression of the pair production amplitude, besides the fact that this factor is not just a phase. In this paragraph, we derive an expression for this vacuum-vacuum amplitude, that depends only on the Feynman propagator. In order to calculate $`|0_{\mathrm{out}}|0_{\mathrm{in}}|^2`$, we start from the well known fact that the corresponding amplitude is the exponential of the sum of vacuum-vacuum diagrams :

$$0_{\mathrm{out}}|0_{\mathrm{in}}=e^{iV},$$
(78)
with

$$iV\text{},$$
(79)
where the coefficients $`1/2`$, $`1/4,\mathrm{}`$ are the symmetry factors of the the corresponding diagrams. Then, the probability of vacuum to vacuum transition is given by

$$|0_{\mathrm{out}}|0_{\mathrm{in}}|^2=e^{2\mathrm{Im}V}.$$
(80)
To proceed, one could use cutting rules202020In order to use cutting rules, it is essential to notice that, since the sign of the energy flowing in the photon lines is not fixed, one cannot exclude contributions where the cut divides the diagram in more than two connected pieces. This is precisely what happens in terms like the second diagram of figure 6. in order to compute the imaginary part of the vacuum-vacuum diagrams. However, this approach makes cumbersome the tracking of symmetry factors. For this reason, it is much simpler to remark that $`V`$ is also given by212121This formula, together with Eq. (78), is well known for the vacuum-to-vacuum amplitude. In the literature on pair production in nuclear collisions, it already appears in  (see their Eq. (45)). However,  does not work out its modulus squared $`|0_{\mathrm{out}}|0_{\mathrm{in}}|^2`$.

$$iV=\mathrm{Tr}\mathrm{ln}(1+ie/AG_F^0),$$
(81)
where the symbol $`\mathrm{Tr}`$ denotes a trace on Dirac’s indices, as well as a trace over space-time. It is trivial to check that the expansion of the logarithm generates the series of diagrams in $`V`$, with the correct symmetry factors and a minus sign for the fermion loop. To begin with, one must write the Lippmann-Schwinger equation for the full Feynman scattering matrix222222Let us recall that $`G_F=G_F^0+G_F^0𝒯_FG_F^0`$. (in presence of the two nuclei). Formally, it reads

$$𝒯_F=ie/Aie/AG_F^0𝒯_F=ie/A𝒯_FG_F^0ie/A,$$
(82)
from which we deduce232323$`𝒯_F^{}\gamma ^0𝒯_F^{}\gamma ^0`$, $`G_F^0{}_{}{}^{}\gamma ^0G_F^0{}_{}{}^{}\gamma _{}^{0}`$, $`/A`$.

$$𝒯_F+𝒯_F^{}=𝒯_F^{}(G_F^0+G_F^0{}_{}{}^{})𝒯_F=𝒯_F^{}(\rho ^{(+)}+\rho ^{()})𝒯_F,$$
(83)
where we denote $`\rho ^{(\pm )}(p)2\pi \theta (\pm p_0)(/p+m)\delta (p^2m^2)`$.

Noticing then that closed loops of retarded propagators are zero, and that $`G_R^0=G_F^0\rho ^{()}`$, we have

$`0`$ $`=\mathrm{Tr}\mathrm{ln}(1+ie/AG_R^0)`$ (84)
$`=\mathrm{Tr}\mathrm{ln}(1+ie/AG_F^0ie/A\rho ^{()})`$
$`=\mathrm{Tr}\mathrm{ln}(1+ie/AG_F^0+(1+ie/AG_F^0)𝒯_F\rho ^{()})`$
$`=iV+\mathrm{Tr}\mathrm{ln}(1+𝒯_F\rho ^{()}).`$
From there, one obtains

$`2\mathrm{I}\mathrm{m}V=i(VV^{})`$ $`=\mathrm{Tr}\mathrm{ln}((1+𝒯_F^{}\rho ^{()})(1+𝒯_F\rho ^{()}))`$ (85)
$`=\mathrm{Tr}\mathrm{ln}((1+𝒯_F^{}\rho ^{()}𝒯_F\rho ^{()}+(𝒯_F+𝒯_F^{})\rho ^{()}))`$
$`=\mathrm{Tr}\mathrm{ln}(1𝒯_F^{}\rho ^{(+)}𝒯_F\rho ^{()}).`$
Therefore, we have

$$|0_{\mathrm{out}}|0_{\mathrm{in}}|^2=e^{\mathrm{Tr}\mathrm{ln}(1𝒯_F^{}\rho ^{(+)}𝒯_F\rho ^{()})}.$$
(86)
From this formula, we know $`0_{\mathrm{out}}|0_{\mathrm{in}}`$, up to an irrelevant phase. This expression shows that the vacuum-vacuum amplitude is completely determined once we have set up some approximation scheme that gives the Feynman propagator. The very reason for this is the nature of the Lagrangian of Eq. (2). Indeed, this Lagrangian implies that Wick’s theorem holds242424This would not be true if we had kept the photon kinetic term $`F_{\mu \nu }F^{\mu \nu }/4`$ in the Lagrangian., so that we can express the full $`2n`$-point Green’s functions in terms of $`G_F`$ only.

### 6.2 Probability to produce $`n`$ pairs

Since there is some confusion in the literature regarding unitarity in ultra-relativistic heavy ion collisions , it is useful to derive the probability of production of $`n`$ pairs. The corresponding transition amplitude is

$$M_n(\{p_i,q_i\})e^+(p_1)\mathrm{}e^+(p_n)e^{}(q_1)\mathrm{}e^{}(q_n){}_{\mathrm{out}}{}^{}|0_{\mathrm{in}}.$$
(87)
This amplitude can be related by a reduction formula similar to Eq. (6) to the Fourier transform of the amputated $`2n`$-point time-ordered correlator (with appropriate spinors for the final states). Again, we pull out the factor $`0_{\mathrm{out}}|0_{\mathrm{in}}`$, in order to get rid of vacuum-vacuum diagrams in the second factor:

$$0_{\mathrm{out}}|0_{\mathrm{in}}\frac{0_{\mathrm{out}}|\mathrm{T}\overline{\psi }(y_1)\mathrm{}\overline{\psi }(y_n)\psi (x_1)\mathrm{}\psi (x_n)|0_{\mathrm{in}}}{0_{\mathrm{out}}|0_{\mathrm{in}}}0_{\mathrm{out}}|0_{\mathrm{in}}G_{2n}(\{x_i,y_i\}).$$
(88)
At this stage, it is important to notice that since the Lagrangian Eq. (2) does not contain any dynamical field that couples to the fermions, Wick’s theorem applies to the second factor $`G_{2n}(\{x_i,y_i\})`$ so that we can write it in terms of the full Feynman propagator:

$$G_{2n}(\{x_i,y_i\})=\underset{\sigma 𝔖_n}{}ϵ(\sigma )G_F(x_1,y_{\sigma (1)})\mathrm{}G_F(x_n,y_{\sigma (n)}),$$
(89)
where $`ϵ(\sigma )`$ is the signature of the permutation $`\sigma `$ (required when permuting fermion fields). This property extends to the amputated correlators and their Fourier transform, so that we can write directly the transition amplitude as252525This formula is equivalent to the Eq. (11) of . However, these authors do not make use of it to calculate the exact $`P_n`$ (see our Eq. (95)).

$`M_n(\{p_i,q_i\})=0_{\mathrm{out}}|0_{\mathrm{in}}{\displaystyle \underset{\sigma 𝔖_n}{}}ϵ(\sigma )[\overline{u}(q_{\sigma (1)})𝒯_F(q_{\sigma (1)},p_1)v(p_1)]\mathrm{}`$
$`\mathrm{}[\overline{u}(q_{\sigma (n)})𝒯_F(q_{\sigma (n)},p_n)v(p_n)],`$ (90)
where $`𝒯_F`$ is the Feynman scattering matrix in presence of the two nuclei.

The integrated probability of producing exactly $`n`$ pairs is then (we have changed $`p_ip_i`$ in the second line)

$`P_n{\displaystyle \frac{1}{n!^2}}{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{d^4p_i}{(2\pi )^4}\frac{d^4q_i}{(2\pi )^4}\stackrel{~}{\rho }^{(+)}(p_i)\stackrel{~}{\rho }^{(+)}(q_i)\left|M_n(\{p_i,q_i\})\right|^2}`$
$`={\displaystyle \frac{|0_{\mathrm{out}}|0_{\mathrm{in}}|^2}{n!}}{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{d^4p_i}{(2\pi )^4}\frac{d^4q_i}{(2\pi )^4}\stackrel{~}{\rho }^{()}(p_i)\stackrel{~}{\rho }^{(+)}(q_i)}`$
$`\times {\displaystyle \underset{\sigma 𝔖_n}{}}ϵ(\sigma )(1)^n{\displaystyle \underset{i=1}{\overset{n}{}}}\left[\overline{u}(q_i)𝒯_F(q_i,p_i)(/p_i+m)𝒯_F^{}(p_i,q_{\sigma (i)})u(q_{\sigma (i)})\right],`$ (91)
where $`\stackrel{~}{\rho }^{(\pm )}(p)2\pi \theta (\pm p_0)\delta (p^2m^2)`$. We can see now that we obtain closed chains like:

$$\mathrm{tr}[𝒯_F(q_1,p_1)(/p_1+m)𝒯_F^{}(p_1,q_2)(/q_2+m)𝒯_F(q_2,p_2)(/p_2+m)𝒯_F^{}(p_2,q_1)(/q_1+m)].$$
(92)
For instance, for $`n=2`$, cycles come in two sizes, illustrated on figure 6. A systematic tool to construct these loops is the (unique) decomposition of permutations in products of disjoint circular permutations. It is a trivial matter of combinatorics to find that the number of $`p`$-cycles in $`𝔖_n`$ is $`n!/p(np)!`$, and that the number of permutations that are made of $`a_1`$ $`1`$-cycles, $`a_2`$ $`2`$-cycles,$`\mathrm{}`$,$`a_n`$ $`n`$-cycles ($`a_1+2a_2+\mathrm{}+na_n=n`$) is:

$$\frac{n!}{a_1!\mathrm{}a_n!}\frac{1}{1^{a_1}}\mathrm{}\frac{1}{n^{a_n}}.$$
(93)
The signature of such a permutation is $`ϵ(\sigma )=(1)^n_i(1)^{a_i}`$. With the following compact notations for “links”262626Note that $`L=0`$ in a background field that cannot produce pairs, like the field of a single nucleus for instance. and “loops”:

$`L(p^{},p){\displaystyle \frac{d^4q}{(2\pi )^4}𝒯_F^{}(p^{},q)\rho ^{(+)}(q)𝒯_F(q,p)\rho ^{()}(p)},`$
$`\mathrm{Tr}L^n{\displaystyle \left[\underset{i=1}{\overset{n}{}}\frac{d^4p_i}{(2\pi )^4}\right]\mathrm{tr}[L(p_1,p_2)L(p_2,p_3)\mathrm{}L(p_n,p_1)]},`$ (94)
the total probability to produce $`n`$ pairs is

$$P_n=|0_{\mathrm{out}}|0_{\mathrm{in}}|^2\underset{a_1+2a_2+\mathrm{}+na_n=n}{}\underset{i=1}{\overset{n}{}}\left[\frac{(1)^{a_i}}{a_i!}\left(\frac{\mathrm{Tr}L^i}{i}\right)^{a_i}\right].$$
(95)
We should emphasize here the fact that the individual $`P_n`$ are functions of the Feynman scattering matrix $`𝒯_F`$, while the average number of pairs $`\overline{n}=_nnP_n`$ has a simple expression in terms of the retarded scattering matrix $`𝒯_R`$ (see Eq. (50)).

### 6.3 Another look at unitarity

We can now check unitarity without making any approximation, by using the previous formula for the $`P_n`$, and Eq. (86) for the vacuum-vacuum transition probability:

$`{\displaystyle \underset{n=0}{\overset{+\mathrm{}}{}}}P_n`$ $`=|0_{\mathrm{out}}|0_{\mathrm{in}}|^2{\displaystyle \underset{n=0}{\overset{+\mathrm{}}{}}}{\displaystyle \underset{a_1+2a_2+\mathrm{}+na_n=n}{}}{\displaystyle \underset{i=1}{\overset{n}{}}}\left[{\displaystyle \frac{(1)^{a_i}}{a_i!}}\left({\displaystyle \frac{\mathrm{Tr}L^i}{i}}\right)^{a_i}\right]`$ (96)
$`=e^{\mathrm{Tr}\mathrm{ln}(1𝒯_F^{}\rho ^{(+)}𝒯_F\rho ^{()})}e^{\mathrm{Tr}\mathrm{ln}(1L)}=1.`$
The compensation of the two exponentials is obvious from the first of Eqs. (94).

It therefore clarifies the long standing “unitarity problem” floating around in the literature related to pair production in nuclear collisions . The problem can be stated as follows: perturbation theory seems to lead to cross sections that are too large in order to comply with unitarity (or to production probabilities larger than $`1`$). It appears that the factor $`0_{\mathrm{out}}|0_{\mathrm{in}}`$ (the modulus of which is smaller than one) has been overlooked272727Two exceptions are  and , which noticed that unitarity is related to the vacuum to vacuum amplitude. However, the authors of  calculated this amplitude by requiring that unitarity is preserved (in addition, they only did that in an approximation that leads to a Poissonian distribution). One needs an independent calculation of $`0_{\mathrm{out}}|0_{\mathrm{in}}`$ in order to claim that this factor restores unitarity. The authors of  give an exact expression of $`0_{\mathrm{out}}|0_{\mathrm{in}}`$ equivalent to our Eq. (81), but unfortunately do not exploit it to check unitarity. in the literature: only the connected piece of Feynman diagrams has been considered, but not the disconnected vacuum-vacuum diagrams. In this paper, we have shown that the factor $`0_{\mathrm{out}}|0_{\mathrm{in}}`$ naturally emerges from the reduction formula for pair production (and therefore that calculating only the connected diagrams leads to an incomplete result), and that this factor restores unitarity.

Moreover, since the factor $`0_{\mathrm{out}}|0_{\mathrm{in}}`$ depends only on the Feynman propagator, the above considerations provide a way to make approximations that preserve unitarity. It is sufficient to use the same approximate $`G_F`$ in the calculation of the $`P_n`$ and in the calculation of the vacuum-vacuum amplitude.

For such an approximation scheme to be consistent, we have also to verify that it gives positive $`P_n`$. Indeed, the fact that the sum of Eq. (96) is $`1`$ leaves open the possibility that some probabilities could come out negative. To prove that all the probabilities are positive (and hence smaller than one, because their sum is one), it is convenient to introduce a “generating function” for the probabilities $`P_n`$:

$$F(x)e^{\mathrm{Tr}\mathrm{ln}(1L)}e^{\mathrm{Tr}\mathrm{ln}(1xL)},$$
(97)
such that $`F(1)=1`$ and $`P_n=F^{(n)}(0)/n!`$. Then, one notices that it can be rewritten as

$$F(x)=e^{\mathrm{Tr}\mathrm{ln}(1+tt^{})}e^{\mathrm{Tr}\mathrm{ln}(1+xtt^{})},$$
(98)
with $`t(q,p)\overline{u}(q)𝒯_F(q,p)v(p)`$. The operator $`tt^{}`$ is positive, and we can write the generating function in terms of its eigenvalues $`\tau _i`$ , which are positive:

$$F(x)=\underset{i}{}\frac{1+\tau _ix}{1+\tau _i}.$$
(99)
We therefore see that all the coefficients in the Taylor expansion of $`F(x)`$ are positive independently of the approximation made for $`𝒯_F`$, which proves that all the $`P_n`$ are positive in this framework.

### 6.4 Nature of the multiplicity distribution

One can also note that Eq. (95) for the probability of producing $`n`$ pairs is not a Poisson distribution. This fact contradicts the solution proposed in  for the unitarity problem, where a Poissonian distribution is obtained. However, one can see from Eq. (95) the nature of the approximation that would lead to such a distribution: a Poisson distribution is obtained if one drops all the $`\mathrm{Tr}L^n`$ for $`n>1`$. Indeed, this drastic (and a priori not justified) simplification leads to

$$P_ne^{(\mathrm{Tr}L)}\frac{(\mathrm{Tr}L)^n}{n!},$$
(100)
a Poisson distribution of average $`\mathrm{Tr}L`$. This explains the somewhat confusing statement of  saying that the “perturbation theory prediction” for the probability to produce one pair (i.e. $`P_1`$ in which one would forget the factor $`|0_{\mathrm{in}}|0_{\mathrm{out}}|^2`$, that is $`\mathrm{Tr}L`$ in our notations) should not be interpreted as a probability (because it can be larger than $`1`$) but instead should be reinterpreted as the average number of pairs produced in a nuclear collision. Modulo the approximation of Eq. (100), we agree with this statement, except for one thing: $`\mathrm{Tr}L`$ is not the perturbation theory prediction for the probability to produce one pair; perturbation theory applied correctly indicates that this probability is $`|0_{\mathrm{out}}|0_{\mathrm{in}}|^2\mathrm{Tr}L<1`$. In fact, the present analysis makes the following clear: even if $`𝒯_F`$ (and hence $`L`$) were known exactly, using $`\mathrm{Tr}L`$ for the probability to produce one pair would still violate unitarity. The unitarity problem does not really come from perturbation theory, but from forgetting the contribution of the vacuum-vacuum amplitude.

Thanks to the generating function introduced above, it is easy to find what the average number of produced pairs is282828One can check that this formula is equivalent to Eq. (50), by using the Lippmann-Schwinger equations for $`𝒯_F`$ and $`𝒯_R`$. One can obtain successively:    $`𝒯_F=[1𝒯_R\rho ^{()}]^1𝒯_R,`$    $`L=[1𝒯_R^{}\rho ^{()}]^1[𝒯_R^{}\rho ^{(+)}𝒯_R\rho ^{()}][1𝒯_R\rho ^{()}]^1,`$    $`[1𝒯_R^{}\rho ^{()}][1𝒯_R\rho ^{()}][𝒯_R^{}\rho ^{(+)}𝒯_R\rho ^{()}]=1,`$  (101)  and finally
$$\overline{n}=\mathrm{Tr}[L(1L)^1]=\mathrm{Tr}[𝒯_R^{}\rho ^{(+)}𝒯_R\rho ^{()}],$$
(102)  which can be cast into Eq. (50) thanks to $`/q+m=_{\mathrm{spin}}u(q)\overline{u}(q)`$ and $`/pm=_{\mathrm{spin}}v(p)\overline{v}(p)`$.:

$$\overline{n}\underset{n=0}{\overset{+\mathrm{}}{}}nP_n=F^{}(1)=\mathrm{Tr}[L(1L)^1]=\underset{n=1}{\overset{+\mathrm{}}{}}\mathrm{Tr}L^n.$$
(103)
By the same method, we can find the variance of the number of pairs:

$$\overline{n^2}\overline{n}^2=\mathrm{Tr}[L(1L)^1]\mathrm{Tr}[L^2(1L)^2].$$
(104)
Since $`\overline{n}\overline{n^2}\overline{n}^2`$, the exact probability distribution cannot be a Poisson distribution. This is a consequence of the permutations in the final state when squaring the amplitude, which can lead to assign to the same pair292929The problem comes from the fact that particles are produced in pairs here. As an example of model in which particles are produced individually, one can consider a toy model in which a scalar field $`\varphi (x)`$ is coupled to a background classical source $`j(x)`$ by $`_{\mathrm{int}}j(x)\varphi (x)`$. This model is exactly solvable, and the production probability is found to be Poissonian (see , pages 163-170). It is also possible to understand from this remark why  obtained a Poisson distribution. Indeed, this paper modeled $`e^+e^{}`$ pairs as elementary fields (quasi-bosons), so that the problem of pairings in the final state never shows up. the electron of one fermion line and the positron from another fermion line. This is illustrated in figure 6 for the case of two pairs, where one can see clearly the origin of the terms $`(\mathrm{Tr}L)^2`$ and $`\mathrm{Tr}L^2`$. One also sees that the term in $`\mathrm{Tr}L^2`$ correlates the emission of the two pairs, and therefore prevents the probability distribution from being Poissonian303030Another argument against a Poisson distribution is that it would contradict the fact that $`\overline{n}`$ can be calculated exactly, while the individual $`P_n`$ cannot.. This explains the fact that a Poisson distribution is obtained only if one neglects the $`\mathrm{Tr}L^n`$ for $`n>1`$.

The authors of  attempted to size the departure from a Poisson distribution in the Magnus model, and found it to be around $`1\%`$, which may justify for most practical purposes to calculate only $`\overline{n}`$ (which can be calculated exactly according to our analysis) and plug it in a Poisson formula for $`P_n`$. This is the approach followed in .  started from the lowest order (only one photon is exchanged with each nucleus, and $`0_{\mathrm{out}}|0_{\mathrm{in}}`$ is approximated by 1) $`P_1`$ in the background field of two nuclei, “reinterpreted” as $`\overline{n}`$.  followed the same approach by starting from the formula given in Eq. (53) of , which was assumed to be the exact formula for the pair production amplitude. If this assumption were correct, its square should be the exact probability to produce one pair, and should therefore be smaller than $`1`$. This is not the case as noted in , who “reinterpreted” this quantity as the average number of pairs produced in a collision in order to save unitarity. These authors did not realize that their remark was directly pointing to the result we have proven in our section 3: that the retarded solution of the Dirac equation does not give the pair production probability, but rather the average number of pairs (without any need to “reinterpret” anything).

## 7 Strategies of approximation for $`G_F(x,y)`$

### 7.1 Limit $`Z_1\alpha 1`$

Since the series in Eq. (77) giving the full Feynman propagator cannot be summed exactly, one must use approximations in order to simplify it. We present in this section two different approximation schemes which lead to closed form expressions for the Feynman propagator.

The simplest approximations one can think of are truncations of the Watson series. Such a simplification is obtained in a natural way if one assumes that only one nucleus has a large electric charge. If $`Z_1`$ is small enough so that we have

$$Z_1\alpha 1Z_2\alpha ,$$
(105)
then we can neglect most of the corrections in $`Z_1\alpha `$. Since we are interested in pair production, the kinematics requires at least one interaction with the first nucleus (otherwise, the pair cannot be produced on-shell). We are therefore going to keep only terms with at most one insertion of $`/A_1`$. In that approximation, $`T_{1F}`$ is just $`ie/A_1`$, and we drop any term with two or more insertions of $`T_{1F}`$ in the Watson series. It is immediate to verify that the only terms left in $`G_F`$ are

$`G_F`$ $`G_F^0`$ (106)
$`+G_F^0(ie/A_1)G_F^0+G_F^0T_{2F}G_F^0`$
$`+G_F^0(ie/A_1)G_F^0T_{2F}G_F^0+G_F^0T_{2F}G_F^0(ie/A_1)G_F^0`$
$`+G_F^0T_{2F}G_F^0(ie/A_1)G_F^0T_{2F}G_F^0.`$
This approximation for $`G_F`$ is in fact the starting point used by  in their approach to the problem of pair production. The terms that participate to pair production (i.e. having at least one interaction with each nucleus) are displayed on figure 7.

The term labeled $`(1)`$ in this figure corresponds to the creation of an on-shell electron and an off-shell positron. The positron is subsequently put on its mass-shell by an additional interaction. For the process $`(2)`$, both the electron and the positron are first created off-shell, and then interact independently to go on-shell.

This approximation, which leads to an analytically tractable expression, suffers however from several limitations. The obvious one is that for the collision of two nuclei like gold for RHIC, the charge of the nuclei is not small enough to justify this approximation, and the inequality (105) is not satisfied313131For the case $`Z_1Z_2`$, a term has to be added were the roles of $`Z_1`$ and $`Z_2`$ are exchanged. The accuracy of the result obtained via this procedure is claimed to be better than the percent level for RHIC energies .. Moreover, this approximation leaves out some terms that seem physically important. In particular, it does not include any term where both the $`e^+`$ and the $`e^{}`$ in the pair interact with the two nuclei. Such a term would indeed have four scattering matrices.

### 7.2 Limitation in the number of intermediate pairs

Another way to approximate the Feynman propagator is to start from the definition

$$G_F(x,y)\frac{0_{\mathrm{out}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}}{0_{\mathrm{out}}|0_{\mathrm{in}}}.$$
(107)
Then, one can use a complete323232Strictly speaking, these states form a basis of the subspace containing only Fock states with zero electric charge. This is sufficient here since they will be contracted with the vacuum. set of states $`\{|n_{\mathrm{in}}\}_{n=0\mathrm{}\mathrm{}}`$ ($`n`$ counts the number of pairs, other continuous indices like the momenta of the particles have not been written explicitly) and insert in the previous equation an identical operator constructed as

$$\mathit{1}=\underset{n=0}{\overset{+\mathrm{}}{}}|n_{\mathrm{in}}n_{\mathrm{in}}|.$$
(108)
Separating the contribution from $`n=0`$, we obtain

$`G_F(x,y)=`$ $`0_{\mathrm{in}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$ (109)
$`+{\displaystyle \underset{n=1}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{0_{\mathrm{out}}|n_{\mathrm{in}}}{0_{\mathrm{out}}|0_{\mathrm{in}}}}n_{\mathrm{in}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}.`$
In this formula, $`n`$ can be seen as the number of extra pairs produced (and then destroyed) in the course of the evolution of the electron. Another way to simplify the Feynman propagator is therefore to truncate the previous sum, and keep only the term obtained with $`n=0`$, that is

$$G_F(x,y)0_{\mathrm{in}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}.$$
(110)
One can note that both the exact Feynman propagator $`G_F`$ and the approximation $`0_{\mathrm{in}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$ are solutions of the Dirac equation $`(i/_xe/A(x)m)G(x,y)=\delta (xy)`$. Their difference is therefore a solution of the homogeneous Dirac equation. It is in fact immediate to verify that all the terms with $`n1`$ in the right hand side of Eq. (109) are solutions of the homogeneous Dirac equation.

The calculation of section 3.2.2 in fact gives the answer for this correlator in terms of retarded and advanced propagators only: we have an exact expression for the Fourier transform $`G_{++}(q,p)`$ of $`0_{\mathrm{in}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$, which reads

$$G_{++}(q,p)=\frac{1}{2}\left[G_R(q,p)+G_A(q,p)+G_S(q,p)\right],$$
(111)
with $`G_S(q,p)`$ given in Eq. (42). We see that contrary to the exact $`G_F(q,p)`$, this correlator can be written in closed form, in terms of the exactly known $`G_R`$ and $`G_A`$:

$`G_{++}(q,p)={\displaystyle \frac{1}{2}}[G_R(q,p)+G_A(q,p)`$
$`+{\displaystyle }{\displaystyle \frac{d^4k}{(2\pi )^4}}2\pi \delta (k^2m^2)G_R^0(q)𝒯_R(q,k)(/k+m)𝒯_A(k,p)G_A^0(p)],`$ (112)
where $`𝒯_R`$ and $`𝒯_A`$ are the retarded and advanced scattering matrices in presence of the two nuclei. Each of them can contain up to two of the 1-nucleus scattering matrices $`T_{R,A}`$ (see Eq. (75)).

A common feature of the two expressions Eqs. (42) and (111) and Eq. (112) for the approximation of $`G_F`$ by $`0_{\mathrm{in}}|\mathrm{T}\overline{\psi }(y)\psi (x)|0_{\mathrm{in}}`$ is that they contain terms involving up to four scattering matrices on individual nuclei. Therefore, this approximation contains a little more than the previous one (based on $`Z_1\alpha 1`$), since it includes a selected subset of the terms of order four in the Watson series. The physical meaning of those terms is the following: an off-shell pair is first created. Then, the electron and the positron independently scatter off each of the two nuclei. No additional pair is created or annihilated. These terms are illustrated in figure 8.

## 8 Conclusions

In this paper, we have studied various aspects of the problem of pair production in the collision of ultra-relativistic heavy ions, focusing in resolving the discrepancies of the existent literature.

By showing that the inclusive cross-section of pair production is related to retarded amplitudes, while the Feynman amplitude gives the exclusive cross-section of single pair production, we found that the discrepancy between two methods used to attack the problem of pair production lies at a deeper level than expected: these two methods do not calculate the same physical quantity, and their result should not be compared directly.

Then, we have studied the propagator of an electron in the electromagnetic background field created by two colliding nuclei. It appeared that the exact retarded propagator can be obtained in closed form if the nuclei are ultra-relativistic. However, such a simplification does not occur for the Feynman propagator, which can only be expressed as an infinite series. This problem arises in any background field that can produce pairs of particles, and can be traced back into the relativistic nature of the wave equation governing the electron field. This observation indicates that in the collision of two ultra-relativistic nuclei, the inclusive cross-section of pair production can be obtained exactly, but not the more exclusive ones. Experimentally, it would therefore be desirable to measure both types of cross-sections.

Within the model of 2, the inclusive cross-section, which is expressed in terms of the exactly known retarded amplitude, seems to contain Coulomb corrections as well. However, the calculation of completely integrated cross-sections is questionable in this model. Indeed, the problem raised by Lee and Milstein occurs at zero momentum transfer (or at infinite impact parameter), precisely where the model is expected to break down. In particular, any infrared cutoff on the momentum transfer, coming either from an improved theoretical model or from experimental cuts, could drastically reduce these corrections.

We have also analyzed the unitarity puzzle, and shown that pair production probabilities satisfy all requirements of unitarity if the factor $`|0_{\mathrm{out}}|0_{\mathrm{in}}|^2`$ is correctly taken into account. In addition, since everything depends on the 2-point Feynman Green’s function, unitarity is preserved if one starts from an approximation of this propagator. A side product of this analysis is that the multiplicity distribution is not Poissonian.

Finally, in addition to the approximation $`Z_1\alpha 1Z_2`$ used in , we have presented a completely different approximation for $`G_F(x,y)`$ that leads also to closed expressions, and seems to contain more of the relevant physics for pair production by two heavy nuclei.

Acknowledgements: This work is supported by DOE under grant DE-AC02-98CH10886. A.P. is supported in part by the A.-v.-Humboldt foundation (Feodor-Lynen program). We would also like to thank R. Venugopalan and C. Bertulani for interesting discussions. We thank U. Eichmann for pointing out the fact that Coulomb corrections may be present in the inclusive cross-section as well.