The inverse problem for Schwinger pair production

The production of electron-positron pairs in time-dependent electric fields (Schwinger mechanism) depends non-linearly on the applied field profile. Accordingly, the resulting momentum spectrum is extremely sensitive to small variations of the field parameters. Owing to this non-linear dependence it is so far unpredictable how to choose a field configuration such that a predetermined momentum distribution is generated. We show that quantum kinetic theory along with optimal control theory can be used to approximately solve this inverse problem for Schwinger pair production. We exemplify this by studying the superposition of a small number of harmonic components resulting in predetermined signatures in the asymptotic momentum spectrum. In the long run, our results could facilitate the observation of this yet unobserved pair production mechanism in quantum electrodynamics by providing suggestions for tailored field configurations.


Introduction
The vacuum breakdown in external electric fields by the emission of the lightest charged particle-antiparticle pairs (Schwinger effect) is an unobserved prediction of quantum electrodynamics (QED) [1,2,3]. Until recently, this pair production mechanism seemed to be a mere academic problem owing to the required electric field strength of the order of E S ∼ 10 18 V/m. Recent efforts have raised the hopes that a direct observation might become feasible at ultra-high intensity laser facilities [4,5,6,7]. Alternatively, there have been suggestions to study the Schwinger effect in graphene [8,9,10] or to quantum simulate it in analogue cold atom systems [11,12].
The pair production problem in QED cannot be solved in full glory but requires an approximate treatment. Assuming ultra-high intensity lasers, it is well-justified to neglect quantum fluctuations of the electromagnetic field and to treat it as a classical, external one. In addition, since the spatio-temporal laser scales are orders of magnitude different from the QED scale, which is set by the electron mass m, the focal region is typically modeled by a spatially homogeneous, time-dependent electric field. In this case, quantum kinetic equations are especially well-suited to study the pair production problem numerically since they can be formulated as a coupled system of ordinary first-order differential equations [13,14], which can be solved very efficiently. For more complicated field configurations including spatial inhomogeneities or non-vanishing magnetic fields, more sophisticated approaches such as worldline methods [15,16,17,18], the Dirac-Heisenberg-Wigner formalism [19,20,21,22] or methods from real-time lattice gauge theory [23,24,25,26] have been employed.
It has been proposed to lower the threshold for pair production by superimposing electric fields with different frequencies Email address: hebenstreit@itp.unibe.ch (F. Hebenstreit) (dynamically-assisted Schwinger effect) [27] and the resulting momentum spectra have been investigated. Several studies revealed a very rich structure of the momentum distribution, owing to the fact that the pair-production problem depends nonlinearly on the applied field profile. In fact, certain characteristic features of the momentum distribution were interpreted in terms of symmetries of the electric field and resonance conditions [28,29,30,31,32,33,34,35,36]. However, it is almost impossible to predict the variations in the momentum spectrum upon changing the field parameters a priori, i. e. without actually solving the dynamic equations. A rare exception are solitonic fields, for which particle production can be excluded at one predetermined momentum mode [37,38].
From an experimental point of view, we have to face the challenge that the number of produced electron-positron pairs is exponentially small and therefore hard to discriminate from background noise. Moreover, spectrometers are only sensitive in a certain momentum window and characterized by a finite momentum resolution. Accordingly, it would be most beneficial to generate electric field configurations such that the momentum distribution is peaked in a characteristic momentum window in order to enhance the detection probability. This situation is reminiscent of the inverse scattering problem in quantum mechanics where asymptotic scattering data is used to reconstruct the scattering potential [39].
In this manuscript, we present a method to solve the inverse problem for Schwinger pair production, i. e. we determine the field parameters resulting in a predetermined momentum distribution, by combining quantum kinetic theory with optimal control theory. We employ optimization techniques which have been utilized previously in performing pulse-shape optimization for pair production [40,41], and which have also proven successful in the closely related field of atomic, molecular and optical physics [42,43,44]. Motivated by advances in high harmonic and attosecond pulse generation [45,46,47,48], we superimpose a small number of harmonic components and optimize their amplitudes such that a predetermined momentum distribution is well-approximated. In the long run, our results could be used for tailoring superpositions of high harmonics in order to maximize the detection probability in a given momentum window. On the other hand, a measured momentum distribution could also serve as a tomograph of the laser field, which is hard to control in an absolute manner at ultra-high intensities.

