Pair production at the edge of the QED flux tube

We investigate the process of Abelian pair production in the presence of strong inhomogeneous and time-dependent external electric fields. The spatial dependence of the external field is motivated by a non-Abelian color flux tube in heavy-ion collisions. We show that the inhomogeneity significantly increase the particle yield compared to that in the commonly used models with a constant and homogeneous field. Moreover our results indicate that in contrast to the latter, most of the particles are produced at the interface of the field profile in accordance with Heisenberg's prediction.


Introduction
Recently, pair production from vacuum is gaining interest from both theorists and experimentalists. In the Abelian (QED) case it is considered to be the final frontier of high-energy laser experiments. While the attainable energy density of today experiments is still orders of magnitudes below the threshold defined by the critical field E cr = m 2 c 3 e ≈ 1.32 · 10 18 V m [1], the scale above which pair production is sizable, the development of technology is remarkably sustaining its exponential growth both of energy density and of the frequency of laser pulses. The recent proposals in this area promise to improve these parameters further in the coming years [2,3]. The high interest of research in this area is signaled by the number of high intensity laser experiments under commissioning, construction and planning, such as the Extreme Light Infrastructure (ELI). For a comprehensive list see Ref. [4].
Another motivation comes from ultrarelativistic heavyion collisions performed at Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC). While the microscopic mechanism of hadron production is still not fully understood there is continued effort to develop and improve models that explain experimental data. The investigation of one of the most precisely measured observable, the transverse momentum spectra of produced hadrons, led to the development of a family of models that center around the concept of chromoelectric flux tubes ('strings'). These tubes connect the quark and diquark constituents of colliding protons [5,6,7]. When the sources of these fluxtubes separate, the field energy increases until the threshold of pair production is reached and new quarkantiquark and diquark-antidiquark pairs are created. Such models can describe experimental data successfully at low p T , p T < 2 − 3 GeV, while at higher p T perturbative QCD-based models work well [8,9,10]. However the insight into the microscopical details is quite limited, because these models usually assume homogeneous and often static approximations of field strength [11,12], while in reality the collision and successive events take place on short timescales, and finite size effects may also play an important role.
During the past decades, kinetic description was formulated to describe pair production in arbitrary spacetime dependent fields in the Abelian case both for fermions [13] and bosons [14,15] and for the non-Abelian case for quarks [16,17,18] and gluons [19]. While remarkable analytic progresses have been made [21,20,22], calculating pair production in arbitrary space-time dependent fields is analytically unmanageable and also numerically very demanding. For this reason, usually only the homogeneous and often only the time independent cases are investigated and used [23,24,25]. Still, even in the most simple cases the interplay of field parameters on particle spectra and yield is so complex [18,26,27,28] that before building full featured simulation frameworks, the independent role of field parameters should be understood.
In the case of laser experiments, the pair production is the shortest process to take place and followed by other processes, like back-reaction, radiation reactions etc. While the laser-plasma interaction simulation packages developed today are focused on the physical reactions attainable with todays' laser energies, they expected to form the core of later full featured simulation environments, as it happened in particle and nuclear physics. In this context it is extremely important to fully understand the interplay of parameters that influence pair production observ-ables so that future facilities can be planned for such measurements based on these simulation frameworks. Also, in heavy-ion physics pair production is one of the early stage processes and followed by multi particle interaction, energy loss, thermalization and fragmentation into hadrons that may be more or less important depending on the particle energy. The common feature in these two scenarios is that it is not clear how important the secondary processes are compared to the initial particle creation.
In this paper, we use the three dimensional Dirac-Heisenberg-Wigner evolution equations for QED to model pair production in an inhomogeneous external electric field. The theoretical details are summarized in Section 2. The numerical method is outlined in Section 3. Our results for the longitudinal and transverse spectra are presented in Sections 4 and 5 and summarized in Section 6.