Inverse problem for Schwinger pair production
Vacuum electron-positron production in a unidirectional time-dependent electric field E(t) = −Ȧ(t) can be described in terms of a quasi-particle distribution function F(q, t), where q ≡ (q , q ⊥ ) is the momentum variable with components along and perpendicular to the field direction, and t denotes time.
can be regarded as the momentum distribution of asymptotic particles [49]. Defining auxiliary functions G(q, t) and H(q, t), its dynamics is governed by the coupled system of ordinary first-order differential equations [14] In the following, we take a superposition of n harmonic components on the time interval t ∈ [−T, T ] such that E(±T ) = 0 and A(−T ) = 0. We then mean to select the field amplitudes E ≡ (E 1 , . . . , E n ) ∈ Ê n such that a predetermined distribution function F 0 (q) is obtained at the final time T . While this search could be achieved by brute force for a very small number of harmonics, this approach becomes practically unfeasible for larger values of n. Hence we set up an optimization problem to solve this inverse problem for Schwinger pair production by defining the cost functional as the least square deviation of F(q, T ) from the objective distribution F 0 (q): The global minimum of the cost functional in the space Ê n ∋ E spanned by the field amplitudes is then defined as For reasons which will become clear later, the cost functional will also be denoted as quality factor in the following. We emphasize that the distribution function F(q, T ) implicitly depends on the field amplitudes E = (E 1 , . . . , E n ) via the dynamical system (1). The constrained optimization problem can be recast in an unconstrained optimization problem by introducing Lagrange multiplier functions λ F,G,H (q, t) [40], which need to obey the adjoint equationsλ The stationary condition 0 = ∇η n [E * ] ∈ Ê n is then a necessary condition for being a local extremum, with Here, ∇ ≡ (∂ E 1 , . . . , ∂ E n ) is the gradient with respect to the field amplitudes E j , and E * denotes a local minimizer configuration. Moreover, all occurring functions are supposed to be solutions of the dynamical equations (1) and (5), respectively, and the field gradients ∇E, ∇A ∈ Ê n are easily determined from (2).
To take full advantage of the gradient information (6), we employ a multi-start method along with a local optimization algorithm. To this end, we generate random trial configurations E 0 which are sampled from a probability distribution in field amplitude space. The corresponding local minimizer configurations E * = lim k→∞ E k are then found iteratively with k ∈ AE 0 . We calculate the search directions d k according to the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, and viable step sizes α k are found via an inexact line search fulfilling the strong Wolfe conditions. For further algorithmic details we refer to, e. g. reference [50]. As a word of caution we emphasize that it is only guaranteed for infinite sample sizes that the global minimumη n ≡ η n [Ẽ], whereẼ denotes the global minimizer configuration, is among the detected stationary points [51]. Moreover, it is also not guaranteed that an exact solution of the inverse problem, which is characterized by a vanishing cost functionalη n = 0, actually exists. In general, we will rather findF(q, T ) F 0 (q) for certain momenta, whereF(q, T ) denotes the momentum distribution corresponding to the global minimumη n . This is in fact a sign that the inverse problem for Schwinger pair production is in general ill-posed [52]. A configuration that approaches the exact solution to the inverse problem is characterized byη n ≪ 1, and we will study its dependence on the number of harmonic components n in the following.

Solving the inverse problem
For simplicity, we study the inverse pair production problem in a 1-dimensional system but emphasize that the general procedure is basically not altered in higher dimensions. To be specific, we predefine a Gaussian model objective distribution and we take as parameters A 0 = 10 −5 , q 0 = 5m and σ 0 = m. We could specify any functional form in principle, however, we will restrict ourselves to this simple and instructive example. Specifically, we consider the harmonic superposition (2) in the time interval 2T = 100/m. We first consider the superposition of n = 2 harmonics and calculate the landscape η 2 [E 1 , E 2 ] by brute force, cf. Fig. 1. We note that the distribution function F 2 (q, T ) is computed for an array of field amplitudes E i ∈ {−0.275E S , . . . , 0.275E S } with ∆E i = 0.005E S , resulting in O(10 4 ) simulated data points. We therefore conclude that the brute force approach becomes practically unfeasible for large values of n. 1 The landscape shows a diamond-shaped region with quality factors η 2 [E] = O(1) which is rather flat, whereas it increases exponentially outside. On the other hand, we also find a small number of distinct local 1 Given the chosen time extent T and requiring sufficient momentum resolution, it takes a few minutes on a standard desktop computer to calculate the momentum distribution F 2 (q, T ) according to (1). Consequently, the computation time on the chosen array of field amplitudes is of the order of T 2 ∼ O(10 3 ) CPU hours. Increasing the number of harmonics n > 2 and keeping the resolution ∆E i unaltered, the computation time grows exponentially T n ∼ 100 (n−2) T 2 , indicating that the brute force approach becomes practically unfeasible for large values of n. minima which are indicated in green in Fig. 1, amongst them the global minimumẼ = (Ẽ 1 ,Ẽ 2 ) = (0.1644E S , −0.1238E S ) withη 2 = 0.1592 > 0. We emphasize that there is only an approximate solution to the inverse problem as the quality factor is non-vanishing.
To illustrate the functionality of the employed optimization method, we first perform the numerical optimization for n = 2 and we display some of the calculated trajectories in Fig. 1. Depending on the randomly chosen initial configuration E 0 , the various trajectories converge to the diverse stationary points in the landscape η 2 [E * 1 , E * 2 ]. For n = 2, the average number of required iteration steps for sufficient convergence -based on N 0 = 200 randomly chosen initial configuration -isk 2 = 17.5.
In order to converge towards the objective momentum distribution F 0 (q), i. e. to further decrease the quality factorη n , we need to increase the number of harmonic components n. Most notably, we then observe quick and substantial improvement of η n . This assertion is quantified in Fig. 2, where we display the quality factorη n as a function of n. Moreover, we display the underlying global minimum distributionsF n (q, T ) for different values of n along with the corresponding electric field configurations E(t) in Fig. 3. Based on our result, we note that the average number of required iteration steps grows likek n ∼ n, indicating that the simulations become more expensive for increasing n. We emphasize, however, that this linear growth of the optimization approach outweighs by far the exponential growth of the brute force approach.
For n = 1 we find an optimal distribution functionF 1 (q, T ) which exhibits a broad central peak accompanied by two narrow side peaks. Notably, the central peak is located around q 1 = 7.4m and is thus substantially shifted from q 0 = 5m. Including one additional harmonic such that n = 2, the optimal distribution functionF 2 (q, T ) becomes peaked around q 2 = 5m = q 0 but is still somewhat broader and smaller than F 0 (q). Notably, we still observe two undesired peaks around q = 0 and q = 10.5m, respectively. Further increasing the number of harmonics n > 2, the optimal distributionF n (q, T ) approximates the objective distribution F 0 (q) better and better. Here, we emphasize that the global minimum configurationF n (q, T ) does not necessarily arise from the global minimum configurations F n−1 (q, T ), meaning that (Ẽ n−1 , 0) ∈ Ê n does not necessarily lie in the basin of attraction ofẼ n ∈ Ê n . For instance, given the n = 2 global minimum configurationẼ 2 = (0.1644E S , −0.1238E S ) withη 2 = η 2 [Ẽ 2 ] = 0.1592 (cf. Fig.1) and taking these field parameters as the initial point for the n = 3 optimization problem, we converge towards the stationary point  Fig. 3).
The presented results indicate that a rather small number of harmonic components n suffices to well-approximate specific predetermined features of the momentum distribution. For instance, the distribution around q 0 = 5m is already extremely well-described for n = 4 whereas there are still comparatively large deviations at other momentum values present. On the other hand, the number of harmonics n needs to be further increased in order to approach the exact solution of the inverse problem for Schwinger pair production. Quantitatively, we find that the quality factor drops fromη 1 = 0.7916 tõ η 10 = 0.9525 · 10 −6 ≪ 1. Based on our results, the functional dependence of the quality factor is reasonably well-described by an exponential functionη n ∼ 4.1(2) × exp (−1.64(4)n). One may therefore speculate that there is an exact solution of the inverse problem for the chosen objective distribution F 0 (q). It is, however, not clear whether this solution then exists for finite n or only in the limit n → ∞.
We already mentioned that the employed multi-start method is guaranteed to converge to the global minimum only for infinite sample sizes. Accordingly, Fig. 2 is to be considered as an upper limit whereas there are, in principle, even better configurations conceivable. In the current study, we have used N 0 = 200 random initial configurations E 0 for all considered values of n. Given this ensemble size, all the numerically determined global minimum distributionF n (q, T ) were encountered several times. We are therefore confident that the identified configurations actually correspond to the global minimumη n or lie, at least, in its close neighborhood.

Conclusions
We demonstrated that quantum kinetic theory along with optimal control theory can be utilized to solve the inverse problem for Schwinger pair production. Based on the instructive example of a Gaussian objective distribution F 0 (q), we found that it suffices to superimpose a comparatively small number of harmonics n = 4 in order to well-approximate specific predetermined features. These findings could substantially facilitate the observation of Schwinger pair production by providing suggestions for tailored field configurations. In the long run, the presented method could also serve as a tomograph for applied laser pulses by reconstructing the applied field after measuring the asymptotic momentum distribution. One has to be aware, however, that this procedure might be non-unique due to the ill-posedness of the inverse problem. On the other hand, it is not clear at this point whether an exact solution of the inverse problem withη n = 0 actually exists and, if so, a finite number of harmonic components n suffices. This issue is beyond the scope of the current investigation but should be further studied in the future.
The used optimization method is especially well-suited for parameter spaces which are not too high-dimensional. In fact, we found that the average number of required iterations grows likek n ∼ n, indicating that the optimization algorithm becomes less efficient for increasing n. Moreover, the ensemble size N 0 of random initial configuration E 0 needs to be enlarged upon further increasing the number of harmonics n, indicating that the multi-start approach becomes less efficient as well. Accordingly, we should also envisage to apply global optimization strategies based on metaheuristic algorithms in the future in order to study the inverse problem for more intricate objective distributions F 0 (q), to investigate the possibility of an exact solution of the inverse problem withη = 0, or to perform elaborate pulse shaping investigations.