The Dirac-Heisenberg-Wigner formalism
Pair production is a quantum phenomena and thus needs a proper quantum description. Also, to account for extreme scales needed for these processes, a relativistic description is necessary. An insightful description in phasespace is available in the form of singe-time relativistic one-particle Wigner function formalism [13]. The Wigner function is a quantum generalization of the classical oneparticle distribution function. However, in contrast to the latter the Wigner function does not posses a probabilistic interpretation in the strict sense. The Wigner function describes the quantum interference of negative energy states in the form of negative values it may posses, but these are restricted to non-connected areas with extent of the order of few . By integrating with a gaussian envelope function that mimics a classical measurement one ends with a classically interpretable distribution function.
While we would like to describe non-Abelian pair production, it was shown already [18] that there is a strong Abelian dominance that enables us to use the U(1) formalism of QED instead of the more complicated SU(N) models.
The definition of the one particle Wigner function in terms of the wave function is: where the Wilson line factor is to guarantee that p is the proper eigenvalue of the kinetic momentum operator. The evolution equations are derived in the temporal gauge in Refs. [13,29]. In this approach, the external field is treated as classical. This approximation is justified by the strong field strength needed for pair production in laser experiments. In high-energy heavy-ion collisions the gluon density is also high [30] thus it can be approximated by its classical expectation value. The evolution equations of the Wigner function contain three non-local differential operators defined as integrals in Fourier space. By virtue of the gradient expansion, the evolution equations can be approximated by series in configuration space that include increasing gradients of external fields multiplied by increasing order of momentum derivates which act on the Wigner function. In our case, we restrict ourselves to electric fields. There is no magnetic contribution to the spatial gradients and to the momentum. We also choose a linearly polarized electric field in the z direction with a transverse (x) spatial dependence. With this electric field the time derivate operator becomes: This is an exact expression for the chosen electric field because it has no component in the direction of inhomogeneity and in the direction of the electric field there is no inhomogeneity, so the higher momentum derivates vanish. This is in contrast to [29,31] where the inhomogeneity is in the direction of the electric field and the authors show that the naive gradient approximation breaks down.
In this case, the evolution equations for the 16 real components of the Wigner function simplify to: where ∇ x × Ø i is understood as the vector (0, −∂ x Ø iz , ∂ x Ø iy ). Only two derivatives remain, namely the momentum derivatives parallel to the electric field and the spatial derivates transverse to the electric field. While it is tempting to use the method of characteristics to decouple the flow term derivates as in the homogeneous case, there is no gain in this case, because the implied momentum will be different at different coordinates and the two directions remain coupled by the spatial gradients. The initial conditions where ω 2 ( p) = m 2 +p 2 x +p 2 y +p 2 z , corresponds to the vacuum in-state. The final observable that one may consider as the 'density of pairs' is the energy density, that is a linear combination of the mass density and the current density: In the homogeneous case the DHW equations reduce to the quantum Vlasov equation with one of its components f describing the density of pairs [26]: with p 2 T = p 2 x + p 2 y , p z = p z0 − eA(t) and the initial conditions f = u = v = 0 for all p. Also, f and ǫ are related by: For the Sauter field E(t) = E 0 cosh(t/τ ) −2 the analytic form of the asymptotic pair density is known (in fact for arbitrary spin particles) [32]: where We use this result to validate our numerics and to assess the effect of inhomogeneity. We chose the following 'plateau' field, that models the cross-section of a chromoelectric flux tube in the Abelian limit, assuming homogeneous field in the middle and exponential decay at the edges: Thus we have two parameters for the spatial dependence: R the width of the plateau (or radius of the flux tube) and r describing the steepness of the gradient at the edge. In the time direction we use the Sauter field so we can expect to reproduce the homogeneous result at x ≈ 0 when R ≫ r. We restrict ourselves to a one-dimensional inhomogeneity in three dimensions. Practically this field is an infinite homogeneous flux plane possessing a finite extent and gradient only in the direction perpendicular to this plane.

Numerical method
We solve the equation system (3-10) by explicit finite difference integration in the time direction with a usual 8th order Runge-Kutta stepper and account for the spatial derivates with pseudo-spectral collocation over the rational Chebyshev polynomial basis [33]. These polynomials resolve doubly infinite ranges with a user defined scale. The reason for the choice of this method is that the Wigner function is free of non-analyticities and for this class of functions the spectral method has superior convergence rates over finite difference techniques. Another advantage is that the collocation points coincide with those optimal for the integration quadrature . Thus the integrations can be carried out with no further efforts in spectral space with a Clenshaw-Curtis quadrature modified for the rational Chebyshev basis set [34].
Within the pseudo-spectral method, our differential operators turn into dense matrices. Because the operators acting on the Wigner-function components are time independent, we can solve the equations in spectral space, the back and forward transformation in each time step is not necessary. We use this method also because it can be easily extended and programmed on graphical processors (GPUs).
The free scale parameters of the rational Chebyshev polynomials should be estimated before calculations. A good choice for the x direction is in the order of 2R + r, while for the scale of the longitudinal momentum p z the integral of the electric field sets the scale:

Longitudinal spectra
At first we restrict ourselves to zero transverse momentum. It is known that the p T only acts as an additional mass term regarding the particle yield; increasing it results in an exponential decay of the pair production.
The 3-d plots on Figs. 1. and 2. we show the asymptotic (t → ∞) pair density f , as defined by Eq. (17), as a function of the transverse coordinate x and longitudinal (z) momentum for two different values τ = 0.3λ c /c and τ = 2λ c /c. We note that the first value corresponds to the local maxima of the homogeneous Sauter model in τ space. The middle of the plateu behaves like the homogeneous solution, decreasing with increasing τ , but the inhomogeneous edge tends to increase the pair density in its vicinity proportional to the pulse width. Next to the edge a negative region develops that signifies the admixture of antiparticles. The less steeper the gradient, the less negative density appears. We note that f is not a component of the Wigner function, therefore there are no bounds on the values it may take.
To have a scalar quantity to compare to the homogeneous case we calculate the particle yield integrated in the transverse coordinate and longitudinal momentum: Pair density (f) at p T =0 Figure 1: Phase space view of the asymptotic pair density with parameters: E 0 = 0.5Ecr, R = 5λc, r = λc, τ = 0.3λc/c. Note that the E field has the steepest gradient at z = 5λc. Pair density (f) at p T =0 Figure 2: Phase space view of the asymptotic pair density with parameters: E 0 = 0.5Ecr, R = 5λc, r = λc, τ = 2λc/c. Note that the E field has the steepest gradient at x = 5λc.  This quantity can be compared to the homogeneous case by calculating it with the same E(x, t) field as in the inhomogeneous one and setting p T = 0 and performing the same integrals in x and p z . First we plot n L as a function of the Sauter pulse width τ , while keeping R and r fixed on Fig. 3. We see that the homogeneous and inhomogeneous results for small temporal widths coincide as expected as there is no time for the particles to be created by the inhomogeneity. The two curves reach a local maximum as in the homogeneous case, but they start to separate. For large τ both curves are expected to be proportional to it, and finally will be approximated with the form of ∝ τ × n static , where n static is the constant static solution (for the homogeneous field it is proportional to the Schwinger formula). We find, that the onset of this approximation happens earlier for the inhomogeneous configuration and thus the observed particle yield is more than a magnitude larger than in the homogeneous case.
If we fix pulse width τ = 2λ c /c and gradient r = λ c and vary R, we again expect a linear proportionality since we are changing the interaction volume. We indeed find a linear relation (Fig. 4) and also the slopes turn out to be the same within 5%. This means that the homogeneous and inhomogeneous solutions are almost identical in the sense of volumetric scaling. This together with Figs. 1-2 implies that inhomogeneous pair production is a surface effect as it was already predicted by Heisenberg in 1934 when he calculated the fluctuation of a charge in QED and found that it was ∝ V 2 3 or the surface of the interaction volume [35]. This is in a disagreement to what one may conclude from the constant-homogeneous string models usually applied in heavy ion physics.

Transverse spectra
The p y momentum is fully conserved in this setup so it is just an additional mass term and therefore set to 0 in the following. We integrate the x coordinate and plot the p z − p x momentum spectra in Fig. 5. In the homogeneous case it would be a simple radially symmetric peak as larger  Fig. 2. Note that the right side of the distribution is wider than the left. The distribution is accelerated to the left. momenta are exponentially suppressed. However in this case particles are also created by the inhomogeneity and they depart from the main peak and this increases the production rate in the transverse direction. The earlier the particles produced the further they get from the center and with the additional effect of accelerating in the E field this gives a triangular shape to the distribution.
The angular distribution of pairs is an important observable. It can be computed by the momentum integral in polar coordinates n(θ) = ∞ 0 dpdx 2π pf (x, p x = p sin(θ), p y = 0, p z = p cos(θ)).
(24) The resulting distribution can be seen in Fig. 6. For the homogeneous case this would be a simple peak as illustrated by the case τ = λ c /c, but the excess particles provided by the inhomogeneity for τ = 2λ c /c and τ = 3λ c /c give rise to two side peaks while the central peak shrinks due to its increasing distance from the origin. The appearance of side peaks is remarkably similar to the predicted bifurcation for squeezed states in quantum optics, see Ref. [36].

Summary
We investigated the Abelian pair production in inhomogeneous external electric field with a shape motivated by the color flux tubes in heavy ion collisions. We found that the number of particles created can be significantly underestimated by the homogeneous models. Moreover we showed evidence that the results obtained with homogeneous string models can be conceptually misleading and a proper description including finite size effects may needed. Such a description may provide microscopical model for the empirical parameters. We also presented the widening of the transverse spectra as a potential discriminant of homogeneous and inhomogeneous processes.