In this article, we review recent theoretical works on the Schwinger mechanism of particle production in external electrical fields. Although the non-perturbative Schwinger mechanism is at the center of this discussion, many of the approaches that we discuss can cope with general time and space dependent fields, and therefore also capture the perturbative contributions to particle production.

Introduction

The Schwinger mechanism [START_REF] Schwinger | On gauge invariance and vacuum polarization[END_REF] is a non-perturbative phenomenon by which electron-positron pairs can be produced by a static electrical field [START_REF] Schwinger | On gauge invariance and vacuum polarization[END_REF][START_REF] Sauter | Uber das Verhalten eines Elektrons im homogenen elektrischen Feld nach der relativistischen Theorie Diracs[END_REF][START_REF] Heisenberg | Consequences of Dirac's theory of positrons[END_REF]. It is the fact that the field is time independent that makes this process genuinely non-perturbative. Indeed, such an electromagnetic field can be viewed as an ensemble of photons of zero frequency, which by energy conservation are unable to produce massive particles. The non-perturbative nature of this phenomenon is visible in the fact that the corresponding transition amplitudes are non analytic in the coupling: since they are proportional to exp(-πm 2 /eE), the coefficients of their Taylor expansion around e 2 = 0 are all zero, and only calculations to all orders in the coupling can make this effect manifest. The non-analytic behavior of these transition probabilities can be understood as a tunneling phenomenon (see for instance [4]) by which a particle from the Dirac sea is pulled into the positive energy states: as illustrated in the figure 1, the external gauge potential tilts the gap between the Dirac sea and the positive energy states, allowing a hole from the sea to tunnel through this gap and materialize as an on-shell positive energy particle on the other side. From this geometrical interpretation, it is clear that the length that needs to be crossed is inversely proportional to the electrical field, leading to a tunneling probability that is exponential in 1/E.

From this pocket formula, we also see that the external field must be very intense in order to lead to a significant probability of particle production: at small coupling, this probability is of order one only for fields that are comparable to a critical field E c inversely proportional to the coupling constant, E ∼ E c ≡ m 2 /e. Even for the lightest charged particle -the electron-, these are extremely intense fields that still surpass by several orders of magnitude the largest fields achievable experimentally. Note that when one studies the interaction of a charged particle with such an external electrical field, the corresponding expansion is in powers of eE. However, if eE is comparable to the squared mass of the particle, then there is no small parameter to control this expansion and one must treat the external field to all orders. This is a prerequisite of any formalism aiming at describing the Schwinger mechanism.

Although the original derivation of the Schwinger mechanism was made in the context of quantum electrodynamics, it can be generalized to any quantum field theory coupled to some external field. In the context of strong interactions and quantum chromo-dynamics, where the coupling constant g is numerically much larger than in QED, the production of pairs by the Schwinger mechanism may be achievable with more moderate chromo-electrical fields, and this phenomenon may play a role in the discussion of particle production in heavy ion collisions [START_REF] Biro | Color Rope Model for Extreme Relativistic Heavy Ion Collisions[END_REF][START_REF] Kajantie | Decay of Strong Color Electric Field and Thermalization in Ultrarelativistic Nucleus-Nucleus Collisions[END_REF][START_REF] Gatoff | The Flux Tube Model for Ultrarelativistic Heavy Ion Collisions: Electrohydrodynamics of a Quark Gluon Plasma[END_REF][START_REF] Kharzeev | Multi-particle production and thermalization in high-energy QCD[END_REF], or in the decay of "hadronic strings" [START_REF] Andersson | Parton Fragmentation and String Dynamics[END_REF] in the process of hadronization. However, the gauge fields that are generated in heavy ion collisions have two important features: (i) since gauge fields have a direct coupling to gluons, the leading contributions for gluon production in heavy ion collisions are tree level contributions, that supersede the 1-loop contributions encountered in the Schwinger mechanism; (ii) these fields are in general space and time dependent [START_REF] Mclerran | Computing quark and gluon distribution functions for very large nuclei[END_REF][START_REF] Mclerran | Gluon distribution functions for very large nuclei at small transverse momentum[END_REF], and one must therefore use a formalism that can cope with the most general type of external field, without assuming any symmetry in its spacetime dependence. Naturally, such fields generally imply that it is no longer possible to obtain closed analytic expressions (such analytical results have been obtained only for very special time dependences, and for spatial dependences that possess a high degree of symmetry). As a consequence, it is usually necessary to resort to numerical studies, and a recurrent concern in this review will be the practicality of various approaches for these numerical simulations. The Schwinger mechanism can also be discussed in the context of a much simpler scalar field theory, coupled either to a scalar or vector external field (as in scalar QED). Despite its lack of connection to possible experimental realizations, this is the simplest example one may think of and it offers a very useful playground for testing new theoretical developments. In this review, we will often use such toy theories to illustrate various novel approaches, because of their didactic or phenomenological interest.

The outline of this review is as follows. In the section 2, we consider a scalar theory coupled to an external source and discuss particle production at tree level and 1-loop. The goal of this section is to present general ideas about these processes, whose range of validity is much broader than the simple toy model used to introduce them. The section 3 is devoted to a discussion of the multi-(anti)particle correlations that exist in the Schwinger mechanism. We first derive their general structure and then work it out in the special case of spatially homogeneous fields (possibly time dependent). The general discussion of the section 2 leads to a formulation of the Schwinger mechanism in terms of a complete basis of mode functions. In the section 4, we relate this representation to other approaches: the method of Bogoliubov transformations, the quantum kinetic equations, and the Wigner formalism, and we discuss lattice numerical implementations of this approach in the section 5. The section 6 is devoted to the worldline formalism, a radically different (but equivalent) formulation of the Schwinger formalism based on Schwinger's proper-time representation of propagators in an external field. Besides new methods for calculating particle production in an external field, this approach provides a great deal of intuition on the spatio-temporal development of the production process. In the section 7, we discuss the idea of dynamically assisted Schwinger mechanism, where one superimposes two fields that have vastly different timescales and magnitudes in order to reach particle yields that are much larger than what would have been achieved with each field separately.

Quantum fields coupled to external sources

In this section, we discuss general aspects of quantum field theories coupled to an external source [START_REF] Gelis | Particle production in field theories coupled to strong external sources[END_REF][START_REF] Gelis | Particle production in field theories coupled to strong external sources. II. Generating functions[END_REF]. Our main goal is to present the general aspects of such theories, focusing on their main differences with field theories where the sole interaction are the self-interaction of the fields. For the sake of simplicity, we consider a scalar field theory with a quartic self-interaction, whose Lagrangian is given by

1 L ≡ 1 2 (∂ µ φ)(∂ µ φ) - 1 2 m 2 φ 2 - g 2 4! φ 4 + jφ , (1) 
where j(x 0 , x) is an unspecified function of space-time. Note that this source is a commuting number-valued object, rather than an operator. In the same way that one usually assumes that self-interactions are adiabatically turned on and off when x 0 → ∓∞, we assume that the external source decreases fast enough when time goes to ±∞. Moreover, in order to preserve the unitarity of the theory, the external source must be real valued. Since this section is devoted to a general discussion, we allow the particles to have self-interactions. This is important for applications to gluon production in heavy ion collisions, and it leads to some complications. In contrast, the production of electron-positron pairs in QED is simpler to study because the electrons can only interact directly with the photon field. Interactions between electrons and positrons can happen indirectly, with the mediation of a photon, but this effect is an extremely small correction in practice, that would arise beyond the order considered here.

Vacuum diagrams

The main feature of such a theory is that it describes an open system: even if the system is initialized in the vacuum state that contains no particles, the external source can -and in general will-produce particles. This is in sharp contrast with the same theory in the vacuum, where the in-vacuum state and the out-vacuum state are related by a unitary transformation Ω(+∞, -∞) :

0 out = Ω(+∞, -∞) 0 in (when j ≡ 0) . (2) 
This property ensures that when j ≡ 0, the vacuum state evolves with probability one into the vacuum state, i.e. that no particle is created. An equivalent statement is that the vacuum to vacuum transition amplitude is a pure phase:

0 out 0 in = e i V with V ∈ R . (3) 
This has some important practical consequences. The perturbative expansion for transition amplitudes generates diagrams that contain disconnected vacuum sub-diagrams, i.e. diagrams that have no external legs, but the above property tells us that the sum of all the vacuum diagrams is a pure phase that, although it appears as a prefactor in every transition amplitude, does not play any role after squaring the amplitudes in order to obtain transition probabilities. Therefore, it is legitimate to ignore from the start all the graphs that contain vacuum sub-diagrams in the j ≡ 0 case. Let us return to the case j = 0 and denote P (α) the transition probability from the initial vacuum state 0 in to an arbitrary final state α out . The unitarity of the theory implies that the sum of these probabilities over all the possible final states is equal to unity, α P (α) = P (0) +

α =0 P (α) = 1 , (4) 
where in the first equality we have separated the vacuum and the populated states. If particles are produced during the time evolution of the system (this is possible only if j = 0), then there is at least one non-empty state α for which P (α) > 0. Since all these probabilities are numbers in the range [0, 1], the previous identity implies that P (0) < 1, i.e. that the vacuum does not evolve into the vacuum state with probability one anymore. A trivial consequence is that the vacuum to vacuum transition amplitude is no longer a pure phase, since its square is strictly smaller than unity. Because of this, it is no longer possible to disregard the vacuum graphs, as illustrated in the figure 2. In general, this complicates Figure 2: Example of diagram contribution to the probability of producing 11 particles (here illustrated in the case of a scalar theory with a cubic coupling to simplify the diagrams).

The black dots denote the insertions of the external source j(x). The vertical line indicates the final state (left of this line: amplitude, right: complex conjugated amplitude).

considerably the diagrammatic expansion when j = 0, but we will see later that inclusive quantities have a diagrammatic expansion made only of connected diagrams. Moreover, in the case of a strong external source (of order 1/g), it becomes hopelessly complicated to calculate exclusive quantities (such as the probability P (α) for an individual final state) while the inclusive ones are much easier to access.

Power counting

In order to discuss ways of organizing the calculation of observables in field theories coupled to an external source, the first step is to assess the order of magnitude of a graph in terms of its topology. Obviously, for a graph with multiple disconnected components such as the one in the figure 2, the order is obtained as the product of the orders of each of its connected subdiagrams. Therefore, it is sufficient to consider only connected graphs in this discussion. A connected diagram is fully characterized by the number of sources n J , the number of propagators n P , the number of vertices n V , the number of loops n L and the number of external particles n E (i.e. propagators endpoints that are not connected to a vertex or to a source). These quantities are not all independent. A first constraint comes from the fact that each propagator has two endpoints and each vertex receives four lines (for the scalar theory described by the Lagrangian of eq. ( 1)):

2n P = 4n V + n J + n E . (5) 
A second identity expresses the number of independent loops n L in terms of the other characteristics of the graph2 :

n L = n P -n E -n J -n V + 1 . (6) 
Thanks to these two formulas, the order of a connected graph can be written as

ω(G) ≡ g 2n V j n J = g -2+n E +2n L (gj) n J . (7) 
In this expression, we have combined one power of the coupling g with each power of the external source j, because the combination gj disappears from this power counting formula in the strong source regime where j ∼ g -1 (i.e. when particle production by the Schwinger mechanism becomes likely). Eq. ( 7) displays a standard dependence on the number of loops and external legs. But it also indicates that in the strong source regime the dependence of the order on n J disappears: therefore, infinitely many graphs contribute at each order. In conjunction with our previous observation that vacuum graphs cannot be disregarded, this poses a serious bookkeeping challenge for the calculation of quantities such as the transition probability shown on the figure 2.

Exclusive and inclusive quantities

Earlier, we have alluded to important differences between inclusive and exclusive observables when it comes to calculating them. Before going into these technical differences, let us first define them more precisely. Loosely speaking, exclusive observables are related to a full measurement of the final state, while inclusive observables involve dropping a lot of information about the final state.

If we assume the initial state to be empty, all the information about the final state is encoded in the transition amplitudes [START_REF] Kharzeev | Multi-particle production and thermalization in high-energy QCD[END_REF] in which the final state is fully specified by giving the list of the momenta of all the final particles (in theories with more structure than the scalar theory under consideration in this section, one would also need to specify other quantum numbers of the final particles). Any observable related to this evolution can be expressed in terms of these amplitudes. For instance, the differential probability for producing n particles can be expressed as

p 1 • • • p nout 0 in
dP n d 3 p 1 • • • d 3 p n = 1 n! 1 (2π) 3 2E p1 • • • 1 (2π) 3 2E pn p 1 • • • p nout 0 in 2 , (9) 
where E p ≡ p + m 2 is the on-shell energy of a particle of momentum p. This quantity is called exclusive because only a single final state can contribute to it, at the exclusion of all others.

The archetype of inclusive observables are the particle spectra, obtained from the above probability distributions by integrating out the phase-space of all particles but a few. The simplest of them is the single particle spectrum, defined as

dN 1 d 3 p ≡ ∞ n=0 (n + 1) d 3 p 1 • • • d 3 p n dP n+1 d 3 pd 3 p 1 • • • d 3 p n . ( 10 
)
This quantity gives the number of particles (thanks to the factor n + 1 included under the sum) produced in a given momentum range, but the information about the distribution of individual final states has been lost. As we shall see, this quantity is much easier to calculate than the exclusive observables. Experimentally, it is also much more accessible since it only requires to make an histogram of the momenta of the final state particles. The single particle spectrum has obvious generalizations: the 2-particle spectrum, the 3-particle spectrum, etc... that provide information about the correlations between the final state particles.

The "vacuum survival probability" often plays a special role in discussions of the Schwinger mechanism. This quantity is nothing but

P 0 ≡ 0 out 0 in 2 . (11) 
A standard result in field theory is that the vacuum transition amplitude 0 out 0 in can be written as the exponential of the sum of the connected vacuum diagrams,

0 out 0 in = e i V , i V ≡ + + + + . . . 1 6 1 8 1 8 (12) 
(Only a few tree level contributions to V have been shown in the illustration.) After squaring this amplitude, we get P 0 = exp(-2 Im V). In the presence of an external source that remains constant over a long interval of time and is sufficiently spatially homogeneous, the imaginary part 2Im V can usually be written as an integral over space-time,

2Im V = d 4 x (• • • ) , (13) 
and the integrand in this formula is often interpreted as the "particle production rate". We will return on this interpretation later, as it not entirely accurate in general and can be a source of confusion.

Cutting rules and Schwinger-Keldysh formalism

Earlier in this section, we have seen that the imaginary part of the connected vacuum graphs, Im V, plays an important role in the expression of the transition probabilities. This imaginary part can be obtained from an extension of the standard Feynman rules known as Cutkosky's cutting rules [START_REF] Cutkosky | Singularities and discontinuities of Feynman amplitudes[END_REF][START_REF] Hooft | Diagrammar[END_REF]. Here, we do not re-derive these rules but simply state them as a recipe, that must be applied to every graph:

i. Divide the nodes (vertices and sources) contained in the graph into a set of + nodes and a set ofnodes, in all the possible ways. It is customary to materialize diagrammatically these two sets of nodes by drawing a line (the "cut") that divides the graph in two subgraphs. Even if the original graph is connected, the sub-graphs on each side of this cut do not have to be connected.

ii. The + vertices give a factor -ig 2 and thevertices give a factor +ig 2 . The + sources give a factor +ij(x) and thesources give a factor -ij(x).

iii. Two nodes of types = ± and = ± are connected by a bare propagator G 0 . In momentum space, these four propagators read

G 0 ++ (p) = i p 2 -m 2 + i , G 0 --(p) = -i p 2 -m 2 -i G 0 +-(p) = 2πθ(-p 0 )δ(p 2 -m 2 ) , G 0 -+ (p) = 2πθ(+p 0 )δ(p 2 -m 2 ) . ( 14 
)
iv. Multiply the outcome of these rules by 1/2.

The subgraph that involves only + labels is obtained by the usual Feynman rules used to calculate transition amplitudes, while thesubgraph is given by the complex conjugation of these rules. The interpretation of these cutting rules is straightforward: the + sector corresponds to an amplitude and thesector to a complex conjugated amplitude, while the +-and -+ propagators provide the phase-space integration for the (on-shell, hence the delta functions) final state particles. In fact, these rules are the perturbative realization3 of the optical theorem, that stems from unitarity. If we write the S matrix as S ≡ 1 + iT , then S † S = 1 can be rewritten as

i(T † -T ) = T † T . ( 15 
)
By taking the expectation value of this identity in the vacuum state and inserting a complete set of states in the right hand side, it leads to

Im 0 in T 0 in V = 1 2 α 0 in T † α in α in T 0 in = 1 2 α +- 0 in α out -sector α out 0 in + sector . (16) 
The factor 1/2 in the right hand side of this identity is the origin of the rule iv above. The Schwinger-Keldysh [START_REF] Schwinger | Brownian motion of a quantum oscillator[END_REF][START_REF] Keldysh | Diagram technique for nonequilibrium processes[END_REF] formalism is essentially equivalent to Cutkosky's cutting rules, but it is usually introduced as a tool for the perturbative calculation of squared matrix elements rather than as a technique for calculating the imaginary part of a scattering amplitude (however, the optical theorem states that the two are closely related).

General remarks on the distribution of produced particles

All the transition amplitudes such as [START_REF] Kharzeev | Multi-particle production and thermalization in high-energy QCD[END_REF] contain in their diagrammatic expansion the disconnected vacuum graphs, as shown in the example of the figure 2. In other words, one can pull out a factor exp(iV) in all these amplitudes, or a factor exp(-2 Im V) in squared amplitudes. From our previous power counting formula [START_REF] Gatoff | The Flux Tube Model for Ultrarelativistic Heavy Ion Collisions: Electrohydrodynamics of a Quark Gluon Plasma[END_REF], we expect Im V to start at the order 1/g 2 in the strong source regime,

2Im V = a g 2 , P 0 = e -a/g 2 , ( 17 
)
where a is an infinite series in powers of g 2n L (gj) n J . At tree level (n L = 0) and in the strong source regime, a is therefore of order unity (but depends non-perturbatively on the source).

Let us now turn to the probabilities of producing particles. Besides the factor exp(iV), the transition amplitude from the vacuum to a state containing one particle is made of graphs that connect sources to a single final particle. From eq. ( 7), these graphs start at the order g -1 , and the leading contribution to P 1 is of the form

P 1 = e -a/g 2 b 1 g 2 , ( 18 
)
where b 1 /g 2 is the sum of the 1-particle cuts (i.e. cuts that cut exactly one propagator) through connected vacuum diagrams. Its diagrammatic expansion starts with the following terms b 1

g 2 = -+ 1 2 + - 1 2 + -+ 1 6 + - + 1 6 + - + 1 6 + + - 1 6 + + - 1 6 + + - 1 6 + • • • . ( 19 
)
The probability P 2 of producing two particles is a bit more complicated. It contains a term (b 1 /g 2 ) 2 /2! corresponding to the independent emission of two particles (the factor 1/2! is a symmetry factor due to the fact that the two particles are indistinguishable) and an additional term b 2 /g 2 in which the two particles are correlated:

P 2 = e -a/g 2 1 2! b 2 1 g 4 + b 2 g 2 . ( 20 
)
The quantity b 2 /g 2 can be obtained as the sum of the 2-particle cuts of connected vacuum diagrams:

b 2 g 2 = - + + - + 1 6 + - + 1 6 + + - 1 6 + + - 1 6 + + - 1 6 + • • • (21) 
Likewise, the probability P 3 of producing 3 particles reads

P 3 = e -a/g 2 1 3! b 3 1 g 6 + b 1 b 2 g 4 + b 3 g 2 , ( 22 
)
where the last term is the sum of the 3-particle cuts of connected vacuum diagrams

b 3 g 2 = 1 8 -+ + 1 8 - + + 1 8 + - -+ + • • • (23) 
The previous examples can be generalized easily into a formula for the probability of producing n particles,

P n = e -a/g 2 n p=1 1 p! r1+•••+rp=n b r1 • • • b rp g 2p . ( 24 
)
In this formula, p is the number of clusters (sets of correlated particles) into which the n particles can be divided. Note that, except for the powers of g 2 that indicate the order of each term in the strong source regime, this formula is completely generic and does not depend on the details of the field theory under consideration: it just expresses the combinatorics of grouping n objects into p clusters. Let us end this subsection by noting that unitarity requires that a = b 1 + b 2 + b 3 + • • • (this identity can be viewed as a consequence of the optical theorem, since a/g 2 is 2 Im (V) and (b 1 + b 2 + b 3 + ...)/g 2 is the sum of all the cuts through V).

Generating functional

So far, our discussion has been at a rather qualitative level. We now turn to a more detailed discussion of what quantities can be calculated and of what tools can do it. For bookkeeping purposes, it is useful to introduce the following generating functional,

F[z(p)] ≡ ∞ n=0 1 n! d 3 p 1 (2π) 3 2E p1 • • • d 3 p n (2π) 3 2E pn z(p 1 ) • • • z(p n ) p 1 • • • p nout 0 in 2 , (25) 
where z(p) is a test function over the 1-particle phase-space. Any observable which is expressible in terms of the transition amplitudes (8) can be obtained from derivatives of F[z(p)]. In particular, the differential probabilities and the single particle spectrum are obtained as:

dP n d 3 p 1 • • • d 3 p n = 1 n! δ n F[z(p)] δz(p 1 ) • • • δz(p n ) z(p)=0 , dN 1 d 3 p = δF[z(p)] δz(p) z(p)=1 . (26) 
As one can see on these two examples, observables in which the final state is fully specified correspond to derivatives evaluated at z(p) = 0 (which eliminates most of the final states from the sum in eq. ( 25)), while inclusive observables -in which all the final states are kept and most of their particles are integrated out-correspond to derivatives evaluated at z(p) = 1. This is in fact a general property. For instance, the 2-particle spectrum is

dN 2 d 3 pd 3 q = δ 2 F[z(p)] δz(p)δz(q) z(p)=1 , (27) 
and this formula has an obvious generalization to the case of the inclusive n-particle spectrum.

As we shall see, inclusive observables are much simpler to calculate than the exclusive ones. To a large extent, this simplification is due to unitarity. The simplest consequence of unitarity is

F[z(p) = 1] = 1 . (28) 
At z(p) = 0, one would have instead obtained F[z(p) = 0] = P 0 , which is a very complicated object. These considerations show that it is much simpler to study the generating functional near the point z(p) = 1 than near the point z(p) = 0. In terms of a diagrams, the identity (28) corresponds to an exact cancellation among an infinite set of diagrams when one evaluates F[z(p)] at z(p) = 1.

The reason why we discussed at length vacuum-vacuum diagrams in the previous subsection is that they play an important role in organizing the calculation of other quantities.

The key observation here is that the sum of the vacuum-vacuum diagrams is nothing but the generating functional for time-ordered Green's functions. More precisely, one has

0 out T φ(x 1 ) • • • φ(x n ) 0 in = δ iδη(x 1 ) • • • δ iδη(x n ) e i V[j+η] η=0 , (29) 
where the notation V[j + η] indicates that one should evaluate the connected vacuum diagrams with a fictitious source η added to the physical source j. The fictitious source is set to zero after having performed the functional differentiations. Then, it is easy to obtain a formal but useful formula for F[z(p)]. Start from the Lehmann-Symanzik-Zimmermann [4] reduction formula for the transition amplitude to a final state with n particles, 

p 1 • • • p nout 0 in = i n d 4 x 1 • • • d 4 x n e i(
×( x1 + m 2 ) • • • ( xn + m 2 ) 0 out T φ(x 1 ) • • • φ(x n ) 0 in .( 30 
)
By plugging eq. ( 29) in this reduction formula and squaring the result, we can write the squared amplitude as

p 1 • • • p nout 0 in 2 = C p1 • • • C pn e iV[j+η+] e -iV * [j+η-] η±=0 , (31) 
where the operator C p is defined by

C p ≡ d 4 xd 4 y e ip•(x-y) ( x + m 2 )( y + m 2 ) δ 2 δη + (x)δη -(y) . ( 32 
)
In the right hand side of eq. ( 31), the factors exp(iV) and exp(-iV * ) come respectively from the amplitude and its complex conjugate. Note that it is essential to keep their arguments distinct-hence the separate η + and η --so that the two derivatives in the operators C p act on different factors. One should set η ± to zero only after all the derivatives have been evaluated. The final step is to substitute eq. ( 31) into the definition [START_REF] Gelis | Color Glass Condensate and Glasma[END_REF] of F[z(p)]. One obtains immediately

F[z(p)] = exp d 3 p (2π) 3 2E p z(p) C p e iV[j+η+] e -iV * [j+η-] η±=0 . ( 33 
)
From this formula, one can show that the generating functional F[z(p)] is the sum of all the vacuum diagrams in a modified Schwinger-Keldysh formalism in which all the off-diagonal (G 0 +-(p) and G 0 -+ (p)) propagators are multiplied by z(p). There is no simple expression for the generating functional F[z(p)] itself, but it turns out that it is much easier to obtain a formula for its first derivative δF[z(p)]/δz(p). Using eq. ( 33) and the explicit form of the operator C p , we can write this derivative as

δ ln F[z(p)] δz(p) = 1 (2π) 3 2E p d 4 xd 4 y e ip•(x-y) ( x + m 2 )( y + m 2 ) × ϕ + (x)ϕ -(y) + G +-(x, y) . ( 34 
)
where ϕ ± (x) and G +-(x, y) are the connected 1-point and 2-point Green's functions in this modified Schwinger-Keldysh formalism, respectively (they implicitly depend on the external source j and on the test function z(p)), as illustrated in the figure 3. The order of magnitude of these objects is easily obtained from our general results for the power counting of connected graphs :

ϕ ± ∼ O(g -1 ) , G +-∼ O(1) . (35) 
Therefore, the first derivative of ln F[z(p)] starts at the order g -2 . Moreover, at this order, only the first term in ϕ + (x)ϕ -(y) contributes. The term in G +-starts contributing only at the next-to-leading order. A further simplification at leading order is that it is sufficient to keep tree level contributions to the 1-point functions ϕ ± . Thanks to this tree structure, the functions ϕ ± at Leading Order are solutions of the classical equation of motion,

( x + m 2 )ϕ ± (x) + U (ϕ ± (x)) = j(x) , (36) 
where U is the first derivative of the interaction potential (e.g. U (ϕ) = g 2 ϕ 3 /6 in the scalar model considered in this section). This equation depends on the external source j(x), but not on the test function z(p). The latter comes only via the boundary conditions

f (+) + (-∞, p) = f (-) -(-∞, p) = 0 , f (-) + (+∞, p) = z(p) f (-) -(+∞, p) , f (+) -(+∞, p) = z(p) f (+) + (+∞, p) . (37) 
In these equations, the boundary conditions have been written in terms of the coefficients of the Fourier decomposition of the fields ϕ ± ,

ϕ (y) ≡ d 3 p (2π) 3 2E p f (+) (y 0 , p) e -ip•y + f (-) (y 0 , p) e +ip•y . (38) 
(The Fourier coefficients are time dependent because ϕ ± are not free fields.) By plugging the Fourier representation of ϕ ± in the general formula [START_REF] Schenke | Fluctuating Glasma initial conditions and flow in heavy ion collisions[END_REF] (and setting G +-= 0 at this order), we get a very simple formula for the first derivative of ln F[z(p)] at leading order:

δ ln F[z(p)] δz(p) LO = 1 (2π) 3 2E p f (+) + (+∞, p) f (-) -(+∞, p) . ( 39 
)
Note that it is in general extremely difficult to solve the non-linear partial differential equation [START_REF] Gelis | High energy factorization in nucleus-nucleus collisions. 3. Long range rapidity correlations[END_REF] with boundary conditions imposed both at x 0 = -∞ and at x 0 = +∞, as in eqs. [START_REF] Voloshin | Transverse radial expansion in nuclear collisions and two particle correlations[END_REF]. Therefore, one should not hope to be able to find solutions of this problem (either analytically or numerically). As we shall see later, the only exception is when the fields under consideration do not have self-interactions but are solely driven by the external source. This is for instance the case for fermion production under the influence of an external electromagnetic field. Nevertheless, this result for the first derivative of the generating functional F[z(p)] is very useful as an intermediate tool for deriving other results, as will be shown in the rest of this section.

2.7

Inclusive quantities at leading order (tree level)

Single particle spectrum

Let us now show how to obtain inclusive moments (for now, at leading order) from eq. ( 39). The simplest one is the single inclusive spectrum. At leading order, it is simply obtained by evaluating eq. ( 39) at the special point z(p) = 1, since F[z(p) = 1] = 1. This means that one must solve the classical equation of motion with boundary conditions [START_REF] Voloshin | Transverse radial expansion in nuclear collisions and two particle correlations[END_REF] in which one sets z(p) = 1. Setting z(p) = 1 in these boundary conditions simplifies them considerably: the two fields ϕ + and ϕ -are identical,

if z(p) ≡ 1 , ϕ + (x) = ϕ -(x) ≡ ϕ(x) , (40) 
and obey the simple retarded boundary condition lim

x 0 →-∞ ϕ(x 0 , x) = 0 , lim x 0 →-∞ ∂ 0 ϕ(x 0 , x) = 0 . (41) 
Thus, the prescription for computing the single inclusive spectrum at leading order is the following:

i. Solve the classical field equation of motion with a null initial condition in the remote past,

ii. At x 0 → +∞, compute the coefficients4 f (±) (+∞, p) of the Fourier decomposition of this classical field,

iii. The single inclusive spectrum is then obtained as:

dN 1 d 3 p LO = 1 (2π) 3 2E p f (+) (+∞, p) 2 . (42) 
In eq. ( 42), we have used the fact that the retarded classical field ϕ is purely real5 . Since in the step i the boundary conditions are retarded, this problem is straightforward to solve, at least numerically. A few important comments are in order here:

• If the fields are not self-interacting, then the single particle spectrum is simply the square of the Fourier coefficients of the source itself.

• If the source has only space-like Fourier modes (this happens if there is frame in which it is time independent), then this is also the case for the solution of the classical equation of motion [START_REF] Gelis | High energy factorization in nucleus-nucleus collisions. 3. Long range rapidity correlations[END_REF], and the single particle spectrum is zero at leading order.

• This leading order contribution can only exist if the field under consideration is directly coupled to the external source. For instance, the direct production of electrons from an electromagnetic current is impossible at this order, but can happen at next-toleading order (the current couples to an electron-positron pair, but this requires an extra coupling constant). This is the form that gluon production takes in the Color Glass Condensate framework [START_REF] Iancu | Nonlinear gluon evolution in the color glass condensate. 1[END_REF][START_REF] Ferreiro | Nonlinear gluon evolution in the color glass condensate. 2[END_REF][START_REF] Iancu | The Color glass condensate: An Introduction[END_REF][START_REF] Iancu | The Color glass condensate and high-energy scattering in QCD[END_REF][START_REF] Weigert | Evolution at small x(bj): The Color glass condensate[END_REF][START_REF] Gelis | High energy scattering in Quantum Chromodynamics[END_REF][START_REF] Gelis | The Color Glass Condensate[END_REF][START_REF] Gelis | Color Glass Condensate and Glasma[END_REF] when applied to heavy ion collisions. In this effective description of high energy nucleus-nucleus collisions, the color gauge fields are coupled to two external currents that represent the color charges carried by the fast partons of the incoming nuclei. At high energy, it is expected that the gluon occupation number in nuclei may reach non-perturbative values of order 1/g 2 , which would correspond to strong sources in the sense used in this section. The spectrum of produced gluons at leading order in this description is obtained from the retarded classical solutions of the Yang-Mills equations,

D µ , F µν = J ν 1 + J ν 2 , (43) 
where J ν 1 and J ν 2 are the color currents of the two projectiles. This approach has been implemented in a number of works [START_REF] Kovner | Gluon production from nonAbelian Weizsacker-Williams fields in nucleus-nucleus collisions[END_REF][START_REF] Kovchegov | Classical gluon radiation in ultrarelativistic nucleus-nucleus collisions[END_REF][START_REF] Krasnitz | Nonperturbative computation of gluon minijet production in nuclear collisions at very high-energies[END_REF][START_REF] Krasnitz | The Initial energy density of gluons produced in very high-energy nuclear collisions[END_REF][START_REF] Krasnitz | Coherent gluon production in very high-energy heavy ion collisions[END_REF][START_REF] Lappi | Production of gluons in the classical field model for heavy ion collisions[END_REF][START_REF] Lappi | Some features of the glasma[END_REF][START_REF] Gelis | Glittering Glasmas[END_REF], and is now included in the IP-glasma model for the matter produced immediately after such a collision [START_REF] Schenke | Fluctuating Glasma initial conditions and flow in heavy ion collisions[END_REF].

Multi-particle spectra

The n-particle inclusive spectrum6 is also obtained from derivatives of the generating functional F[z(p)] evaluated at z(p) = 1,

dN n d 3 p 1 • • • d 3 p n = δ n F[z(p)] δz(p 1 ) • • • δz(p n ) z(p)=1 , (44) 
that can equivalently be written as :

dN n d 3 p 1 • • • d 3 p n = n i=1 δ ln F δz(p i ) z(p)=1 + i<j δ 2 ln F δz(p i )δz(p j ) z(p)=1 k =i,j δ ln F δz(p k ) z(p)=1 + • • • (45) 
The terms we have not written explicitly contain increasingly high order derivatives (but less and less factors), up to a single factor with an n-th derivative. However, these terms are not needed. Indeed, we already know that at leading order ln F is of order g -2 since it is a sum of connected vacuum-vacuum diagrams. Therefore, in the right hand side of this equation, the first term is of order g -2n , the second term is of order g -2(n-1) , etc... The leading contribution is thus the first term, and all the subsequent terms are subleading 7 . We see that, at leading order, the n-particle inclusive spectrum is simply the product of n single particle spectra:

dN n d 3 p 1 • • • d 3 p n LO = n i=1 dN 1 d 3 p i LO . ( 46 
)
d 3 p 1 • • • d 3 pn dNn d 3 p 1 • • • d 3 pn = ∞ N =n N (N -1) • • • (N -n + 1) P N ,
where P N is the total probability of producing exactly N particles. 7 The second term will play a role in the next-to-leading order corrections.

Any deviation from this factorized result is a subleading effect. Note also that at leading order, there is no difference between the factorial moments N (N -1) • • • (Nn + 1) and the ordinary moments N n . Moreover, at this order, the multiplicity distribution cannot be distinguished from a Poisson distribution.

In the Color Glass Condensate framework at leading order, the correlations among the produced particles can either originate from correlations that pre-exist in the distribution of the color sources that produce the gauge field [START_REF] Gelis | High energy factorization in nucleus-nucleus collisions. II. Multigluon correlations[END_REF][START_REF] Gelis | High energy factorization in nucleus-nucleus collisions. 3. Long range rapidity correlations[END_REF], or be built up at a later stage of the evolution of the system through collective motion of the produced particles (e.g. radial hydrodynamical flow [START_REF] Voloshin | Transverse radial expansion in nuclear collisions and two particle correlations[END_REF][START_REF] Shuryak | On the origin of the 'Ridge' phenomenon induced by jets in heavy ion collisions[END_REF]). In heavy ion collisions, strong correlations have been observed between pairs of hadrons [START_REF] Adare | Dihadron azimuthal correlations in Au+Au collisions at s(NN)**(1/2) = 200-GeV[END_REF][START_REF] Alver | System size dependence of cluster properties from two-particle angular correlations in Cu+Cu and Au+Au collisions at s(NN)**(1/2) = 200-GeV[END_REF][START_REF] Abelev | Long range rapidity correlations and jet production in high energy nuclear collisions[END_REF][START_REF] Alver | High transverse momentum triggered correlations over a large pseudorapidity acceptance in Au+Au collisions at s(NN)**1/2 = 200 GeV[END_REF], characterized by a ridge shape, very elongated in the relative rapidity of the two particles, and peaked in their relative azimuthal angle. By causality, the correlations in rapidity have to be created in the very early stages of the collisions [START_REF] Dumitru | Glasma flux tubes and the near side ridge phenomenon at RHIC[END_REF], and they can be simply understood as a consequence of the near boost invariance of the sources of the incoming nuclei. In contrast, the azimuthal correlations can be produced at any time, and are easily explainable by the hydrodynamical flow that develops in the later stages of the collision process [START_REF] Alver | Collision geometry fluctuations and triangular flow in heavy-ion collisions[END_REF][START_REF] Alver | Triangular flow in hydrodynamics and transport theory[END_REF].

Exclusive quantities at leading order

Let us now consider exclusive quantities. This discussion will be very short, and its purpose is only to illustrate the fact that the calculation of exclusive quantities is considerably more difficult than that of inclusive quantities. Let us consider as an example the calculation of the differential probability for producing exactly one particle. It may be obtained from F[z(p)] by the formula

dP 1 d 3 p = δF[z(p)] δz(p) z(p)=0 = F[z(p) = 0] P 0 δ ln F[z(p)] δz(p) z(p)=0 . ( 47 
)
There are two major differences compared to the inclusive spectra studied in the previous section :

i. The derivative of ln F[z(p)] must be evaluated at the point z(p) = 0. At leading order, it can still be expressed in terms of the Fourier coefficients of a pair of solutions of the classical equation of motion, via eq. [START_REF] Adare | Dihadron azimuthal correlations in Au+Au collisions at s(NN)**(1/2) = 200-GeV[END_REF]. However, because we must now set z(p) = 0 in the boundary conditions [START_REF] Voloshin | Transverse radial expansion in nuclear collisions and two particle correlations[END_REF] for these classical fields, they are not retarded fields anymore 8 , and there is no practical way to calculate them.

ii. The quantity F[z(p) = 0] appears as a prefactor in front of all the exclusive quantities. This prefactor is nothing but the probability P 0 for not producing anything, i.e. the vacuum survival probability. Calculating P 0 directly is a very difficult task. However, if one were able to calculate the second factor for all the probabilities P 1 , P 2 , • • • , one could then obtain P 0 from the unitarity condition ∞ n=0 P n = 1. These difficulties, observed here on the example of dP 1 /d 3 p, are in fact generic for all exclusive quantities. Note however that this is to a large extent an academic problem, since exclusive quantities -where one specifies in minute detail the final state-are not very interesting for the phenomenology of processes in which the final state has typically a very large number of particles, parametrically of order g -2 . Indeed, in this context, the probability of occurrence of a given fully specified final state is exponentially suppressed, like e -c/g 2 .

Particle production at next-to-Leading order (one loop)

In situations where the particles under consideration couple directly to a time dependent external source, these particles can be produced by the leading order mechanism described in the previous subsections, and the moments of the particle distribution are expressible in terms of the classical field generated by this source. However, when this direct coupling does not exist or when the external source is static, the production of particles is impossible at this order and can at best happen at the next-to-leading order in the coupling. This is in particular the case in the traditional setup of the Schwinger mechanism, where an electromagnetic current couples indirectly to fermion pairs via an electromagnetic field. Therefore, our goal in this subsection is not to calculate the subleading corrections to particle spectra in general, but to calculate the first nonzero contributions they receive when the leading order contribution vanishes.

In this subsection, we discuss the new production mechanisms that arise at NLO, by considering the single particle spectrum at 1-loop. This spectrum is given by the first derivative of the generating functional F[z(p)], for which a general formula was given in eq. [START_REF] Schenke | Fluctuating Glasma initial conditions and flow in heavy ion collisions[END_REF]. At NLO, i.e. at the order g 0 , it involves two quantities [START_REF] Gelis | High energy factorization in nucleus-nucleus collisions[END_REF]:

i. The 1-loop corrections β ± to the 1-point functions ϕ ± , ii. The 2-point function G +-at tree level.
However, since we are interested in the NLO corrections only in cases where the leading order is zero, we do not need to evaluate the terms that contain β ± . Indeed, β ± being a 1-loop correction to ϕ ± , it also leads to a vanishing contribution for static sources or if the particles of interest cannot couple directly to the external source. Therefore, only the second term is needed. This term corresponds to the production of a pair of particles from the external source (in the case of the single particle spectrum, one of the two produced particles is integrated out). Moreover, since we are not going to calculate further derivatives with respect to z(p), it is sufficient to evaluate this quantity at the point z(p) = 1 -which simplifies considerably the calculation. This NLO correction to the single inclusive spectrum reads

dN 1 d 3 p NLO = 1 (2π) 3 2E p d 4 xd 4 y e ip•(x-y) ( x + m 2 )( y + m 2 ) G +-(x, y) . (48) 
The 2-point function G +-(x, y) at tree level obeys the following integral equation :

G (x, y) = G 0 (x, y) -i η=± η d 4 z G 0 η (x, z) U (ϕ(z)) G η (z, y) , (49) 
and therefore mixes with the other three components, G -+ , G ++ and G --. Here, -iU (ϕ(z)) is the general form for the insertion of a background field on a propagator in a theory with a potential U (ϕ). From these equations 9 , one can obtain the following formulas :

G +-(x, y) = d 3 p (2π) 3 2E p a * p (x)a p (y) , G -+ (x, y) = d 3 p (2π) 3 2E p a p (x)a * p (y) , (50) 
where the functions a p (x) are mode functions defined by

x + m 2 + U (ϕ(x)) a p (x) = 0 , lim x 0 →-∞ a p (x) = e ip•x . (51) 
Thus, the problem of finding the Schwinger-Keldysh propagators in a background field can be reduced to determining how plane waves propagate on top of the classical background field (and are distorted by this field). In the previous formulas, we have used the plane wave basis, but other choices are possible, as long as one chooses a properly normalized complete basis. Note also that the 4-dimensional space-time integrals in eq. ( 48) can in fact be rewritten as purely spatial integrals on a constant-t surface (with t → +∞),

dN 1 d 3 p NLO = 1 (2π) 3 2E p d 3 k (2π) 3 2E k lim t→+∞ d 3 x e -ip•x (∂ t -iE p ) a k (x) 2 . ( 52 
)
Before going into more specialized subjects, let us summarize this section by diagrammatic representations of the LO and NLO contributions to the single particle spectrum,

dN 1 d 3 p = p LO + p + p + • • • NLO ( 53 
)
In this representation, the line represented in bold is the on-shell propagator of momentum p corresponding to the final particle one is measuring. In this illustration, we show only one typical graph representative of each class. Therefore, each of the trees that appear in these graphs should be understood as the infinite series that comprises all such trees, and an arbitrary number of such trees can be attached around the loop for the NLO graphs.

Let us also mention a few issues that arise in some special cases:

i. Static sources : The LO contribution is zero unless the external source is time dependent, so that the classical field it generates has nonzero time-like Fourier components. Indeed, this contribution is just the square of the Fourier transform of the classical field generated by the external source. Likewise, the first of the two NLO contributions contributes only when the external source is time dependent.

The second NLO contribution is also zero for a static source if one inserts a finite number of trees around the loop. This means that this 1-loop graph does not contain any contribution analytic in the external field. However, it also contains a non-perturbative contribution that exists even for a static source. This non-perturbative contribution is non-analytic in the coupling constant. More precisely, all its Taylor coefficients at zero coupling vanish, which is why it cannot be seen with any finite order in the external field. This 1-loop diagram contains the contribution that one usually calls the Schwinger mechanism.

ii. Time dependent sources : In situations where the external source is time dependent (especially when it is slowly varying compared to the natural frequency scale set by the mass of the particles), there may be a competition between the perturbative (analytic) contributions and the non-perturbative (non-analytic) ones, and one should avoid considering them in isolation if one wishes to describe the transition between static and time dependent fields. The general formulation that we have adopted in this presentation is well suited to this case. Indeed, since the external field is treated to all orders, it naturally "packages" the perturbative and non-perturbative contributions in the same formulas.

iii. Particles not directly coupled to the external source : The structure illustrated pictorially in eq. ( 53) is completely generic, despite the fact that we have used the example of a scalar field theory in this section. Of course, certain graph topologies may be impossible in certain theories. For instance, in QCD, if the external source is a color current and if we are interested in the quark spectrum, then the first two topologies of eq. ( 53) cannot exist because a single quark field does not couple directly to a color field. For the quark spectrum, the second 1-loop graph therefore constitute the lowest order contribution [START_REF] Gelis | Quark-antiquark production from classical fields in heavy ion collisions: 1+1 dimensions[END_REF][START_REF] Gelis | Chemical thermalization in relativistic heavy ion collisions[END_REF][START_REF] Gelis | Quark production in heavy ion collisions I. Formalism and boost invariant fermionic light-cone mode functions[END_REF]. In contrast, if we are interested in the gluon spectrum, then all the three topologies contribute.

iv. Gauge theories : Even if there is no (chromo)electromagnetic field at asymptotic times, it may happen that the gauge potential is a nonzero pure gauge at t → ±∞.

In the presence of such a pure gauge background, one should not use the vacuum mode functions (i.e. plane waves) in the Fourier decomposition of the fields in order to calculate the particle spectrum. Instead, one should first apply to the mode functions the gauge rotation that transforms the null gauge potentials into the pure gauge of interest.

3 Correlations in the Schwinger mechanism

Generating functional at one loop

In the section 2.7.2, we have observed that when particles can be produced directly at leading order from an external source, their multiplicity distribution is a Poisson distribution at this order. The interpretation of this absence of correlations is that, when particles are produced at leading order by a strong (time dependent) source, they come from 1-point functions and the graphs that produce two or more particles can be factorized into subgraphs in which a single particle is produced. In this case, deviations from a Poisson distribution only arise in subleading corrections. For instance, the second NLO graph of eq. ( 53) contains the production of particle pairs, which introduces correlations among the final state particles. In order to discuss this effect in more detail, let us consider a simple scalar QED model in which we disregard the self interactions of the charged scalar fields [START_REF] Fukushima | Multiparticle correlations in the Schwinger mechanism[END_REF]. Its Lagrangian reads

L ≡ (D µ φ)(D µ φ) * -m 2 φφ * , (54) 
where D µ ≡ ∂ µ -igA µ is the covariant derivative. In this section, the electromagnetic potential A µ is assumed to be non-dynamical. In other words, it is a purely classical field imposed by some external action, and we disregard the feedback (such as screening effects) of the produced charged particles on the electromagnetic field. In this model, the lowest order for the production of charged scalars is at 1-loop, because there is no direct coupling involving A µ 's and a single field φ. Since these particles are charged, it is interesting to keep track separately of particles and antiparticles. For this, we generalize the generating functional introduced in eq.( 25) into

F[z, z] ≡ ∞ m,n=0 1 m!n! m i=1 d 3 p i (2π) 3 2E pi z(p i ) n j=1 d 3 q j (2π) 3 2E qj z(q j ) × p 1 • • • p m particles q 1 • • • q n antiparticles out 0 in 2 , (55) 
where z(p) and z(q) are two independent functions. Unitarity trivially implies

F[1, 1] = 1.
By differentiating with respect to z or z, we obtain the single particle and antiparticle spectra,

dN + 1 d 3 p = δF[z, z] δz(p) z=z=1 , dN - 1 d 3 q = δF[z, z] δz(q) z=z=1 , (56) 
and second derivatives give the two particles and two antiparticles spectra:

dN ++ 1 d 3 p 1 d 3 p 2 = δ 2 F[z, z] δz(p 1 )δz(p 2 )) z=z=1 , dN -- 1 d 3 q 1 d 3 q 2 = δ 2 F[z, z] δz(q 1 )δz(q 2 )) z=z=1 , (57) 
as well as a mixed spectrum

dN +- 1 d 3 pd 3 q = δ 2 F[z, z] δz(p)δz(q)) z=z=1 . ( 58 
)
At the lowest non-zero order (i.e. one loop), it is possible to obtain a compact expression of this generating functional. This will provide complete information about the production of charged particles10 by an external field at this order, and make a contact with other approaches to this problem. From the general result that the generating functional is the sum of the vacuum diagrams in the Schwinger-Keldysh formalism, and the fact that this sum is the exponential of the subset of the connected vacuum diagrams, we can first write:

ln F[z, z] = constant + + - + + - + - + + -+ - + - + ... , (59) 
where the unwritten constant is independent of z and z (its value should be adjusted to satisfy unitarity, i.e. ln F[1, 1] = 0). In eq. ( 59), the objects attached to the loop already resum infinite sequences of + orvertices respectively,

T + ≡ + = + + + + + + + + + ... T -≡ - = - + -- + --- + ... (60) 
Therefore, these objects do not contain any G 0 +-or G 0 -+ propagators, and are therefore independent of z or z. The z and z dependence is carried by the propagators that appear explicitly in the loops in eq. ( 59). In the diagrammatic representation used in eqs. [START_REF] Smolyansky | Dynamical derivation of a quantum kinetic equation for particle production in the Schwinger mechanism[END_REF], the dotted lines are a shorthand for the sum of the two possible interactions that exist in scalar QED:

= + . (61) 
(The wavy lines terminated by a circle denote the external electromagnetic potential.)

The sum in eq. ( 59) can be written in the following compact form,

ln F[z, z] = constant + ∞ n=1 1 n Tr T + (zG 0 +-) T -(zG 0 -+ ) n = constant -Tr ln 1 -T + (zG 0 +-) T -(zG 0 -+ ) , (62) 
where we have made explicit that the factors z and z come along with the off-diagonal propagators G 0 +-and G 0 -+ (the functions z and z carry the same momentum as the propagator they are attached to). The trace denotes an integration over all the spacetime coordinates of the vertices around the loop. The factor 1/n is a symmetry factor, absorbed in the second line in the Taylor expansion of a logarithm. By using the relationship between timeordered and retarded propagators, as well as unitarity, the argument of the logarithm can be rearranged and expressed in terms of a retarded scattering amplitude T R :

ln F[z, z] = constant -Tr ln 1 + M[z, z] , (63) 
where

M[z, z] is a compact notation for M[z, z] p,q ≡ d 4 k (2π) 4 (1 -z(k)z(-q)) T R (p, k) G 0 +-(k) T † R (k, q) G 0 -+ (q) = 2πθ(q 0 )δ(q 2 -m 2 ) d 3 k (2π) 3 2E k (1 -z(-k)z(-q)) T R (p, -k) T * R (q, -k) . (64) 
In these formulas, T R (p, k) is the retarded scalar propagator in the external field, amputated of its external lines, with an incoming momentum k and outgoing momentum p (in the second line, we have changed k → -k, so that k 0 > 0.). The equations ( 63) and ( 64) provide the generating functional, from which one can in principle extract all the information about the distribution of produced particles at one loop. Moreover, the retarded scattering matrix it contains can be computed by solving the wave equation in the background field under consideration. More precisely, it can be expressed in terms of mode functions as follows:

T R (p, -k) = lim x0→+∞ d 3 x e ip•x (∂ x0 -iE p ) a k (x) (D 2 + m 2 ) a k (x) = 0 , lim x0→-∞ a k (x) = e ik•x . ( 65 
)
In other words, in order to obtain T R (p, -k), one should start in the remote past with a negative energy plane wave of momentum k, evolve over the background field until late times, and project it on a positive energy plane wave of momentum p.

General structure of the 1-loop particle correlations

Let us denote by z one of the z(p)'s or one of the z(q)'s. From eq. ( 63), we obtain

δF δz = -F Tr (1 + M) -1 δM δz . ( 66 
)
At z(p) = z(q) ≡ 1, we have M = 0 and F = 1, so that this derivative simplifies into

δF δz z=z=1 = -Tr δM δz z=z=1 . ( 67 
)
This is the general structure of the single (anti)particle spectrum. From eq. ( 66), we can take one more derivative to obtain

δ 2 F δz 1 δz 2 z=z=1 = δF δz 1 δF δz 2 z=z=1 -Tr δ 2 M δz 1 δz 2 - δM δz 1 δM δz 2 z=z=1 . (68) 
From this equation, we see quite generally that the two (anti)particle spectra contain the product of the two corresponding single particle spectra, plus a single trace term that contains the non-trivial correlations. In particular, this extra term encodes possible deviations from a Poisson distribution. One can further differentiate with respect to z or z in order to obtain expressions for higher moments of the particle distribution. Starting with the moment of order 3, an additional simplification arises due to the fact that M is proportional to zz. These successive differentiations lead to expressions in terms of traces of products of first or second derivatives of M. Particles that appear in the same trace are correlated, while they are not correlated if their momenta appear in two different traces. The n-th moment contains a term with a single trace, that correlates all the n particles. This term provides the genuine n-particle correlation, since it is not reducible into smaller clusters.

This formulation can be turned into an algorithm for calculating the single particle spectrum, the 2-particle spectrum, etc.. in the presence of an external field. Indeed, after having discretized space (and, consequently, momentum space), M, δM/δz, δM/δz and δ 2 M/δzδz can be viewed as (very large) matrices, and the evaluation of the right side of eqs. ( 67) and ( 68) is just a matter of linear algebra. In order to compute these building blocks, one would have to first obtain the retarded scattering amplitude T R , which can be done by solving the partial differential equation [START_REF] Levai | Nonperturbative enhancement of heavy quark-pair production in a strong SU(2) color field[END_REF] for each mode k. Although this method can in principle provide all the (or more realistically the first few) moments of the distribution of produced particles, it faces a serious computational difficulty in the general case where the background field is a completely generic function of space-time, with no particular symmetry: computing T R would require computations that scale as the square of the number of lattice spacings (one power comes from the size of the spatial domain in eq. ( 65), and another power from the number of modes k). In the subsection 5.2, we describe a method of statistical sampling that considerably reduces this cost, at the expense of a reduced accuracy due to statistical errors.

Correlations in the case of a homogeneous field

The above procedure could in principle be implemented on a lattice, and would provide the answer for a general background field. In the case of a spatially homogeneous background field, momentum conservation allows us to simplify considerably the structure of the generating functional and to completely uncover the n particle moments. Firstly, the scattering matrix becomes diagonal in momentum:

T R (p, -k) = -2iE p (2π) 3 δ(p + k) β p , (69) 
and all the information about the background field is contained in the coefficients β p . This leads to

M[z, z] p,q = 2πθ(q 0 )δ(q 2 -m 2 ) 2E p (2π) 3 δ(p -q) (1 -z(p)z(-q)) β p 2 . ( 70 
)
Since this object is diagonal in momentum, the same momentum p runs in the loop in the trace of eq. ( 63). This implies that only correlations of particles with the same momentum 11 are possible, or correlations between particles of momentum p and antiparticles of momentum -p (because of the respective arguments of the functions z and z). Using this in eq. ( 63), we obtain

ln F[z, z] = constant -V d 3 k (2π) 3 ln 1 -(z(k)z(-k) -1) f k , (71) 
where we denote

f k ≡ |β k | 2 .
The prefactor V is the overall volume of the system, that results from a (2π) 3 δ(0) in momentum space. By differentiating this formula with respect to z and/or z, we obtain the following results for the 1 and 2 particle spectra 12 :

dN + 1 d 3 p = n p ≡ V (2π) 3 f p dN ++ 2 d 3 p 1 d 3 p 2 - dN + 1 d 3 p 1 dN + 1 d 3 p 2 = δ(p 1 -p 2 ) n p1 f p1 dN +- 2 d 3 pd 3 q - dN + 1 d 3 p dN - 1 d 3 q = δ(p + q) n p (1 + f p ) . (72) 
In the limit where f p 1, the right hand side of the two particle correlation would simplify into a form consistent with a Poisson distribution. In contrast, when the occupation number f p is not small, deviations from a Poisson distribution arise due to Bose-Einstein correlations.

Since there is no correlation except for particles with the same momentum or antiparticles with opposite momentum, one can also derive the probability distribution P k (m, n) to produce m particles of momentum k and n antiparticles of momentum -k,

P k (m, n) = δ m,n 1 1 + f k f k 1 + f k n , (73) 
i.e. a Bose-Einstein distribution. The main difference between such a distribution and a Poisson distribution is the existence of large multiplicity tails, that are due to stimulated emission.

From these semi-explicit formulas, we can also clarify the difference between the vacuum survival probability and the exponential of the particle multiplicity. By evaluating the generating functional at z = z = 0, one obtains the following expression for the vacuum-tovacuum transition probability (i.e. the vacuum survival probability):

P 0 = exp -V d 3 k (2π) 3 ln(1 + f k ) , (74) 
while the total particle multiplicity is

N + = V d 3 k (2π) 3 f k . ( 75 
)
Only when the occupation number is small in all modes, we can expand the logarithm in eq. ( 74) and obtain P 0 ≈ exp(-N + ), but this relationship is not exact in very strong fields. This is a limitation of the methods that give only P 0 (for instance by providing a way to calculate the imaginary part of the effective action in a background field): in general, the knowledge of P 0 is not sufficient to obtain the momentum dependence of the spectrum of produced particles (proportional to f k ).

Constant electrical field in scalar QED

For further reference, let us quote a useful approximation 13 for the occupation number and vacuum survival probability in the case of a constant and spatially homogeneous electrical background field E in the z direction. The occupation number f k is independent of the position x and reads:

f k (x 0 ) ≈ θ(k z ) θ(eEx 0 -k z ) exp - π(m 2 + k 2 ⊥ ) eE . ( 76 
)
The k z dependence is a consequence of the acceleration of the particles by the electrical field after they have been produced (we assume that the electrical charge is positive, e > 0). Therefore, we have also

ln(1 + f k ) ≈ θ(k z ) θ(eEx 0 -k z ) ∞ n=1 (-1) n-1 n exp - nπ(m 2 + k 2 ⊥ ) eE , (77) 
and after insertion in eq. ( 74) we obtain

P 0 ≈ exp - V x 0 8π 3 (eE) 2 ∞ n=1 (-1) n-1 n 2 e -nπm 2 /(eE) . ( 78 
)
Note that if we had neglected the Bose-Einstein correlations, i.e. the higher order terms in f k in the expansion of ln(1 + f k ), we would have obtained instead

P 0 ≈ f k 1 exp - V x 0 8π 3 (eE) 2 e -πm 2 /(eE) . ( 79 
)
We should close this section with some words about practical applications of these formulas. In all experimental situations where such formulas are used, one is still very far from the critical electrical field E c = m 2 /e that would make pair production an event that has a probability of order 1. Consequently, the above discussion about Bose-Einstein correlations and the related issue of reconstructing the particle spectrum from the sole knowledge of the vacuum persistence probability P 0 are mostly academic. The formula ( 79) is a very good approximation in these realistic situations, and from P 0 one can read directly the number of produced pairs.

Equivalent formulations of the Schwinger mechanism

In the previous sections, we have exposed a general formulation of the particle production in the presence of strong external sources/fields based on a resummed perturbation theory approach. There is an equivalent derivation based on the canonical quantization of the field operators, in which particle production is described via Bogoliubov transformations [START_REF] Narozhnyi | The Simplist processes in the pair creating electric field[END_REF][START_REF] Gitman | Processes of Arbitrary Order in Quantum Electrodynamics with a Pair Creating External Field[END_REF][START_REF] Soffel | Stability and decay of the Dirac vacuum in external gauge fields[END_REF][START_REF] Ambjorn | Canonical Quantization in Nonabelian Background Fields. 1[END_REF][START_REF] Gavrilov | Vacuum instability in external fields[END_REF][START_REF] Gavrilov | Density matrix of a quantum field in a particle-creating background[END_REF][START_REF] Tanji | Dynamical view of pair creation in uniform electric and magnetic fields[END_REF]. This section is devoted to a presentation of this alternative approach, as well as a few other related methods that are often used in the literatures. In this section, we use spinor QED in an external classical gauge field in order to illustrate these approaches and their relationships.

Bogoliubov transformations

General case

In the method of canonical quantization of the field operators, one can describe the particle production in classical background fields by a Bogoliubov transformation. A fermion field operator ψ(x) coupled to a classical gauge field A µ (x) obeys the following Dirac equation:

[iγ µ D µ -m] ψ(x) = 0 , (80) 
where

D µ = ∂ µ + ieA µ . (81) 
We assume that there is no electromagnetic field at asymptotic times t → ±∞, and take the null initial condition for the gauge field 14 :

lim t→-∞ A µ (x) = 0 . (82) 
Furthermore, we impose the temporal gauge condition

A 0 = 0 , (83) 
which largely simplifies the description of the time-evolution of the system. After this gauge fixing, there is still a residual invariance under all the spatially dependent gauge transformations. Because the Dirac equation is linear in the field operator, its solution can be expanded in normal modes:

ψ(x) = s=↑,↓ d 3 p (2π) 3 ψ in+ p,s (x) a in p,s + ψ in- p,s (x) b in † p,s , (84) 
where a in p,s and b in p,s are annihilation operators for a particle and an antiparticle of momentum p and spin s, respectively. The ψ in± p,s (x) are c-number solutions of the Dirac equation, analogous to the mode functions introduced in the previous sections. The superscript + anddistinguish the positive energy and negative energy modes, while the superscript 'in' specifies the initial condition that the mode functions satisfy are free ones at t → -∞ :

lim t→-∞ ψ in+ p,s (x) = u(p, s) e -ip•x , lim t→-∞ ψ in- p,s (x) = v(p, s) e +ip•x . (85) 
The momentum space free spinors are normalized by

u(p, s)γ µ u(p, s ) = v(p, s)γ µ v(p, s ) = 2p µ δ s,s . (86) 
With the inner product defined by

(ψ|χ) ≡ d 3 x ψ † (t, x)χ(t, x) , (87) 
the mode functions should be normalized as follows

ψ in+ p,s ψ in+ p ,s = ψ in- p,s ψ in- p ,s = 2E p (2π) 3 δ(p -p ) δ s,s , (88) 
ψ in+ p,s ψ in- p ,s = ψ in- p,s ψ in+ p ,s = 0 . (89) 
One can easily confirm that the initial condition (85) satisfies these conditions, and the inner product is conserved 15 by unitary time evolution. The previous orthonormality conditions are therefore satisfied at arbitrary times for any real background gauge field. By using the orthonormality conditions, we can extract the creation and annihilation operators from the field operator as follows:

a in p,s = 1 2E p ψ in+ p,s ψ , b in † p,s = 1 2E p ψ in- p,s ψ . ( 90 
)
From the canonical anti-commutation relation for the field operators

ψ(t, x), ψ † (t, y) = δ(x -y) , (91) 
the following anti-commutation relation for the creation and annihilation operators can be derived:

a in p,s , a in † p ,s = b in p,s , b in † p ,s = 1 2E p (2π) 3 δ(p -p ) δ s,s . (92) 
Instead of the in-solutions, one could have considered the out-solutions, that satisfy the free boundary condition at t → +∞. Although we have assumed that the electromagnetic field is vanishing at asymptotic times, the gauge field A i (x) can generally be nonzero at t → +∞. This non-zero asymptotic gauge field must be a pure gauge since the field strength is zero. Therefore, the out-solutions must approach free spinors that are gauge rotated:

lim t→+∞ ψ out+ p,s (x) = U † (x) u(p, s) e -ip•x , lim t→+∞ ψ out- p,s (x) = U † (x) v(p, s) e +ip•x , ( 93 
)
15 It is essential that the external gauge potential be real for this property to be true.

where U (x) is a gauge factor defined as

U (t, x) ≡ exp ie x dz•A(t, z) . ( 94 
)
The integration path should be contained in the constant-t plane, and its starting point can be chosen arbitrary (this residual arbitrariness amounts to multiplying the spinors by a constant phase). Note that when the background field is a pure gauge, the gauge link U (t, x) depends only on the endpoints of the line integral, but not on the shape of this path.

The field operator can also be expanded in terms of the out-mode functions,

ψ(x) = s=↑,↓ d 3 p (2π) 3 ψ out+ p,s (x) a out p,s + ψ out- p,s (x) b out † p,s . (95) 
Since the out-mode functions satisfy the same orthonormal condition as the in-mode solutions, the creation and annihilation operators for the out-particles obey the same anticommutation relations as those for the in-particles [START_REF] Tanji | Quark pair creation in color electric fields and effects of magnetic fields[END_REF]. At this point, we have two different definitions of a "particle"; one is based on the inbasis and the other on the out-basis. If a nontrivial background gauge field exists at some point of the evolution of the system, these two definitions are in general different. This difference is nothing but the consequence of the particle production from the vacuum under the influence of the background field. If we assume that the initial state is the vacuum, the spectrum of particles observed at the asymptotic time t → +∞ is represented by the in-vacuum expectation value of the out-particle number operator:

dN s d 3 p = 2E p (2π) 3 0 in |a out † p,s a out p,s |0 in . (96) 
The in-vacuum |0 in is defined by a in p,s |0 in = b in p,s |0 in = 0. In order to calculate this spectrum, we need to find the relationship between the creation and annihilation operators of the in-basis and of the out-basis. By substituting the expansion (84) into

a out p,s = 1 2E p ψ out+ p,s ψ , b out † p,s = 1 2E p ψ out- p,s ψ , (97) 
we obtain the following relationship between the two bases,

a out p,s = 1 2E p s d 3 p (2π) 3 ψ out+ p,s ψ in+ p ,s a in p ,s + ψ out+ p,s ψ in- p ,s b in † p ,s , (98) 
b out † p,s = 1 2E p s d 3 p (2π) 3 ψ out- p,s ψ in+ p ,s a in p ,s + ψ out- p,s ψ in- p ,s b in † p ,s . (99) 
We can modify the spinors by a constant phase (which is irrelevant to physical observables) so that the following relations are fulfilled:

ψ out- -p,s ψ in- -p ,s = ψ out+ p,s ψ in+ p ,s * , (100) 
ψ out- -p,s ψ in+ p ,s = -ψ out+ p,s ψ in- -p ,s * . (101) 
Then, eqs. ( 98) and (99) take the form of a Bogoliubov transformation,

2E p a out p,s = s d 3 p (2π) 3 2E p α(p, s; p , s ) a in p ,s + β(p, s; p , s ) b in † -p ,-s , ( 102 
) 2E p b out † -p,s = s d 3 p (2π) 3 2E p α * (p, s; p , s ) b in † -p ,-s -β * (p, s; p , s ) a in p ,s , (103) 
where the coefficients α and β are defined by

α(p, s; p , s ) ≡ 1 2 E p E p ψ out+ p,s ψ in+ p ,s , (104) 
β(p, s; p , s ) ≡ 1 2 E p E p ψ out+ p,s ψ in- -p ,s . (105) 
These equations ensure that the number of anti-particles having momentum -p equals the number of particles having momentum +p:

0 in |a out † p a out p |0 in = 0 in |b out † -p b out -p |0 in . (106) 
By substituting eq. ( 98) into eq. ( 96), we can express the momentum spectrum of produced particles in terms of the mode functions as follows:

dN s d 3 p = 1 (2π) 3 2E p s d 3 p (2π) 3 2E p ψ out+ p,s ψ in- p ,s 2 . ( 107 
)
If we evaluate the inner product in the right hand side at t → +∞, using the out-mode functions given by eq. ( 93), the spectrum reads

dN s d 3 p = 1 (2π) 3 2E p s d 3 p (2π) 3 2E p lim x 0 →+∞ u † (p, s) d 3 x e -ip•x U (x) ψ in- p ,s (x) 2 . ( 108 
)
This equation is the QED analogue of the scalar formula given in eq. ( 52).

Although the particle number can be defined unambiguously only in the asymptotic region where the background electromagnetic field vanishes, it is informative to define quasiparticles at intermediate times when there is nonzero background field. A time-dependent spectrum can be heuristically defined simply by removing the limit of x 0 → +∞ from eq. ( 108),

dN s d 3 p = 1 (2π) 3 2E p s d 3 p (2π) 3 2E p u † (p, s) d 3 x e -ip•x U (t, x)ψ in- p ,s (t, x) 2 . ( 109 
)
This generalization is equivalent to computing the expectation value of a time-dependent particle number operator:

dN s d 3 p = 2E p (2π) 3 0 in |a † p,s (t)a p,s (t)|0 in , (110) 
where an instantaneous quasi-particle definition is introduced by the expansion

ψ(x) = s=↑,↓ d 3 p (2π) 3 ψ (t)+ p,s (x) a p,s (t) + ψ (t)- p,s (x) b † p,s (t) (111) 
with

ψ (t)+ p,s (x) = U † (t, x) u(p, s) e -ip•x , ψ (t)- p,s (x) = U † (t, x) v(p, s) e +ip•x . ( 112 
)
This definition of a time-dependent spectrum naturally interpolates between the zero particle state at t → -∞ and the final state at t → +∞. At intermediate times when the gauge field is not a pure gauge, the gauge link U (x) can depend on the path chosen to define the line integral. Therefore, one must keep in mind that the spectrum evaluated in a region where the background is not a pure gauge suffers from this unavoidable ambiguity of the particle definition.

Uniform electrical field

Let us now restrict ourselves to a spatially homogeneous electric field which can be given by a gauge field that depends only on time:

A 0 = 0 , A i = A i (t) . (113) 
Since the background gauge field has no spatial dependence, the spatial dependence of the mode functions can be trivially factorized as

ψ in± p,s (t, x) = ψ in± p,s (t) e ±ip•x . ( 114 
)
The gauge factor U (x) is path-independent and it simply reads

U (x) = e -ieA(t)•x . ( 115 
)
The x-integration in the inner products between the free mode functions at the time t and the in-mode functions can be performed analytically, resulting in a delta function of momentum of the form δ(p + eAp ). The Bogoliubov transformation between the in-particles and the quasi-particles at time t reads

a p,s (t) = E p+eA E p s α(t; p, s, s ) a in p+eA,s + β(t; p, s, s ) b in † -p-eA,s , (116) b 
† -p,s (t) = E p+eA E p s α * (t; p, s, s ) b in † -p-eA,s -β * (t, p, s, s ) a in p+eA,s , (117) 
with the following time-dependent Bogoliubov coefficients

α(t; p, s, s ) ≡ e iEpt 2 E p E p+eA u † (p, s) ψin+ p+eA,s (t) , (118) 
β(t; p, s, s ) = e iEpt 2 E p E p+eA u † (p, s) ψin- -p-eA,s (t) . (119) 
Note that the momentum label p of the mode functions ψin± ±p,s (t) is shifted by the gauge field eA because of the insertion of the gauge factor. This shift amounts to changing from the canonical momentum to the kinetic momentum, and it is crucial in order to describe properly the acceleration of the produced particles by the electrical field.

From the time independence of the anti-commutation relations, it follows that the following bilinear combination of the Bogoliubov coefficients is also constant:

σ [α(t; p, s, σ)α * (t; p, s , σ) + β(t; p, s, σ)β * (t; p, s , σ)] = δ s,s . (120) 
In terms of the Bogoliubov coefficients, the momentum spectrum of produced particles can be expressed as

dN s d 3 p = V (2π) 3 s |β(t; p, s, s )| 2 , (121) 
where V = (2π) 3 δ(0) is the volume of the system. Thanks to eq. ( 120), the occupation number f p is always smaller than unity

f p ≡ (2π) 3 V dN s d 3 p = s |β(t; p, s, s )| 2 ≤ 1 , (122) 
as required by the Pauli exclusion principle.

Sauter potential

As an explicit example to illustrate the method of Bogoliubov transformations, we derive the spectrum of particles produced by the Sauter-type pulsed electrical field 16 :

E z (t) = E cosh 2 (t/τ ) , (123) 
where τ stands for the pulse duration. The temporal profile of this field is shown in the figure

t 2τ

Figure 4: The temporal profile of the Sauter electrical field given in eq. ( 123).

4. Since the electrical field of eq. ( 123) vanishes exponentially at t → ±∞, a proper particle definition is available without any ambiguity at asymptotic times. Under the influence of this electrical field, the Dirac equation is analytically solvable [START_REF] Narozhnyi | The Simplist processes in the pair creating electric field[END_REF], and one can compute the spectrum in closed form:

(2π) 3 V dN d 3 p = sinh [π(λ + µ -ν)] sinh [π(λ -µ + ν)] sinh (2πµ) sinh (2πν) , (124) 
where we have defined

µ ≡ τ 2 m 2 + p 2 ⊥ + (p z -2eEτ ) 2 , ( 125 
) ν ≡ τ 2 m 2 + p 2 ⊥ + p 2 z , (126) 
λ ≡ eEτ 2 . ( 127 
)
Note that p in these equations is the physical kinetic momentum, while in some works the spectrum is given in terms of a gauge-dependent canonical momentum.

In the figure 5, the spectrum of eq. ( 124) is plotted as a function of p z , for a fixed transverse mass (m 2 + p 2 ⊥ )/eE = 0.1, and for various values of τ . The spectrum has a width of order 2eEτ in the p z direction. This can be understood as follows: particles are produced with a nearly zero momentum because the spatially homogeneous electrical field carries no momentum. After being produced, they are accelerated by the electrical field following the classical equation of motion p z = t eE z (t )dt eEt. Since the particle production and the acceleration mostly happen in the time interval t ∈ [-τ, +τ ], most of the particles are distributed in the range p z ∈ [0, 2eEτ ]. The figure 6 shows the p z -spectrum for √ eEτ = 4, and for various values of the transverse mass m ⊥ ≡ m 2 + p 2 ⊥ . The factor exp(-πm 2 ⊥ /(eE)) is also indicated by thin black lines, which shows that the peak value of the spectrum agrees well with this exponential factor. In this range of parameters, the exponential m ⊥ -dependence of the spectrum indicates that the particle production is dominated by the non-perturbative Schwinger effect. As we will discuss in the section 7.1, there is another range of parameters where the perturbative particle production is the dominant one.

Quantum kinetic equation

In the previous subsection, we have shown that the particle spectrum can be expressed in terms of the Bogoliubov coefficients, which can be computed by solving the Dirac equation for the mode functions. When the electrical field is uniform, it is possible to derive equations that directly describe the evolution of the Bogoliubov coefficient [START_REF] Popov | e + epair production in an alternating electric-field[END_REF]. Furthermore, if the direction of the electrical field is time-independent, one can derive an equation for the distribution function, which is called quantum kinetic equation [START_REF] Smolyansky | Dynamical derivation of a quantum kinetic equation for particle production in the Schwinger mechanism[END_REF][START_REF] Schmidt | A Quantum kinetic equation for particle production in the Schwinger mechanism[END_REF][START_REF] Kluger | The Quantum Vlasov equation and its Markov limit[END_REF].

Let us derive first the equations for the Bogoliubov coefficients in the presence of an uniform and time-dependent gauge field [START_REF] Schmidt | Determinant calculations using random walk worldline loops[END_REF], whose direction may vary in time. The Bogoliubov coefficients are given by eqs. [START_REF] Popov | Pair production in a variable external field (quasiclassical approximation)[END_REF] and [START_REF] Kleinert | Electron-Positron Pair Production in Spaceor Time-Dependent Electric Fields[END_REF]. Note that the momentum p appearing in these equations is the kinetic momentum, which is related to the canonical momentum q by p = q -eA .

Under the influence of a spatially uniform gauge field such as [START_REF] Schmidt | Determinant calculations using random walk worldline loops[END_REF], the canonical momentum is a constant of motion, while the kinetic momentum depends on time due to the acceleration by the electrical field. When we take the time derivative of the Bogoliubov coefficients, the kinetic momentum p must be regarded as a time dependent quantity such that ∂p ∂t = eE. In order to simplify the notations, let us treat the Bogoliubov coefficients as a 2 × 2 matrix whose indices are the spin indexes (s, s ). By taking the time derivative of the Bogoliubov coefficients given in eqs. [START_REF] Popov | Pair production in a variable external field (quasiclassical approximation)[END_REF] 

where we have defined

(t) ≡ eE •p E p t , (131) 
P (t) ≡ eE z E 2 p -(eE •p) p z 2E 2 p m ⊥ 1 + i eE x p y -eE y p x 2E p m ⊥ σ 3 + i m 2E p m ⊥ (eE x σ 2 + eE y σ 1 ) , (132) 
and

S(t) ≡ m(eE x σ 2 + eE y σ 1 ) 2E p (E p -p z ) + (eE x p y -eE y p x ) σ 3 2E p (E p -p z ) . (133) 
(The σ i are the Pauli matrices.) The matrix S(t) represents the precession of the spin in the electrical field, while the matrix P (t) describes the rate of mixing between the coefficients α and β, which is closely related to the rate of pair production. The detailed form of P (t) and S(t) depends on the choice of a spin basis. Here we have used the free spinors whose spin basis diagonalizes the interaction with a uniform electrical field along the z-direction.

If instead we use spinors whose spin basis is given by the eigenstates of σ 3 in the rest frame, the form of these equations becomes essentially equivalent to those in ref. [START_REF] Popov | e + epair production in an alternating electric-field[END_REF]. Of course, spin-averaged observables do not depend on the choice of the spin basis.

Let us now restrict ourselves to the case where the direction of the electrical field is fixed to be the z-direction, and derive the kinetic equation for the distribution function. In this case, the matrices P (t) and S(t) receive important simplifications, namely P (t) is proportional to the unit matrix and S(t) = 0. Consequently the Bogoliubov coefficients are also proportional to the unit matrix. By taking the time derivative of the distribution function (occupation number) defined by

f p (t) ≡ (2π) 3 V 1 2 s dN s d 3 p = |β(t, p)| 2 , (134) 
one obtains

d f p (t) dt = eE z m ⊥ E 2 p Re (g p (t)) , (135) 
where we have introduced the anomalous distribution:

g p (t) ≡ e -2iEpt α(t, p)β(t, p) . (136) 
Eq. ( 135) must be supplemented by an equation for the anomalous distribution,

d g p (t) dt = -2iE p g p (t) + eE z m ⊥ 2E p [1 -2 f p (t)] . (137) 
The latter equation admits the following formal solution

g p (t) = t dt eE z (t )m ⊥ 2E 2 p (t ) [1 -2 f p (t )] e -i(θ(t)-θ(t )) , (138) 
where we denote

θ(t) ≡ 2 t dτ E p (τ ) . (139) 
The factor 1 -2 f p (t ) encodes the effect of Pauli blocking 17 . By substituting eq. ( 138) into eq. ( 135), we obtain a closed equation for the distribution function itself,

d f p (t) dt = eE z (t)m ⊥ 2E 2 p (t) t dt eE z (t )m ⊥ E 2 p (t ) [1 -2 f p (t )] cos (θ(t) -θ(t )) . ( 140 
)
This equation is called quantum kinetic equation [START_REF] Smolyansky | Dynamical derivation of a quantum kinetic equation for particle production in the Schwinger mechanism[END_REF][START_REF] Schmidt | A Quantum kinetic equation for particle production in the Schwinger mechanism[END_REF][START_REF] Kluger | The Quantum Vlasov equation and its Markov limit[END_REF]. From its derivation, it should be obvious that this formalism is equivalent to solving the Dirac equation for the mode functions, as long as the background electrical field is uniform and its direction is fixed. For practical purposes in numerical calculations, it is easier to solve eqs. ( 135) and ( 137) as associated equations rather than eq. ( 140) which is non-local in time. The non-locality in time of eq. ( 140), obtained after eliminating g p (t) to get a closed equation for f p (t), is reminiscent of the quantum nature of the process under consideration. Similar equations to eqs. ( 135) and (137) for the particle production by parametric resonance in a φ 4 scalar theory can be found in ref. [START_REF] Micha | Turbulent thermalization[END_REF]. In this seemingly unrelated problem, the large zero mode of the field acts as a time dependent background field for the non-zero modes. Because of this analogy, this problem is amenable to a treatment which is very similar to that of the Schwinger mechanism.

Wigner formalism

The Wigner formalism is another approach which has been applied to studies of particle production by the Schwinger mechanism [START_REF] Bialynicki-Birula | Phase space structure of the Dirac vacuum[END_REF][START_REF] Levai | Nonperturbative enhancement of heavy quark-pair production in a strong SU(2) color field[END_REF][START_REF] Hebenstreit | Schwinger pair production in space and time-dependent electric fields: Relating the Wigner formalism to quantum kinetic theory[END_REF][START_REF] Hebenstreit | Strong field effects in laser pulses: the Wigner formalism[END_REF][START_REF] Hebenstreit | Particle self-bunching in the Schwinger effect in spacetime-dependent electric fields[END_REF][START_REF] Berényi | Pair production at the edge of the QED flux tube[END_REF]. This approach shares the same spirit as the quantum kinetic approach in the sense that equations for one-particle distributions are obtained instead of an equation for the elementary fields, though it is more general than the quantum kinetic approach as is applicable to inhomogeneous background fields. Under the influence of inhomogeneous backgrounds, the momentum distribution function is not a well-defined quantity. Instead one can consider the Wigner function defined by

W(t; x, p) ≡ - 1 2 d 3 s e -ip•s 0 in e -ie 1/2 -1/2 dλ A(x+λs)•s ψ(t, x + s/2), ψ(t, x -s/2) 0 in .
(141) A Wilson line factor is inserted between the two field operators in order to preserve the invariance under space-dependent gauge transformations (the temporal gauge condition A 0 = 0 is assumed). Here, we take a straight segment to connect the two points. As already mentioned earlier, if the background is not a pure-gauge potential (i.e. if there are non-zero electrical or magnetic fields), the Wilson line depends on the path one chooses. Therefore, one should keep in mind that a certain amount of arbitrariness is present here (but nothing worse than our earlier attempts to define a particle spectrum in a non-pure gauge background).

Note that since the uncertainty principle does not allow the simultaneous measurement of the position and momentum of a particle, the Wigner distribution W(t; x, p) is not a proper probability distribution. But its integrals over x or p are probability distributions in p or x, respectively. In fact, even the positivity of the Wigner function is not guaranteed. However, it is generally the case that the support in phase-space of the negative values of W is of order . After a coarse graining of phase-space into cells of size or more, these negative regions usually disappear and a probabilistic interpretation becomes plausible.

From the Dirac equation for the fermion field operator, one can derive an evolution equation for the Wigner distribution. If the gauge field is a quantum field operator, the equation is not closed and it depends on higher order correlation functions. Since here we regard the gauge field as a classical background field, a closed equation can be derived:

D t W = - 1 2 D x • γ 0 γ, W -im γ 0 , W -iP • γ 0 γ, W , (142) 
where D t , D x , and P are non-local operators defined by

D t ≡ ∂ t + e 1/2 -1/2 dλ E(t, x + iλ∇ p )•∇ p , (143) 
D x ≡ ∇ x + e 1/2 -1/2 dλ B(t, x + iλ∇ p ) × ∇ p , (144) 
P ≡ p -ie 1/2 -1/2 dλ λ B(t, x + iλ∇ p ) × ∇ p . (145) 
Since eq. ( 142) is derived from the Dirac equation ( 80) without any approximation, this approach is equivalent to directly solving the Dirac equation (e.g. via the mode function method). An advantage of this approach over the mode function method is that eq. ( 142) makes gauge invariance more manifest because it depends on the electric and magnetic fields but not on the gauge field (but keep in mind the caveat mentioned in the paragraph following eq. ( 141)). Although eq. ( 142) is valid for arbitrary space-dependent electromagnetic fields, the numerical implementation of the non-local operators is difficult to achieve. Effects of a small spatial inhomogeneity have been studied with a derivative expansion in ref. [START_REF] Hebenstreit | Schwinger pair production in space and time-dependent electric fields: Relating the Wigner formalism to quantum kinetic theory[END_REF].

When the electric field is uniform with a fixed direction and the magnetic field is absent, eq. ( 142) can be reduced to the quantum kinetic equation ( 140) [START_REF] Hebenstreit | Schwinger pair production in space and time-dependent electric fields: Relating the Wigner formalism to quantum kinetic theory[END_REF].

5 Numerical methods on the lattice

Real-time lattice numerical computations

As discussed in the section 4.1, any observable may in principle be computed once we have obtained the mode functions ψ ± p,s (x) by solving the Dirac equation. This method is applicable to completely general space and time dependent background fields. Although we have mainly discussed the momentum spectrum of the produced particles in the previous sections, other observables like the charged current and the energy-momentum tensor can also be expressed in terms of the mode functions. For example, the vacuum expectation value of a fermion bilinear operator ψ † (x)M ψ(y), where M is a matrix in the spinor and coordinate spaces, can be represented in terms of the mode functions as18 

0 in ψ † (x)M ψ(y) 0 in = s d 3 p (2π) 3 2E p ψ in- † p,s (x) M ψ in- p,s (y) . ( 146 
)
Also the momentum spectrum ( 109) can be expressed in this form. In the following, we will omit the index 'in', since the out-mode functions do not appear.

Our problem is thus reduced to solving the Dirac equation for the mode functions with the initial condition [START_REF] Kluger | Fermion pair production in a strong electric field[END_REF]. For general space and time dependent background fields, it is impossible to solve the Dirac equation analytically, and we therefore need to resort to numerical computations. In this subsection, we briefly explain a possible lattice setup for this numerical implementation in SU(N c ) gauge theory [START_REF] Aarts | Real time dynamics with fermions on a lattice[END_REF][START_REF] Aarts | Particle production and effective thermalization in inhomogeneous mean field theory[END_REF]. We assume the temporal gauge condition A 0 = 0, and treat the time variable t as a continuum variable. We divide the 3-dimensional volume, that we take of size of L x × L y × L z , into N x × N y × N z lattice sites. The space coordinates are labeled by integers n i (i = x, y, z) as follows (x, y, z) = (n x a x , n y a y , n z a z ) ,

where the numbers a i = L i /N i are the lattice spacings. It is common practice to use periodic boundary condition (in space) for the fields.

On the lattice, it is more convenient to consider the link variables

U i (x) ≡ e iaigAi(x) , (148) 
that are Wilson lines spanning one elementary edge of the lattice, as the fundamental variables, instead of the gauge fields A i (x). After the gauge fixing by the temporal gauge condition, there is a residual invariance under gauge transformations that depend only on the spatial coordinates. It is highly desirable to preserve exactly this residual invariance through the discretization. Under such a gauge transformation, the fermion fields and the link variables are transformed as

ψ(t, x) -→ Ω(x) ψ(t, x) , (149) 
and

U i (t, x) -→ Ω(x) U i (t, x) Ω † (x + î) . (150) 
Therefore, a natural definition for the covariant derivative applied to the fermion field reads

D i ψ(x) ≡ 1 2a i U i (x) ψ(x + î) -U † i (x -î) ψ(x -î) , (151) 
where we have used a centered difference in order to preserve the unitarity of the theory. This covariant derivative transforms as expected under the residual gauge transformations. With this definition, the Dirac equation on the lattice reads

iγ 0 ∂ 0 + iγ i D i -m ψ(x) = 0 . ( 152 
)
If we regard the background gauge fields as generated from some given sources, we need to solve the Yang-Mills equations on the lattice in addition to the Dirac equation. The lattice Yang-Mills equation can be derived from the lattice Hamiltonian for the SU(N c ) gauge fields,

H gauge = 3 i=1 tr E i E i + N c g 2 3 j=1 (j =i) 1 a 2 i a 2 j 1 - 1 N c Re tr U i,j (x) , (153) 
where the E i (x) are the electrical fields, and the variables U ij (x), called plaquettes, are defined by

U i,j (x) ≡ U i (x) U j (x + î) U † i (x + ĵ) U † j (x) . (154) 
(Plaquettes are Wilson loops spanning an elementary square on the lattice.) The Hamilton equations then read

∂ t U i (x) = iga i E i (x) U i (x) , (155) 
and

∂ t E i (x) = - 1 g a i 3 j=1 (j =i) 1 a 2 j Im [ U i,j (x) + U i,-j (x)] -(trace) , (156) 
where

U i,-j (x) ≡ U i (x) U † j (x + î -ĵ) U † i (x -ĵ) U j (x -ĵ) (157) 
and where, for an element X of the fundamental representation of the SU(N c ) algebra, the notation "X-(trace)" means

X -(trace) = X - 1 N c tr X . (158) 
For a U(1) theory such as QED, the -(trace) term must of course be ignored. Because we employ the centered difference in eq. ( 151), the lattice momentum for the fermion fields suffers from the problem of doublers. In d-dimensional momentum space, there is a 2 d -fold degeneracy. If the background field is spatially homogeneous, it does not carry any nonzero momenta. Therefore, fermions in different momentum modes do not interact with each other, and thus all these degenerated modes remain independent. In this case, we can simply eliminate the effect of the doublers by dividing expectation values by the number of degenerated modes, as long as all the degenerated modes contribute to the expectation value equally 19 . However, when the background gauge field is inhomogeneous, it carries nonzero momenta that cause interactions between fermions in different momentum modes. In this case, a non-trivial mixing happens between degenerate doubler modes, and the fermionic observables may be contaminated by the doublers. One way to suppress the doubler modes is to add the Wilson term to the Dirac equation, which now reads

iγ 0 ∂ 0 + iγ i D i -m ψ(x) + r 2 3 i=1 1 a i U i (x) ψ(x + î) -2ψ(x) + U † i (x -î) ψ(x -î) = 0 , (159) 
where r is a positive parameter. Thanks to the addition of the Wilson term, the doubler modes acquire a heavy mass of the order of the lattice ultraviolet cutoff 1/a, and therefore they decouple from the dynamics when a → 0. Note that the numerical results presented in the subsections 5.3 and 7.2 are computed without the Wilson term since the background gauge fields considered there are uniform.

Statistical sampling method

Formulation

In the previous subsection, we have presented a lattice formulation of the Dirac and Yang-Mills equations. Each mode function ψ ± p,s (x) obeys the Dirac equation

iγ 0 ∂ 0 + iγ i D i -m ψ ± p,s (x) = 0 . (160) 
This equation must be solved for each momentum modes (and other quantum numbers as well). Therefore, the numerical cost to solve the equation is proportional to N 2 latt N t , where N latt ≡ N x N y N z is the number of lattice sites and N t is the number of time steps. This cost has a very unfavorable scaling with the size of the lattice, especially in 3 dimensional space. If the system has some spatial symmetry, this rather expensive cost can be reduced. For instance, if the system is completely uniform, the x-dependence of the mode functions is analytically known to be e ip•x . In this case, we do not need to treat the space-dependence of the field numerically, and thus the cost is reduced to N latt N t .

In the absence of any such symmetries, a way to reduce the numerical cost is to replace the exhaustive listing of the momentum modes by a Monte Carlo sampling method, which is applicable to any space and time dependent background. However, instead of picking randomly a subset of the set of momenta, one can exploit the linearity of the Dirac equation and use random spinors that are linear superpositions (with random coefficients) of all the momentum modes [START_REF] Gelis | Quark production in heavy ion collisions I. Formalism and boost invariant fermionic light-cone mode functions[END_REF][START_REF] Gelis | Formulation of the Schwinger mechanism in classical statistical field theory[END_REF]. This approach provides a better sampling of the entire momentum space.

Let us introduce a stochastic field which is a linear superposition of the mode functions with random number coefficients:

ψ - c (x) = 1 √ V k,s 1 √ 2E k c k,s ψ - k,s (x) , (161) 
where the random numbers are Gaussian distributed with zero mean value and the following variance:

c k,s c * k ,s ens = δ k,k δ s,s . (162) 
• • • ens denotes the ensemble average over these random numbers. Using this stochastic field, we can express the expectation value of a bilinear operator [START_REF] Nayak | Non-perturbative quark-antiquark production from a constant chromo-electric field via the Schwinger mechanism[END_REF] as follows:

0 in ψ † (x) M ψ(y) 0 in = ψ - † c (x) M ψ - c (y) ens . ( 163 
)
Because the Dirac equation is linear, the stochastic fields (161) obey the same Dirac equation as the mode functions. Therefore, expectation values at time t can be computed by the following procedure: i. At an initial time t = t 0 , draw random numbers following the variance [START_REF] Gavrilov | Dirac fermions in strong electric field and quantum transport in graphene[END_REF], and compute a stochastic field (161) by using the initial value for the mode functions ψ - p,s (t 0 , x). In most practical cases, the mode functions at t 0 are known analytically.

ii. Solve the Dirac equation for the stochastic field until the time t, and compute the contribution of this field to the observable.

iii. Repeat the steps i and ii in to perform the ensemble average until satisfactory statistical errors are reached.

The numerical cost of the step i is proportional to N 2 latt since the sum over the momentum modes must be done at each spatial coordinate 20 . This cost scales unfavorably with N latt . However, this must be computed only at the initial time. The cost for the time evolution is lower than that of the direct mode function method, since it is proportional to N latt N t for each configuration of the stochastic field. If we compute N conf configurations for the stochastic field, the total cost is roughly proportional to (N 2 latt + N latt N t )N conf . This cost is lower than that of the direct mode function method provided that

N conf N latt and N conf N t ,
both of which are easily satisfied in practice.

Statistical errors

Eq. ( 163) is exact only if the number of configurations N conf is infinite. In practice, we have to use a finite N conf , which implies nonzero statistical errors. In this subsection, we describe a method to evaluate the statistical errors [START_REF] Gelis | Quark production in heavy ion collisions I. Formalism and boost invariant fermionic light-cone mode functions[END_REF]. We consider the following vacuum expectation:

0 in Ô 0 in ≡ X,f 0 in ψ † (x) M ψ(y) 0 in , (164) 
where the sum over X denotes possible summation over x and y coordinates (for instance in a Fourier transform to compute a spectrum), and the sum over f is over some other quantum numbers that the matrix M may carry. For instance, in the case of the momentum spectrum [START_REF] Gies | Quantum diffusion of magnetic fields in a numerical worldline approach[END_REF], the matrix M is

M spec = 1 (2π) 3 2E p u(p, s)u † (p, s) U † (t, x) U (t, y) e -ip•(x-y) , (165) 
and X stands for {x, y}. If we compute the spin-averaged spectrum, f would be the final spin s.

With a finite number N conf of configurations, the statistical evaluation of the expectation value ( 164) can be expressed as

O N conf = 1 V X,f J,J 1 2 √ E J E J C N conf (J, J ) ψ - † J (x) M ψ - J (y) , ( 166 
)
where J is a shorthand notation for (k, s). The coefficient C N conf contains the average over N conf samples of the Gaussian random numbers c

(n) J :

C N conf (J, J ) ≡ 1 N conf N conf n=1 c (n) * J c (n) J . ( 167 
)
C N conf itself is a random number, whose fluctuations determine the statistical error of the evaluation [START_REF] Fukushima | Schwinger Mechanism with Stochastic Quantization[END_REF]. From the variance of the random coefficients c

(n)

J , we obtain immediately the following mean value for C N conf (J, J )

C N conf (J, J ) = δ J,J . (168) 
Furthermore, if we assume for simplicity that the c

(n) J
are Gaussian distributed (with no correlations if n = n ), we obtain the following variance

C N conf (J, J )C N conf (K, K ) -C N conf (J, J ) C N conf (K, K ) = 1 N conf δ J,K δ K,J , (169) 
which decreases as one increases the number of configurations. By using this equation, we can compute the variance of the statistical evaluation (166) by

O 2 N conf -O N conf 2 = 1 N conf 1 V 2 J,J 1 4E J E J X,f ψ - † J (x) M ψ - J (y) 2 . ( 170 
)
The square root of the right hand side provides an estimate of the statistical error. We can compute it by the statistical method in the following procedure:

i. Prepare two stochastic fields

ψ - 1,2 (x) = 1 √ V J 1 √ 2E J c (1,2) J ψ - J (x) , (171) 
with uncorrelated random numbers c

(1)

J and c

(2) J . ii. Evolve these two stochastic fields by solving their Dirac equation.

iii. Compute

X,f ψ - † 1 (x) M ψ - 2 (y) 2 . ( 172 
)
iv. Repeat the steps i-iii in order to average over the random numbers c

(1) J and c

(2)

J . Because this is merely for an error estimate, one does not need a large number of samples.

v. At the last step, take the square root and divide the result by √ N conf .

Since the summand of the sum over X, f in eq. ( 170) is generally a complex number, phase cancellations may happen in the summation over X and f , making the statistical error smaller. For example, the transverse spectrum, which is obtained by summing over f ≡ (p z , s), should contain smaller statistical error compared with the full momentum spectrum that depends on p z and p ⊥ .

Relation to the low-cost fermion method

In real-time lattice simulations for fermion fields, another method called the low-cost fermion method [START_REF] Borsanyi | Low-cost fermions in classical field simulations[END_REF] has been used in several works [START_REF] Saffin | Real-time Fermions for Baryogenesis Simulations[END_REF][START_REF] Saffin | Dynamical simulations of electroweak baryogenesis with fermions[END_REF][START_REF] Berges | Amplified Fermion Production from Overpopulated Bose Fields[END_REF][START_REF] Kasper | Fermion production from real-time lattice gauge theory in the classical-statistical regime[END_REF]. In the low-cost fermion method, instead of using one stochastic field [START_REF] Zubkov | Schwinger pair creation in multilayer graphene[END_REF], one employs two kinds of stochastic fields called "male" and "female" fields:

ψ M (x) ≡ 1 √ 2V k,s 1 √ 2E k c k,s ψ + k,s (x) + d k,s ψ - k,s (x) , (173) 
ψ F (x) ≡ 1 √ 2V k,s 1 √ 2E k c k,s ψ + k,s (x) -d k,s ψ - k,s (x) , (174) 
where c k,s and d k,s are independent random numbers which have the same variance as [START_REF] Gavrilov | Dirac fermions in strong electric field and quantum transport in graphene[END_REF]. Combining these two fields, one can obtain the vacuum expectation value of a symmetrized bilinear operator by an ensemble average as follows:

0 in 1 2 ψ † (x), M ψ(y) 0 in = 1 2V k,s 1 2E k ψ - † k,s (x) M ψ - k,s (y) -ψ + † k,s (x) M ψ + k,s (y) = -ψ † M (x) M ψ F (y) ens . (175) 
The two kinds of fields are necessary in order to obtain the minus sign in front of the second term in the right hand side of the first line, which originates from the anti-commutation relation for the fermionic operators. For bosonic fields, one would need only one kind of stochastic field [START_REF] Gelis | Formulation of the Schwinger mechanism in classical statistical field theory[END_REF]. By using the completeness relation

1 V k,s 1 2E k ψ + k,s (t, x)ψ + † k,s (t, y) + ψ - k,s (t, x)ψ - † k,s (t, y) = V δ x,y 1l (176) 
(1l is the unit matrix in the spinor space), we can relate the quantities evaluated in the simple statistical method [START_REF] Klimchitskaya | Creation of quasiparticles in graphene by a time-dependent electric field[END_REF] and that in the low-cost fermion method (175) by

ψ - † c (x) M ψ - c (y) ens = -ψ † M (x) M ψ F (y) ens + V 2 tr (M δ x,y ) . (177) 
Therefore, the two methods provide the same result for the vacuum expectation value up to a zero-point vacuum contribution. However, the method ( 163) has an advantage over the low-cost fermion method: it leads to smaller statistical errors for the evaluation of the spectrum if the value of the spectrum is much smaller than one. In the statistical method ( 163), the spectrum is directly obtained by the statistical ensemble. On the other hand, in the low-cost fermion method [START_REF] Hashimoto | Electromagnetic instability in holographic QCD[END_REF], one gets instead a direct access to 1 2f p (f p being the fermion occupation number), and the zero-point occupation number 1/2 must be subtracted. Because this vacuum 1/2 also contains statistical errors, the low-cost fermion method suffers from relatively larger statistical errors when the occupation number is much smaller than 1/2.

Numerical example

To demonstrate the efficiency of the statistical method, we use it in order to compute the momentum spectrum of particles produced in the Sauter electrical field [START_REF] Dunne | Multidimensional Worldline Instantons[END_REF]. For this background field, one can use the analytic expression [START_REF] Dumlu | Complex Worldline Instantons and Quantum Interference in Vacuum Pair Production[END_REF] to monitor the accuracy of the numerical evaluation. In the figure 7, the p z -spectrum computed by the statistical method based on eq. ( 161) is compared with the analytical result. The lattice parameters used for this computation are 21 In the figure 8, the spectrum computed by the statistical method ( 161) is compared to the result of the low-cost fermion method [START_REF] Hashimoto | Electromagnetic instability in holographic QCD[END_REF]. Firstly, we see that both methods agree well with the exact result. As we have pointed out in the section 5.2.3, the statistical errors for the single-field stochastic method are smaller than those of the male-female method in the region where the spectrum is smaller than 1/2. This difference becomes crucial when we compute the particle production in the weak field (or high momentum) regime such that eE < m 2 ⊥ . In such weak fields, the occupation number of produced particle is much smaller than one. In order to resolve these small occupation numbers, we need to ensure that the statistical errors are smaller than the occupation number itself. This is easily attainable if we use the single stochastic field method. In the figure 9, we display the spectrum computed by the single stochastic field method in a weak field eE = 0.25 m 2 , √ eE τ = 25.5. The lattice parameters are N z = 256, ma z = 0.048, and N conf = 48. The occupation number of the order of 10 -6 that one obtains with these parameters is accurately reproduced within statistical errors. If one had used the male-female method, an extremely large number of configurations of the order of 10 12 would have been necessary in order to achieve a similar statistical accuracy, because large statistical errors arise from the vacuum 1/2. 21 In fact, the values of the transverse (x and y) lattice parameters are irrelevant for the pz-spectrum with fixed transverse momentum p ⊥ = p 2

x + p 2 y , since there is no dynamics in the direction transverse to the electrical field [START_REF] Dunne | Multidimensional Worldline Instantons[END_REF]. 

Back reaction

In the previous subsection, we have shown the momentum spectra of particles produced in a gauge field whose evolution is controlled by a given external current. At the level of accuracy at which this calculation was done, the total energy of the system (the electromagnetic field and the produced particles) is not conserved. In order to conserve the total energy, we must take into account the back reaction of the produced particles on the gauge field. When a particle-antiparticle pair is produced, the gauge field is weakened because it looses the energy that the pair takes away and eventually decays [START_REF] Gavrilov | One-loop energy-momentum tensor in QED with electric-like background[END_REF][START_REF] Gavrilov | Consistency restrictions on maximal electric field strength in QFT[END_REF]. In other words, when a particleantiparticle pair is produced, the background gauge field is screened by the produced charges. This screening effect is crucial for the description of string breaking in the color flux tube model [START_REF] Casher | Vacuum polarization and the absence of free quarks[END_REF][START_REF] Glendenning | Creation of anti-Q Q pair in a chromoelectric flux tube[END_REF][START_REF] Hebenstreit | Real-time dynamics of string breaking[END_REF]. In addition to static charges, the screening can also be caused by moving charges (conduction current) [START_REF] Gatoff | The Flux Tube Model for Ultrarelativistic Heavy Ion Collisions: Electrohydrodynamics of a Quark Gluon Plasma[END_REF]. These effects can be taken into account if the gauge field is coupled to a current induced the by produced particles:

∂ µ F µν (x) = Ĵν (x) , (178) 
where Ĵν (x) is charged current operator, and • • • denotes the quantum expectation value by an initial state. In the case of the lattice formulation, the current must be added to the right hand side of eq. ( 156).

The description of the back reaction by eq. ( 178) can be justified in the mean field approximation based on the large-N f expansion [START_REF] Cooper | Quantum Back Reaction in Scalar QED as an Initial Value Problem[END_REF], or by the classical-statistical approximation for strong gauge fields [START_REF] Kasper | Fermion production from real-time lattice gauge theory in the classical-statistical regime[END_REF]. Numerical computations including the back reaction problem have been conducted in several studies [START_REF] Tanji | Dynamical view of pair creation in uniform electric and magnetic fields[END_REF][START_REF] Kluger | The Quantum Vlasov equation and its Markov limit[END_REF][START_REF] Kasper | Fermion production from real-time lattice gauge theory in the classical-statistical regime[END_REF][START_REF] Kluger | Pair production in a strong electric field[END_REF][START_REF] Kluger | Fermion pair production in a strong electric field[END_REF][START_REF] Kluger | Pair production in a strong electric field: An Initial value problem in quantum field theory[END_REF][START_REF] Cooper | Particle production in the central rapidity region[END_REF][START_REF] Dawson | Fermion particle production in semi-classical Boltzmann-Vlasov transport theory[END_REF][START_REF] Asakawa | Dilepton production from a nonequilibrium quarkgluon plasma in ultrarelativistic nucleus-nucleus collisions[END_REF][START_REF] Vinnik | Plasma production and thermalization in a strong field[END_REF][START_REF] Bloch | Pair creation: Back reactions and damping[END_REF][START_REF] Tanji | Quark pair creation in color electric fields and effects of magnetic fields[END_REF][START_REF] Hebenstreit | Simulating fermion production in 1+1 dimensional QED[END_REF]. 

0 ) 3/2 √ eE 0 t E z (t) J z (t)
Figure 10: The time evolution of the electrical field E z (t) scaled by its initial value E 0 (solid line, left axis), and the induced charge current J z (t) divided by (eE 0 ) 3/2 (dotted line, right axis). The parameters are e = 0.3, m/

√ eE 0 = 0.1, N x = N y = 48, N z = 512, √ eE 0 a x = √ eE 0 a y = 0.62, √ eE 0 a z = 0.029.
As an illustration of the effects of back reaction, we show numerical results obtained in real-time lattice computations. We have solved (the lattice version of) the Maxwell equation [START_REF] Kim | Schwinger Pair Production in dS(2) and AdS(2)[END_REF] with a uniform electrical field, E z (0) = E 0 , as initial condition. For the fermionic sector, the mode functions are directly computed by solving the Dirac equation, and the charge current appearing in the right hand side of eq. ( 178) is calculated at every time step. The parameters used in the computation are e = 0.3, m/ √ eE 0 = 0.1, N x = N y = 48, N z = 512, √ eE 0 a x = √ eE 0 a y = 0.62, √ eE 0 a z = 0.029. In the figure 10, the time evolution of the electrical field E z (t) and the induced current density J z (t) are shown. As particles are produced and accelerated by the electrical field, a positive charge current is induced, causing the reduction of the electrical field strength through the Maxwell equation [START_REF] Kim | Schwinger Pair Production in dS(2) and AdS(2)[END_REF]. At some point ( √ eE 0 t 45), the field strength vanishes. However, since particles are already accelerated to high momenta, the decrease of the electrical field does not stop and its direction is flipped. Then, the particles start to be accelerated in the opposite direction. The repeat of such processes results in plasma oscillations. Thanks to the back reaction, the total energy is indeed conserved (see the figure 11). The p z -spectrum of produced particles is shown in the figure 12. In the top figure, the spectrum at early times before the first minimum of the electrical field happens is plotted. In this time range, the behavior of the spectrum is very similar to that obtained in the Sauter electrical field, that we have shown in the section 4.1.3: Particles are produced with a nearly zero longitudinal momentum and then accelerated by the electrical field. Their occupation number is given by exp(-πm 2 ⊥ /(eE)). After the direction of the electrical field is flipped, non-trivial quantum interference phenomena occur. Because the direction of the electrical field is negative from √ eE 0 t ∼ 60 to ∼ 100, the produced particles are accelerated in the negative z direction, and thus the p z -spectrum is shifted towards the negative momentum region. While doing this, they cross the zero momentum point p z = 0, where particles are still being produced by the electrical field (for a spatially homogeneous field, the particles are created at nearly zero momentum). Pauli blocking and interferences with the pre-existing particles lead to a spectrum that displays distinct features in the region p z < 0: (i) its magnitude is significantly smaller than exp(-πm 2 ⊥ /(eE)) and (ii) it shows rapid oscillation in p z (see ref. [START_REF] Tanji | Quark pair creation in color electric fields and effects of magnetic fields[END_REF]). Because of this interference phenomenon that repeats with each oscillation of the field, the spectrum at √ eE 0 t = 200 is not smooth, in contrast to the spectrum at early times shown in the top figure. Even though this system is closed and the equations of motion are symmetric under time-reversal, an apparent irreversibility of the energy flow from coherent fields to fluctuating quantum modes emerges [START_REF] Habib | Dissipation and decoherence in mean field theory[END_REF]. In the figure 11, we can indeed see that the oscillation amplitude of the electrical field is slightly decreasing in time. 

Worldline formalism

Up to this point, our discussion of the Schwinger mechanism was following a fairly standard quantum field theoretical approach based on diagrammatic expansions and field equations of motion. The non-perturbative nature of the phenomenon resulted in the necessity of summing infinite sets of Feynman graphs, the sum of which may be re-expressed in terms of solutions of some partial differential equations (classical field equations of motion at Leading Order, or the equation of motion for the mode functions over a background field at one-loop).

In the past ten years, a radically different approach -the worldline formalism-has been applied to the study of the Schwinger mechanism. The worldline formalism first emerged from ideas around the limit of infinite string tension in string theory [START_REF] Metsaev | On loop corrections to string theory effective actions[END_REF][START_REF] Bern | A New Approach to One Loop Calculations in Gauge Theories[END_REF][START_REF] Bern | Efficient calculation of one loop QCD amplitudes[END_REF] but its potential for performing efficiently 1-loop calculations in quantum field theories was soon realized. This formalism was later establish of purely field theoretical grounds in ref. [START_REF] Strassler | Field theory without Feynman diagrams: One loop effective actions[END_REF], in an approach reminiscent of Schwinger's proper time method. Reviews on this formalism can be found in refs. [START_REF] Schubert | An Introduction to the worldline technique for quantum field theory calculations[END_REF][START_REF] Schubert | Perturbative quantum field theory in the string inspired formalism[END_REF]. Its main application, of direct interest to us in this review, if the calculation of 1-loop effective actions in a background field [START_REF] Reuter | Constant external fields in gauge theory and the spin 0, 1/2, 1 path integrals[END_REF][START_REF] Schmidt | On the calculation of effective actions by string methods[END_REF][START_REF] Schmidt | Multiloop calculations in the string inspired formalism: The Single spinor loop in QED[END_REF][START_REF] Gies | Quantum effective actions from nonperturbative worldline dynamics[END_REF], with special emphasis on pair creation in refs. [START_REF] Dunne | Worldline instantons and pair production in inhomogeneous fields[END_REF][START_REF] Dunne | Worldline instantons. II. The Fluctuation prefactor[END_REF], but is has also been used in studies of the Casimir effect [START_REF] Gies | Casimir effect on the worldline[END_REF][START_REF] Gies | Worldline algorithms for Casimir configurations[END_REF].

Although the result for the Schwinger mechanism should not depend on the method used to obtain it, the worldline formalism offers very interesting new insights about the space-time development of the particle production process, and provides completely novel numerical approaches to these calculations. In this section, we will present several aspects of the worldline approach by returning to the simple example of scalar QED in an external field for illustration purposes.

Worldline formalism for pair production

As we have discussed in the section 2, the vacuum to vacuum transition amplitude can be written as an exponential,

0 out 0 in = e i V , (179) 
where i V is the sum of all the connected vacuum diagrams. The possibility of particle production is intimately related to the imaginary part of V, since the total probability of producing particles reads

∞ n=1 P n = 1 -P 0 = 1 -e -2 Im V . (180) 
In scalar QED, the graphs made of one scalar loop embedded in a background electromagnetic field lead to the following contribution to V,

V 1 loop = ln det g µν D µ D ν + m 2 , (181) 
where D µ is the covariant derivative in the background field and g µν is Minkowski's metric tensor. For technical reasons related to the definition of the worldline formalism, we should introduce the Euclidean analogue of this quantity,

V E ,1 loop ≡ ln det -D i D i + m 2 = Tr ln -D i D i + m 2 , (182) 
where the index i runs over the values 1, 2, 3, 4 (4 being the Euclidean time). The worldline formalism can be derived from several points of view. A simple starting point is Schwinger's proper time representation of inverse propagators,

-D i D i + m 2 -1 = ∞ 0 dτ exp -τ -D i D i + m 2 , (183) 
that leads to the following formula for the logarithm ln

-D i D i + m 2 = - ∞ 0 dτ τ exp -τ -D i D i + m 2 , (184) 
up to some irrelevant integration constant. Using standard manipulations, the trace of the exponential can then be rewritten as a path integral over closed loops in Euclidean spacetime:

V E ,1 loop = - ∞ 0 dτ τ e -m 2 τ x i (0)=x i (τ ) Dx i (τ ) exp - τ 0 dτ ẋi ẋi 4 + ie ẋi A i (x) .
(185) This formula involves a double integration: a path integral over all the worldlines x i (τ ), i.e. closed paths in Euclidean space-time parameterized by the fictitious time τ ∈ [0, τ ], and an ordinary integral over the length τ of these paths. The dot denotes a derivative with respect to the fictitious time. The sum over all the worldlines can be viewed as a materialization of the quantum fluctuations in space-time, and the prefactor exp(-m 2 τ ) suppresses the very long worldlines that explore regions of space-time that are much larger than the Compton wavelength of the particles (making obvious the role of the mass as an infrared regulator). In contrast, the ultraviolet properties of the theory under consideration are controlled by the short worldlines, i.e. the limit τ → 0. Note that going to Euclidean space-time was necessary in order to obtain convergent integrals in eq. ( 185). In the vacuum, i.e. when there is no background field, one can use the following standard result in d space-time dimensions,

x i (0)=x i (τ ) Dx i (τ ) exp - τ 0 dτ ẋi ẋi 4 = 1 (4πτ ) d/2 = d=4 1 (4πτ ) 2 . (186) 
In eq. ( 185), the path integral can be factored into an integral over the barycenter X i of the worldline and the position r i (τ ) about this barycenter,

x i (τ ) ≡ X i + r i (τ ) , τ 0 dτ r i (τ ) = 0 . (187) 
After this separation, all the information about the background field contained in eq. ( 185) comes via a Wilson line,

W X r(τ ) ≡ exp -ie τ 0 dτ ṙi (τ )A i (X + r(τ )) , (188) 
averaged over all closed loop of length τ ,

W X τ ≡ r i (0)=r i (τ ) Dr i (τ ) W X r(τ ) exp - τ 0 dτ ṙi ṙi 4 x i (0)=x i (τ ) Dx i (τ ) exp - τ 0 dτ ẋi ẋi 4 . (189) 
In this path average, the exponential containing the integral of the squared velocity tends to suppress the long paths. Therefore, the average is dominated by an ensemble of loops (sometimes called a loop cloud ) localized around the barycenter X i , and W X τ encapsulates the local properties of the quantum field theory in the vicinity of X i (roughly up to a distance of order τ 1/2 ). In terms of this averaged Wilson loop, the 1-loop Euclidean connected vacuum amplitude reads

V E ,1 loop = - 1 (4π) 2 d 4 X ∞ 0 dτ τ 3 e -m 2 τ W X τ . (190) 
(In this formula, the prefactor and power of τ in the measure assume 4 spacetime dimensions.) Note that this expression should in principle be supplemented by some 1-loop counterterms determined by the renormalization conditions of the parameters of the bare Lagrangian. The imaginary part of V E ,1 loop comes from poles located at real values of the fictitious time τ , and can be written as

Im (V E ,1 loop ) = π (4π) 2 d 4 X Re poles τn e -m 2 τn τ 3 n Res W X τn , τ n . (191) 

Constant electrical field

It is instructive to reconsider first the well known case of a constant electrical field E. Since one can choose a gauge potential which is linear in the coordinates X i , r i , the path integral that gives the average Wilson loop is Gaussian and can therefore be performed in closed form, leading to

W X τ = eEτ sin(eEτ ) . (192) 
This expression has an infinite series of single poles along the positive real axis, located at

τ n = nπ/(eE) (n = 1, 2, 3, • • • )
, that give the following expression for the imaginary part:

Im (V E ,1 loop ) = V 4 16π 3 (eE) 2 ∞ n=1 (-1) n-1 n 2 e -nπm 2 /(eE) . (193) 
In this formula, V 4 is the volume in space-time over which the integration over the barycenter X i is carried out. This formula is identical to the standard result for the vacuum survival probability P 0 = exp(-2 Im V), that we have already recalled in eq. [START_REF] Gavrilov | One-loop energy-momentum tensor in QED with electric-like background[END_REF]. By comparing the origin of the term of index n in the subsection 3.4 and in the present derivation in the worldline formalism, one sees that the Bose-Einstein correlations (i.e. higher orders in the occupation number) are encoded in the poles τ n that are more distant from τ = 0, while the first pole τ 1 only contains information about the uncorrelated part of the spectrum. This observation gives some substance to the intuitive image that quantum fluctuations and correlations are encoded in the fact that the worldlines explore an extended region around the base point X i . From this correspondence, we see that increasingly intricate (the index n is the number of correlated particles) quantum correlations come from worldlines that explore larger and larger portions of space-time.

Numerical approaches

The worldline formalism can be implemented numerically [START_REF] Gies | Quantum diffusion of magnetic fields in a numerical worldline approach[END_REF][START_REF] Gies | Loops and loop clouds: A Numerical approach to the worldline formalism in QED[END_REF] in order to study situations, such as inhomogeneous background fields, for which no analytical solution is readily available.

Let us first present a simple algorithm specific to the case of a constant and uniform background field. Although this algorithm cannot be readily generalized to inhomogeneous background fields, its discussion is useful in order to illustrate the difficulties awaiting us in the general case. For an electrical field E in the z direction, we can choose a gauge in which the Euclidean gauge potential reads 

A i = (0, 0, 0, -iEx 3 ) , (194) 
(Since the loop is closed, there is no term in X 3 left in the exponential.) The quantity A is the projected area of the loop on the plane 34 in Euclidean space-time. Indeed, by Stokes' theorem, the abelian Wilson loop can be rewritten as the integral of F µν dσ µν over a surface whose boundary is the loop under consideration, where F µν is the strength of the background field and dσ µν is the measure on this surface. Note that A is an algebraic area, i.e. it is weighted by the winding number of the loop (and consequently it can also be negative).

Since the probability distribution of the worldlines is Gaussian in ṙ (see eq. ( 189)), it is easy to perform the path integral to obtain the probability distribution for the area A,

P τ (A) = π 4τ 1 cosh 2 πA 2τ . (196) 
Thus, for a given τ , the typical worldlines have an area that reaches values up to A ∼ τ .

It is therefore natural to rescale the area A by a factor τ by defining I ≡ τ A (this can be done by rescaling the coordinates about the barycenter, r i by r i ≡ √ τ y i ). Introducing also a rescaled imaginary fictitious time s ≡ -iτ /eE, we have

V E ,1 loop = eE 4π 2 d 4 X ∞ 0 ds s 3 e -i(m 2 /(eE))s e -isI . (197) 
After this rearrangement of the formula, all the dependence on the external field is contained in the quantity m 2 /(eE), while the factor e -isI is an (background field-independent) average over the rescaled worldlines that can be calculated once for all as a function of s.

A Fourier transform of this quantity then gives the imaginary part of the effective action for all values of m 2 /(eE), at the expense of a very small computational cost. The main difficulty of this algorithm is at small field strength (i.e. m 2 /(eE) 1), for which the effective action probes the s dependence of e -isI at small values of s. This region is dominated by large worldlines, with |I| m 2 /(eE). If the worldline average is obtained by an unbiased Monte-Carlo sampling over the ensemble of loops, very few large loops will be probed, thereby severely limiting the accuracy for weak fields eE m 2 . In the figure 13, we show a comparison of the exact 1-loop result and the result of such a numerical calculation. One can clearly see the increase of the Monte-Carlo statistical errors at small values of eE/m 2 .

In order to try to evade this limitation, let us return to the real valued fictitious time τ . Of course, using the probability distribution (196) and eq. ( 195) leads easily to the analytical result of eq. (192). However, in view of applications to more general situations, let us assume that we want to perform this integral numerically. Thus, we need to worry about its domain of convergence. The probability distribution of the area A decreases as exp(-π|A|/τ ) for large areas, and is weighted by the Wilson loop given in eq. (195). Convergence is therefore guaranteed only for |τ | < π/(eE). This result is of course not surprising: the domain of convergence cannot extend beyond the first pole τ 1 of the integral. For larger values of τ , the integral over A needs to be obtained by an analytic continuation, but this is possible only if the probability P τ (A) is known analytically. Unfortunately, in the case of an inhomogeneous background field, such an analytical knowledge is not available. It is still possible to parameterize the Wilson loop as Figure 13: Comparison of the exact 1-loop result and the result of a numerical worldline computation based on the Fourier transform method. From ref. [START_REF] Gies | Pair production in inhomogeneous fields[END_REF].

W X [r(τ )] = e -eE(X)τ I , (198) 
provided that we generalize the definition of I as follows

I ≡ i τ E(X) τ 0 dτ ṙi (τ ) A i (X + r(τ )) . (199) 
Note that A ≡ τ I is no longer simply related to the geometrical area of the worldline when the background field is not homogeneous. The probability distribution of I now depends on the base point X and the details of the background field around this point. In order to obtain an approximate analytical expression for the distribution of I, it is possible to make an ansatz that generalizes eq. ( 196),

P X (I) = N 1 cosh 2ν π αI 2 , ( 200 
)
where α and ν are X-dependent modifications of the width and shape of the distribution (the prefactor N is not independent, since it is fixed by normalization). These parameters can be fitted from a large enough ensemble of loops generated by Monte-Carlo. Once these parameters have been determined, the integral over I can be done analytically, +∞ -∞ dI P X (I) e -eE(X)τ I = N 4 ν πα

Γ(ν + eE(X)τ πα )Γ(ν -eE(X)τ πα ) Γ(2ν) . (201) 
This expression provides the answer for all values of the background field and all real fictitious times, thus solving the analytical continuation problem that one would face when doing the integral over I by Monte-Carlo. The poles that are responsible for the imaginary part of the effective action are given by the second gamma function in the numerator,

τ n = πα(n + ν) eE(X) , with residue Res W X , τ n = N 4 ν πα Γ(2ν + n) Γ(2ν) (-1) n l! d dτ (ν -eE(X)τ πα ) τn .
(202) This algorithm has been applied in ref. [START_REF] Gies | Pair production in inhomogeneous fields[END_REF] to the following x 1 dependent Sauter potential,

A 0 = -a tanh(kx 1 ) , E 1 = ak cosh 2 (kx 1 ) , (203) 
inverse propagator by the following discrete formula, 2 d

a 2 -D i D i + m 2 -1 = ∞ n=0 1 - a 2 -D i D i + m 2 2 d n , (204) 
where we have defined d ≡ d + 1 2 m 2 a a . Note that this is an exact formula. Likewise, the logarithm can be written as

ln a 2 2 d -D i D i + m 2 = - ∞ n=1 1 n 1 - a 2 -D i D i + m 2 2 d n . ( 205 
)
The trace that must be applied to the logarithm of the propagator is best calculated in position space, as a sum over all the lattice sites,

Tr O = sites x x O x . ( 206 
)
Let us define a probability distribution P 0 on the lattice, localized at y = x, i.e. P 0 (y, x) ≡ δ x,y . From P 0 , we can define a sequence of distributions P n defined iteratively by

P n+1 (y, x) = sites z 1 - a 2 -D i D i + m 2 2 d y,z P n (z, x) . (207) 
In terms of these P n , the trace of the logarithm of the inverse propagator reads Tr ln

a 2 2 d -D i D i + m 2 = - ∞ n=1 1 n sites x P n (x, x) . ( 208 
)
On the lattice, the covariant derivatives have a simple expression as finite differences weighted by link variables that encode the background gauge potential. For instance, for the direction x 1 :

[-D 2 1 f ] i = 2f i -U 1,i f i+1 -U -1 1,i-1 f i-1 a 2 , (209) 
where U 1,i is the link variable in the direction x 1 from the point i to the point i+1 (therefore, U -1 1,i-1 can be viewed as a Wilson line oriented in the opposite direction, from the point i to the point i -1). Using this definition, one obtains the following explicit formula for P n+1 in terms of P n :

P n+1 (i • • • , x) = 1 2 d U 1,i••• P n (i + 1 • • • , x) + U -1 1,i-1••• P n (i -1 • • • , x) + • • • , (210) 
where only the first coordinate has been written explicitly. Thus, P n+1 (y, x) is obtained from the values of P n (y , x) on the nearest neighbors y of the point y, with weights determined by the link variables that start at the point y (divided by 2 d). From this formula, one obtains a simple geometrical interpretation of P n (y, x):

P n (y, x) = 1 (2 d) n γ∈Γn(y,x) ∈γ U , (211) 
where Γ n (y, x) is the set of all the paths of length n on the lattice that start at the point x and end at the point y, and ∈γ U is the product of all the link variables encountered along the path γ. Therefore, one arrives at the following expression for the trace of the logarithm of the inverse lattice propagator:

Tr ln -D i D i + m 2 = const - ∞ n=0 1 n 1 (2 d) n sites x γ∈Γn(x,x) ∈γ U . ( 212 
)
The sum is over all the closed paths on the lattice, and the products of the link variables along such a path forms a Wilson loop 22 , which is the discrete analogue of the Wilson loop encountered in the continuous version of the worldline formalism exposed earlier in this section. A few of these closed paths are shown in the illustration of the figure 15. Note that the number of hops n must be even for closed paths on a cubic lattice, n = 2m. The total number of paths of length n on a d-dimensional lattice is (2d) n . The mass introduces an exponential suppression of the weight of long paths in the sum over the length n in eq. ( 212), since we can bound the summand by (d/ d) n = (1 + m 2 a 2 2d ) -n . The number of closed paths of length n = 2m is given by the following formula γ∈Γ2m(x,x)

1 = n1+•••+n d =m (2m)! n 1 ! 2 • • • n d ! 2 , ( 213 
)
and although it is considerably smaller than the total number of paths, it grows too quickly for an exhaustive enumeration23 of all the closed paths to reach the cutoff imposed by the mass. In refs. [START_REF] Schmidt | Determinant calculations with random walk worldline loops[END_REF][START_REF] Schmidt | Determinant calculations using random walk worldline loops[END_REF][START_REF] Schmidt | Matter determinants in background fields using random walk world line loops on the lattice[END_REF], it was proposed to depart from an exhaustive listing of the loops, and to use instead a statistical sampling that spans a much larger range of lengths. The ensemble of loops needs to be generated only once for a given set of lattice parameters, and can be reused for any background field on this lattice. Note however that, it may be advantageous to choose the ensemble of loops according to the physical length scale relevant for the problem under consideration (which in turn depends on the strength of the background field).

Worldline instanton approximation

Let us finally discuss a semi-classical approximation 24 of the worldline expression of the 1-loop effective action [START_REF] Dunne | Worldline instantons and pair production in inhomogeneous fields[END_REF][START_REF] Dunne | Worldline instantons. II. The Fluctuation prefactor[END_REF][START_REF] Dunne | Pair creation in inhomogeneous fields from worldline instantons[END_REF][START_REF] Dunne | Multidimensional Worldline Instantons[END_REF][START_REF] Dumlu | Complex Worldline Instantons and Quantum Interference in Vacuum Pair Production[END_REF][START_REF] Ilderton | Pair production from residues of complex worldline instantons[END_REF], which bears some resemblance with instantons in the fictitious time. For this reason, these semi-classical solutions are called worldline instantons. The starting point is the worldline representation of eq. ( 185) for the 1-loop effective action. By defining τ ≡ τ u and m 2 τ = s, we can firstly rewrite it as

V E ,1 loop = - ∞ 0 ds s e -s x i (0)=x i (1) Dx i (u) exp - m 2 4s 1 0 du ẋ2 4 +ie 1 0 du ẋi A i (x) , (214) 
where the path integration is now over closed loops of period 1. Since the rescaled fictitious time does not appear any longer as the integration bound inside the exponential, the integral over s can be performed in closed form, yielding a Bessel function:

V E ,1 loop = -2 x i (0)=x i (1) Dx i (u) K 0   m 1 0 du ẋ2   exp -ie 1 0 du ẋi A i (x) . (215)
In the regime where

m 2 1 0 du ẋ2 1 , (216) 
whose physical significance will be clarified later, it is possible to replace the Bessel function by the first term of its asymptotic expansion,

K 0 (z) ≈ z 1 π 2 e -z √ z . (217) 
This approximation leads to

V E ,1 loop ≈ - 2π m x i (0)=x i (1)
Dx i (u)

1 0 du ẋ2 -1/4 × exp -m 1 0 du ẋ2 + ie 1 0 du ẋi A i (x) . (218) 
The remaining path integral can be approximated by a stationary phase approximation, if we define the action

S ≡ m 1 0 du ẋ2 + ie 1 0 du ẋi A i (x) . (219) 
The stationarity condition is fulfilled for a closed path x i (u) that satisfies

m ẍi 1 0 du ẋ2 = ieF ij ẋj , (220) 
where F ij is the strength of the background field. Since F ij must be evaluated at the point x(u), this equation is a very non-trivial set of ordinary differential equations. By contracting eq. ( 220) with ẋi , the right hand side vanishes thanks to the antisymmetry of F ij and we therefore learn that ẋi

ẋi = const ≡ v 2 . ( 221 
)
The condition of eq. ( 216) is thus equivalent to mv 1. For each solution of the stationarity condition (220), the corresponding extremal value of the action gives a contribution to the imaginary part of the effective action

Im V E ,1 loop ∼ e -Sextremum . (222) 
Note that a more complicated calculation, involving the calculation of the determinant of the Gaussian fluctuations around these stationary solutions, is required in order to determine the prefactor in front of the exponential.

In order to illustrate this method, let us consider the case of a spatially dependent Sauter field

E(x 3 ) ≡ E cosh 2 (kx 3 ) , (223) 
that can be derived from the following gauge potential

A 4 = -i E k tanh(kx 3 ) . ( 224 
)
The equations of motion for the stationary solutions are

ẋ3 = v 1 -γ -2 tanh 2 (kx 3 ) , ẋ4 = -γ -1 v tanh(kx 3 ) , (225) 
where γ ≡ mk/(eE). There is a countable infinity of closed paths that satisfy these equations can be parameterized by

x 3 (u) = m eE 1 γ arcsinh γ 1 -γ 2 sin(2πn u) x 4 (u) = m eE 1 γ 1 -γ 2 arcsin (γ cos(2πn u)) , (226) 
where n is a strictly positive integer and u ∈ [0, 1]. The index n corresponds to how many times the closed path is traveled by the instanton solution.

In the limit of an extended field, γ → 0, these closed paths become circular (one could have checked directly that the closed paths that extremalize the action in the case of a constant and homogeneous electrical field are indeed circles). In contrast, when γ increases, these orbits become more and more elongated in the x 4 direction (see the figure 16) and they become infinitely large when γ reaches unity (because of the prefactor 1/ 1γ 2 in the second of eqs. (226)). In this limit, the value of the stationary action becomes infinite, and the rate of particle production goes to zero. This result was of course known from the exact solution to particle production problem in the Sauter field. Physically, it reflects the fact that a field must exist over a spatial domain as large as one Compton wavelength in order to be able to produce particles. This simple academic example is an illustration of a more general phenomenon, namely that spatial inhomogeneities tend to decrease particle production (this reduction can become total when the field becomes incoherent over scales comparable to the Compton wavelength).

7 Dynamically assisted Schwinger mechanism

Comparison of perturbative and non-perturbative particle production

As indicated by the form of the factor exp(-πm 2 /(eE)), the Schwinger mechanism in a constant field is completely nonperturbative, even in a weak field regime eE m 2 . This is because the field is time-independent, and therefore no process initiated by a finite number of photons can produce the energy of a pair. If the electrical field has a time dependence, a perturbative component can also contribute particle production [START_REF] Brezin | Pair production in vacuum by an alternating field[END_REF]. The transition between the perturbative regime and the non-perturbative regime can be characterized by the following dimensionless parameter,

γ ≡ m eEτ , (227) 
where τ is a characteristic time-scale of the time-dependence of the electrical field. This parameter has first been introduced by Keldysh in the context of the ionization of atoms [START_REF] Keldysh | Ionization in the field of a strong electromagnetic wave[END_REF], and is therefore called Keldysh parameter 25 . In this subsection, by using an analytically solvable time-dependent field, we explicitly show how the transition between the perturbative and the nonperturbative regimes happens [START_REF] Narozhnyi | The Simplist processes in the pair creating electric field[END_REF][START_REF] Taya | Finite pulse effects on e + epair creation from strong electric fields[END_REF]. Firstly, let us consider the distribution function obtained in a perturbative expansion up to the first order in a general background field. Since the distribution function is expressed by the mode function ψ - p,s (x), we need to find the perturbative solution of the equation for the mode function,

[iγ µ D µ -m] ψ - p,s (x) = 0 . (228) 
It can be easily obtained from the Green's formula

ψ(x) = y 0 =t0 d 3 y D 0 R (x, y) γ 0 ψ(y) -ie y 0 >t0 d 4 y D 0 R (x, y) / A(y) ψ(y) , (229) 
where D 0 R (x, y) is the free retarded fermion propagator:

D 0 R (x, y) = d 4 q (2π) 4 i(/ q + m) q 2 -m 2 + iq 0 e -iq•(x-y) .
(230)

By taking the limit of t 0 → -∞ and using the boundary condition [START_REF] Kluger | Fermion pair production in a strong electric field[END_REF], one can see that the first term in the right hand side of eq. ( 229) is simply the free spinor, propagated from t 0 to t. Therefore, the Green's formula also reads

ψ - p,s (x) = ψ free- p,s (x) -ie d 4 y D 0 R (x, y) / A(y) ψ - p,s (y) . (231) 
This equation has an iterative solution

ψ - p,s (x) = ∞ n=0 ψ (n)- p,s (x) , (232) 
where

ψ (n)- p,s (x) = (-ie) n d 4 y 1 • • • d 4 y n D 0 R (x, y 1 ) / A(y 1 ) • • • D 0 R (y n-1 , y n ) / A(y n ) ψ free- p,s (y n ) . (233 
) This formula describes the correction to the spinor due to n scatterings off the external field.

From the iterative solution (232), we can easily obtain the spectrum (108) at the first non-zero order in e. At this order, the gauge rotation factor U (x) has no effect (it alters the spectrum only in higher orders). Therefore the spin-averaged spectrum at the lowest-order reads If the Fourier-transformed gauge field Ã(q 0 + E p , p + p ) contains no pole in the q 0 -plane and is bounded when Im (q 0 ) → -∞, the q 0 -integration can be performed by the theorem of residues by picking the pole of the prefactor (q 0 -E p + i ) -1 (note that the exponential factor forces us to close the contour in lower half of the q 0 plane). After summing over the initial and final spin states, we have

dN (1) d 3 p = e 2 (2π) 3 2E p 1 2 s,s d 3 p (2π) 3 2E p × lim x 0 →+∞ ū(p, s) dq 0 2π e -iq 0 x 0 q 0 -E p + i / A(q 0 + E p , p + p )v(p , s ) 2 , (234 
dN (1) d 3 p = 2e 2 (2π) 3 2E p d 3 p (2π) 3 2E p p µ p ν + p ν p µ -(p•p + m 2 )g µν × Ãµ (E p + E p , p + p ) à * ν (E p + E p , p + p ) , (236) 
which is gauge invariant.

As an example, let us consider a sinusoidal time-dependent electrical field, E z (t) = E cos(ωt), for which the occupation number of the produced electrons is

(2π) 3 V dN (1) d 3 p = e 2 E 2 16E 2 p 1 - p z E p 2 2πδ(ω -2E p ) T , (237) 
where T ≡ 2πδ(0) stands for the time duration of the electrical field. The proportionality of the result to this time T is nothing but a manifestation of Fermi's golden rule. This lowestorder spectrum corresponds to the process γ → e + e -, where a single photon decays into an electron-positron pair. From the delta function, we see that this production mechanism is possible only if the frequency ω is above the threshold 2m, as expected from energy momentum conservation.

To compare the lowest-order perturbative result with the non-perturbative all-order result, we use the Sauter-type pulsed field of eq. ( 123), which can be obtained from the gauge potential

A 3 (t) = E τ tanh t τ . (238) 
The photons that make up this electrical field typically carry an energy ω ∼ τ -1 . In this gauge field background, the lowest-order perturbative spectrum reads

(2π) 3 V dN (1) d 3 p = e 2 E 2 1 - p z E p 2 π 2 τ 4 sinh 2 (πE p τ ) . (239) 
Because the gauge field (238) has a continuous spectrum in Fourier space, the perturbative production does not display any threshold behavior. Unlike with the sinusoidal electrical field, the particle production happens for any value of the pulse characteristic timescale τ . However, the particle production is most efficient if the photon energy is near the energy of the produced pair, 2E p . Indeed, as a function of τ , the spectrum has a peak around E p τ 0.61. The exact result for the spectrum in the same background, that includes the interaction with the background gauge field to all orders, is given by eq. ( 124). When changing the pulse duration τ , this spectrum has a clear transition between the perturbative regime and the non-perturbative regime [START_REF] Taya | Finite pulse effects on e + epair creation from strong electric fields[END_REF]. Since we have three dimensionful quantities, m, eE, and τ , the system is governed by two dimensionless parameters,

γ ≡ m eEτ , λ ≡ eEτ 2 . ( 240 
)
Taking the limit τ → +∞ (limit of constant electrical field) while keeping m and eE constant corresponds to γ 1 and λ 1. One can confirm that the exact spectrum [START_REF] Dumlu | Complex Worldline Instantons and Quantum Interference in Vacuum Pair Production[END_REF] goes to the spectrum due to a constant electrical field in this limit26 :

(2π) 3 V dN d 3 p -------→ γ 1, λ 1 exp - πm 2 ⊥ eE . (241) 
At the other extreme, let us consider the short pulse limit, τ → +∞, while keeping m and eE fixed. In terms of the dimensionless parameters γ and λ, this limit corresponds to γ 1 and λ 1. Note that taking the limit eE → 0 while keeping m and τ fixed would also provide the same behavior of γ and λ. Therefore, we expect that the lowest-order perturbative result (239) is recovered by this limit. Indeed, one can confirm that

(2π) 3 V dN d 3 p -------→ γ 1, λ 1 e 2 E 2 1 - p z E p 2 π 2 τ 4 sinh 2 (πE p τ ) . (242) 
In order to illustrate this behavior, we plot the exact result ( 124) and the lowest-order result (239) for several values of the Keldysh parameter in the figure 17. The p z -distributions at p ⊥ = 0 are shown in this plot. For large enough Keldysh parameter (γ 20), the two results show a good agreement. Conversely, as γ gets smaller, the discrepancy becomes larger. The transition between the perturbative regime and the non-perturbative regime is further illustrated in the figure 18. The maximum values of the momentum spectra (located at p = 0 for the perturbative spectrum, and at p = (0, 0, eEτ ) for the all-order one) are compared as a function of mτ , for several values of eE/m 2 . The constant field Schwinger result exp(-πm 2 /(eE)) is also indicated by thin black lines for reference. The top figure shows the subcritical case eE/m 2 < 1, and the bottom figure shows the supercritical case eE/m 2 > 1. For small τ (large γ), the all-order result and the lowest-order one agree well. For large τ (small γ), the all-order result converges to the constant field Schwinger result, while the perturbative computation breaks down. What is remarkable in the left plot (subcritical case) is that the spectrum can be several orders of magnitude larger than the constant field Schwinger result in the regime mτ ∼ 1. This means that, in subcritical electrical fields, eE/m 2 < 1, a time-dependence of the electrical field dramatically amplifies the particle production, especially if the typical energy carried by the time-dependent electrical field is near the perturbative threshold 2m.

Lattice numerical results

In the previous subsection, we have observed that the particle production in a time-dependent electrical field can be described by the lowest order perturbative calculation if the time-dependence is fast enough, such that γ 1 and λ 1. An important observation was that the particle production by the perturbative process γ → e + e -is much stronger than the non-perturbative Schwinger result if the field is subcritical eE < m 2 and the typical frequency of the time-dependent electrical field is near the threshold energy 2m.

However, if the particle production is a purely perturbative process, it is not much interesting since perturbative processes are well understood theoretically and have been verified experimentally in many circumstances. Our main interest in this review is the particle production in the non-perturbative regime. As we shall see, a similar enhancement mechanism is realized also in the non-perturbative regime if the electrical field is a superposition of a slowly varying and strong (but still subcritical) field and a fast and weak field [START_REF] Schutzhold | Dynamically assisted Schwinger mechanism[END_REF]. This phenomenon, called dynamically assisted Schwinger mechanism, has been the subject of many recent studies [START_REF] Piazza | Barrier control in tunneling e+ -e-photoproduction[END_REF][START_REF] Dunne | Catalysis of Schwinger Vacuum Pair Production[END_REF][START_REF] Monin | Semiclassical Calculation of Photon-Stimulated Schwinger Pair Creation[END_REF][START_REF] Orthaber | Momentum Spectra for Dynamically Assisted Schwinger Pair Production[END_REF][START_REF] Fey | Momentum dependence in the dynamically assisted Sauter-Schwinger effect[END_REF][START_REF] Jansen | Strongly enhanced pair production in combined high-and low-frequency laser fields[END_REF][START_REF] Augustin | Nonperturbative Bethe?Heitler pair creation in combined high-and low-frequency laser fields[END_REF][START_REF] Otto | Lifting shell structures in the dynamically assisted Schwinger effect in periodic fields[END_REF], since it opens up the possibility that the non-perturbative electron-positron pair production may be observed more easily than expected in experiments [START_REF] Piazza | Extremely high-intensity laser interactions with fundamental quantum systems[END_REF][START_REF] Otto | Dynamical Schwinger process in a bifrequent electric field of finite duration: survey on amplification[END_REF]. To be specific, we consider a superposition of two Sauter time-dependent electrical fields

E z (t) = E 1 cosh 2 t τ1 + E 2 cosh 2 t τ2 , (243) 
where E 1 E 2 and τ 1 τ 2 so that the first term represents a strong and slow field, while the second is a fast and weak field. In terms of the Keldysh parameter, these parameters are chosen so that γ 1 = m eE1τ1 1 and γ 2 = m eE2τ2 1. Therefore, the first pulse causes only non-perturbative particle production, while the second field can produce particles perturbatively. By superposing these two fields that have very different scales, an interplay between the non-perturbative and the perturbative physics occurs, leading to an important enhancement of the particle production.

Before presenting numerical evidence of this effect, let us qualitatively explain the mechanism at play. This is illustrated in the figure 19 (to be compared with the figure 1). In the standard Schwinger effect, a hole excitation must tunnel through the gap that separates the Dirac sea from the band of positive energy states, and there is therefore an exponential suppression in the length of this gap. As cartooned in the figure 19, a photon from the short weak field E 2 can slightly raise the energy of a hole excitation. But this effect alone is too weak to reach the positive energy band and produce an on-shell electron-positron pair. However, once it has been brought to this slightly higher energy, the hole excitation also has a shorter length to tunnel through in order to reach the positive energy band by quantum tunneling. Since the tunneling probability is exponentially sensitive to this length, even a moderate change in the length can produce an important increase in the yield. This is the essence of the dynamically assisted Schwinger mechanism. e , τ 1 = 10 -4 eV -1 and eE 2 = 0.025 m 2 e , τ 2 = 7 × 10 -6 eV -1 . p z -dependence at p ⊥ = 0 is shown. The curve labeled E 1 + E 2 shows the particle spectrum resulting from the sum of the two fields, while the curves labeled E 1 E 2 show the particle spectrum that results from the two fields considered separately. Following ref. [START_REF] Orthaber | Momentum Spectra for Dynamically Assisted Schwinger Pair Production[END_REF], we show the momentum spectra of particles produced in the electrical field (243). While the quantum kinetic approach was used in ref. [START_REF] Orthaber | Momentum Spectra for Dynamically Assisted Schwinger Pair Production[END_REF], we present the same results computed by real-time lattice calculations in the mode function method. In the figure 20, the momentum spectrum with the parameters eE 1 = 0.25 m 2 e , τ 1 = 10 -4 eV -1 = 510/m e and eE 2 = 0.025 m 2 e , τ 2 = 7 × 10 -6 eV -1 = 3.57/m e is plotted. The Keldysh parameters for each pulse are γ 1 0.078, and γ 2 11.2, respectively. Therefore, the first pulse is in the non-perturbative regime, while the second is in the perturbative regime. For comparison, the spectrum produced by the strong and slow pulse (curve labeled E 1 ) alone and that by the single weak and fast pulse (curve labeled E 2 ) alone are also shown. Although the single pulse E 2 is extremely weak and the spectrum it produces by itself is several orders of magnitude smaller than the spectrum produced by the pulse E 1 , superposing this weak and fast fields causes a dramatic enhancement of the production yield, by a factor around 4.

E 1 E 1 + E 2
In the lowest-order perturbative production studied in the previous subsection, the particle production is the most intense when the typical energy carried by photons constituting the electrical field is near to the threshold value 2m. In the case of the single Sautertype pulse, the particle yield by the lowest-order production takes maximum value when mτ 0.61. Also in the case of the dynamically assisted Schwinger mechanism in the double pulse (243), the maximum enhancement occurs when the typical energy carried by the weak pulse is near to 2m. In Fig. 21, the p z -spectra are plotted for several values of m e τ 2 . Other parameters are fixed to eE 1 = 0.25 m 2 e , m e τ 1 = 510, and eE 2 = 0.025 m 2 e . Within this parameter set, the enhancement is the most prominent at m e τ 2 0.61. Of course, the degree of the enhancement largely depends not only on τ 2 but also on other parameters and the profile of the electrical field as well. Optimization of the field profile to gain the maximum particle yield has been investigated in refs. [START_REF] Kohlfurst | Optimizing the pulse shape for Schwinger pair production[END_REF][START_REF] Li | Enhanced pair production in strong fields by multiple-slit interference effect with dynamically assisted Schwinger mechanism[END_REF][START_REF] Hebenstreit | Optimization of Schwinger pair production in colliding laser pulses[END_REF]. 

Summary and conclusions

In this review, we have covered some of the recent works related to the study of the Schwinger mechanism. Our main motivation in this review has been to present methods that have been designed in order to cope with more general backgrounds than those for which exact analytical solutions exist27 , allowing (numerical) studies of pair production in backgrounds with arbitrary space-time dependences (for smooth enough dependences, one may instead use a derivative expansion of the effective action [START_REF] Gusynin | Derivative expansion of the effective action for QED in (2+1)-dimensions and (3+1)-dimensions[END_REF][START_REF] Gusynin | Derivative expansion for the one loop effective Lagrangian in QED[END_REF]).

Two main classes of methods have been discussed in this review. The first one, that follows more conventional techniques of quantum field theory, is based on the possibility of expressing any one-loop amplitude in a background field in terms of mode functions. Loosely speaking, such one-loop amplitudes can be generated by small Gaussian fluctuations around a given background, and are therefore fully determined by the knowledge of a basis of small field perturbations about this background. This approach has a number of variants, that differ in the formulation but share common origins and are rigorously equivalent. The method of Bogoliubov transformations is directly connected to the mode functions by the fact that the Bogoliubov coefficients are obtained from the decomposition of the mode functions in Fourier modes. For spatially homogeneous electrical fields, it is possible to express the time evolution directly at the level of the Bogoliubov coefficients themselves, via quantum kinetic equations. These can be formulated either as two local in time equations for a pair of coupled distributions, one of which is the usual occupation number and the other is an anomalous distribution, or as a single equation, non-local in time, for the occupation number. When the background field is inhomogeneous, quantum kinetic equations can be generalized into an equation for the Wigner distribution of the system (i.e. the Wigner transform of the density operator).

The second family of methods that we have discussed in this review revolves around the worldline formalism, a formulation in which one-loop amplitudes in a background field are expressed as a path integration over all the closed paths in Euclidean space-time, parameterized by a extra fictitious time. Although it was originally derived from ideas in string theory, this representation has also some connections with Schwinger's proper time formulation for propagators. For a constant and homogeneous background, this approach leads easily to the well known analytical one-loop answer. In other cases, such as the Sauter background electrical field, one can obtain approximations valid in the weak field regime by finding extrema of a 5-dimensional action, called worldline instantons due to their resemblance with ordinary instanton solutions. For even more general backgrounds, numerical methods have been developed, in which one samples statistically the ensemble of worldlines. Although completely equivalent to more conventional methods, the worldline formalism emphasizes the space-time development of the production process, and brings some useful intuition over the relevant phenomena (e.g. the fact that in weak fields large worldlines play a very important role).

Besides these developments in the technical tools available for theoretical studies of the Schwinger mechanism, a number of works have focused on finding the type of external fields that may maximize the yield of produced particles, following the simple observation that the superposition of two fields results in a particle spectrum which is not simply the sum of the spectra that these two fields would produce individually (because the Schwinger mechanism is highly non-linear in the field). A simple yet spectacular such effect is the dynamically assisted Schwinger mechanism, obtained by superimposing a slow but intense field with a fast and much weaker field, both of which are much lower than the critical field. Intuitively, the short pulse raises the energy of a hole excitation, thereby shortening the tunneling length it has to overcome in order to be produced on-shell and considerably enhancing the yield. These phenomena could be used in order to greatly facilitate the experimental production of pairs by the Schwinger mechanism.

Let us finish this section by listing several topics related to the Schwinger mechanism that were not discussed in any detail in this review. Firstly, all the discussion of pair production by external fields in electrodynamics can be carried through to strong interactions in quantum chromodynamics. The most direct extension obviously concerns quark production in a chromo-electrical field [START_REF] Tanji | Quark pair creation in color electric fields and effects of magnetic fields[END_REF][START_REF] Cox | Q anti-Q pair creation: a field theory approach[END_REF][START_REF] Nayak | Non-perturbative quark-antiquark production from a constant chromo-electric field via the Schwinger mechanism[END_REF][START_REF] Iwazaki | Quark Pair Production in Expanding Glasma[END_REF][START_REF] Gyulassy | Quark and gluon pair production in SU(N) covariant constant fields[END_REF], but the Schwinger mechanism can also be relevant for gluon production [START_REF] Ambjorn | Canonical Quantization in Nonabelian Background Fields. 1[END_REF][START_REF] Yildiz | Vacuum Behavior in Quantum Chromodynamics[END_REF][START_REF] Nayak | Soft-gluon production due to a gluon loop in a constant chromo-electric background field[END_REF][START_REF] Tanji | Schwinger mechanism enhanced by the Nielsen-Olesen instability[END_REF][START_REF] Iwazaki | Anomalous Gluon Production and Condensation in Glasma[END_REF]. In heavy-ion collisions at high energies, the system is approximately boost-invariant in the longitudinal beam direction. The Schwinger mechanism in such a boost invariant expanding geometry has been studied in refs. [START_REF] Cooper | Particle production in the central rapidity region[END_REF][START_REF] Tanji | Pair creation in boost-invariantly expanding electric fields and two-particle correlations[END_REF]. In this review, we have only considered the case where the system is initially in the vacuum state. This can be generalized to situations where the system is initially non-empty, such as a thermal system [START_REF] Gavrilov | Density matrix of a quantum field in a particle-creating background[END_REF][START_REF] Gies | QED effective action at finite temperature: Two loop dominance[END_REF][START_REF] Hallin | Fermionic and bosonic pair creation in an external electric field at finite temperature using the functional Schrodinger representation[END_REF][START_REF] Kim | Schwinger pair production at finite temperature in scalar QED[END_REF][START_REF] Kim | Schwinger Pair Production at Finite Temperature in QED[END_REF][START_REF] Kleinert | Electron-positron pair productions in classical electric field and electromagnetic wave[END_REF], where quantum statistical effects will alter the production rate of particles by the external field. In QCD, it has also been argued that the Schwinger mechanism may lead to the dynamical generation of a gluon mass, that could explain some features of the Landau gauge propagators and verticein the soft sector [START_REF] Aguilar | Massless bound-state excitations and the Schwinger mechanism in QCD[END_REF].

Besides the production of elementary particles by external fields in vacuum, the Schwinger mechanism can also lead to the production of quasi-particle excitations in more exotic materials such as graphene [START_REF] Allor | The Schwinger mechanism and graphene[END_REF][START_REF] Zubkov | Schwinger pair creation in multilayer graphene[END_REF][START_REF] Gavrilov | Dirac fermions in strong electric field and quantum transport in graphene[END_REF][START_REF] Klimchitskaya | Creation of quasiparticles in graphene by a time-dependent electric field[END_REF][START_REF] Fillion-Gourdeau | Time-dependent pair creation and the Schwinger mechanism in graphene[END_REF], where it may be easier to achieve experimentally. The possibility of using ultracold atoms in an optical lattice as a simulator for the Schwinger mechanism was considered in ref. [START_REF] Kasper | Schwinger pair production with ultracold atoms[END_REF].

Among the theoretical tools used to study the Schwinger mechanism, let us mention a proposal to use stochastic quantization in ref. [START_REF] Fukushima | Schwinger Mechanism with Stochastic Quantization[END_REF], and some of the many works where holographic (e.g. based on the AdS/CFT correspondence) setups were considered in order to investigate particle production by external fields in the strong coupling regime [START_REF] Semenoff | Holographic Schwinger Effect[END_REF][START_REF] Ambjorn | Remarks on Holographic Wilson Loops and the Schwinger Effect[END_REF][START_REF] Bolognesi | Comments on Critical Electric and Magnetic Fields from Holography[END_REF][START_REF] Sato | Holographic description of the Schwinger effect in electric and magnetic fields[END_REF][START_REF] Sato | Potential Analysis in Holographic Schwinger Effect[END_REF][START_REF] Sonner | Holographic Schwinger Effect and the Geometry of Entanglement[END_REF][START_REF] Hashimoto | Vacuum Instability in Electric Fields via AdS/CFT: Euler-Heisenberg Lagrangian and Planckian Thermalization[END_REF][START_REF] Dietrich | Worldline holographic Schwinger effect[END_REF][START_REF] Hashimoto | Electromagnetic instability in holographic QCD[END_REF][START_REF] Kawai | A holographic description of the Schwinger effect in a confining gauge theory[END_REF]. Let us also finally mention another recent review on the Schwinger mechanism more focused on the backreaction and applications in astrophysics [START_REF] Ruffini | Electron-positron pairs in physics and astrophysics: from heavy nuclei to black holes[END_REF], as well as some works on the closely related question of pair creation in a curved space-time [START_REF] Kim | Schwinger Pair Production in dS(2) and AdS(2)[END_REF][START_REF] Garriga | Observer dependence of bubble nucleation and Schwinger pair production[END_REF][START_REF] Frob | Schwinger effect in de Sitter space[END_REF][START_REF] Kim | Geometric Origin of Pair Production by Electric Field in de Sitter Space[END_REF][START_REF] Cai | One-Loop Effective Action and Schwinger Effect in (Anti-) de Sitter Space[END_REF][START_REF] Fischler | Holographic Schwinger effect in de Sitter space[END_REF].

Figure 1 :

 1 Figure 1: Schematic picture of the tunneling process involved in the Schwinger mechanism. The white band is the gap between the anti-electron Dirac sea and the positive energy electron continuum, tilted by the potential V (x) = -Ex in the presence of an external electrical field E.

Figure 3 :

 3 Figure 3: Schematic diagrammatic representation of the two contributions to the derivative of the generating functions. Left: term in ϕ + ϕ -. Right: term in G +-.

Figure 5 :

 5 Figure 5: The spectrum (124) as a function of p z for various values of τ . The transverse mass is fixed to m ⊥ / √ eE = 0.1. A thin black horizontal line denotes exp(-πm 2 ⊥ /(eE)).

Figure 6 :

 6 Figure 6: The spectrum (124) as a function of p z for various values of m ⊥ . The pulse duration is fixed to √ eEτ = 4. Thin black horizontal lines denote exp(-πm 2 ⊥ /(eE)).

Figure 7 :

 7 Figure 7: Comparison of the spectrum computed by the statistical method and the analytic result (black lines) in the Sauter electrical field with √ eEτ = 2. This plot shows the p zspectra for fixed values of m ⊥ . The lattice parameters are N z = 128, √ eEa z = 0.16, and N conf = 256.

40 -Figure 8 :

 408 Figure8: Comparison of the two statistical methods; one using the single kind of stochastic field[START_REF] Zubkov | Schwinger pair creation in multilayer graphene[END_REF] and the other using the male-female fields[START_REF] Hashimoto | Electromagnetic instability in holographic QCD[END_REF].√ eE τ = 2 and m ⊥ / √ eE = 0.1. The lattice parameters are the same as those in the figure7.

Figure 9 :

 9 Figure 9: Comparison of the spectra computed by the statistical method and the analytic result (black line) in a weak field regime, eE = 0.25 m 2 , √ eE τ = 25.5. The lattice parameters are N z = 256, ma z = 0.048, and N conf = 48.

Figure 11 :

 11 Figure 11: The time dependence the energy density divided by E 0 = 1 2 E 2 0 . The energy density of the electrical field and that of produced particles are individually plotted as well as total energy density.

Figure 12 :

 12 Figure 12: Plots of the p z -spectrum with fixed transverse momentum p ⊥ = 0 with back reaction taken into account. The spectrum at several times is plotted (Top: early times, Bottom: later times). The parameters are the same as those in the figure 10.

  and the Wilson loop reads W X [r(τ )] = e -eEA , where we denote A ≡ τ 0 dτ ṙ4 (τ )r 3 (τ ) .

Figure 15 :

 15 Figure 15: Illustration of the lattice worldline formalism. The gray dot indicates the base point x.

Figure 16 :

 16 Figure16: Closed paths in the (x 3 , x 4 ) plane of the instanton solutions in the case of an x 3 -dependent Sauter electrical field, for various values of the parameter γ ≡ mk/(eE). γ = 0 corresponds to a circle of radius 1, and the closed paths become more and more elongated in the x 4 direction as γ increases.

  )where we have introduced õ (p)≡ d 4 x A µ (x) e ip•x . (235)This equation corresponds to computing the following 1-loop diagram:

Figure 17 :

 17 Figure 17: p z -dependence of the spectrum at p ⊥ = 0. The all-order exact result (124) (solid lines) and the lowest-order perturbative result (239) (dashed lines) are compared for several values of the Keldysh parameter γ. The ratio eE/m 2 = 1 is fixed. For γ = 20, the two lines are almost undistinguishable.

Figure 18 :

 18 Figure 18: Plots of the maximum value of the momentum spectrum as a function of mτ for several values of eE/m 2 . Top: subcritical field strength eE/m 2 < 1. Bottom: supercritical field strength eE/m 2 > 1. The all-order exact result (124) (solid lines) and the lowest-order perturbative result (239) (dashed lines) are compared. The constant field Schwinger result exp(-πm 2 /(eE)) is plotted as thin black lines.

Figure 19 :

 19 Figure 19: Schematic illustration of the dynamically assisted Schwinger mechanism.

Figure 20 :

 20 Figure 20: The momentum spectrum of electrons produced by the field (243) with the parameters eE 1 = 0.25 m 2e , τ 1 = 10 -4 eV -1 and eE 2 = 0.025 m 2 e , τ 2 = 7 × 10 -6 eV -1 . p z -dependence at p ⊥ = 0 is shown. The curve labeled E 1 + E 2 shows the particle spectrum resulting from the sum of the two fields, while the curves labeled E 1 E 2 show the particle spectrum that results from the two fields considered separately.

Figure 21 :

 21 Figure 21: The momentum spectrum of electrons produced by the field (243) for several values of m e τ 2 . The other parameters are fixed to eE 1 = 0.25 m 2 e , m e τ 1 = 510 and eE 2 = 0.025 m 2 e . The result with the single pulse E 1 is plotted as a black dashed line for comparison. The lattice parameters are the same as those in the figure 20.

For extra simplicity, we choose a potential whose absolute minimum is located at φ = 0, so that the perturbative expansion around the vacuum also corresponds to small field fluctuations. The mass is also important: when m = 0, producing a particle costs some energy, which is crucial to avoid infrared singularities.

This can be proven from Euler's formula for a graph, #nodes -#edges + #faces = 2 -2#holes, and #nodes = n V + n J + n E , #edges = n P and n L = #faces + 2#holes -1.

In the scalar theory under consideration in this section, the cutting rules realize the optical theorem at the level of single graphs. In non-abelian gauge theories, where ghosts cancel the unphysical gluon polarizations, similar cutting rules provide a realization of the optical theorem for groups of graphs that form a gauge invariant set.

Since the fields ϕ + and ϕ -are equal, there is no need to keep a subscript ± for these coefficients.

Its initial condition is real, and its equation of motion involves only real quantities.

Note that the n-particle spectrum defined in this way gives the expectation value of N (N -1) • • • (Nn + 1) when integrated over the momenta p 1 to pn :

It is precisely because in exclusive observables the final state is constrained that the boundary conditions for the fields cannot be purely retarded.

The simplest method is to write the equations of motion obeyed by G +-and its boundary conditions, and then to check a posteriori that eq. (50) satisfies both. Alternatively, it is possible to perform explicitly the summation implied by eq. (49), which avoids having to guess the structure of the answer.

This does not contradict what was said in the section 2.8. In this subsection, the difficulty was due to the fact that we were considering a strong source that couples directly to the fields we want to produce, and that these fields had self-interactions. Here, we are considering a much simpler problem since the fields we want to produce have no self-interactions.

Obviously, this result is true only for an homogeneous external field.

Our definition of the two particle spectrum dN ++ 2 /d 3 p 1 d 3 p 2 correspond to pairs of distinct particles. Therefore, its integral over p 1 and p 2 leads to the expectation value N + (N + -1) . This explains why the right hand side of the second equation contains a factor fp 1 instead of 1 + fp 1 .

The exact formula has a more complicated kz dependence, but this formula captures the main features of the result if the time x 0 is large, and in the strong field limit eE m 2 + k 2 ⊥ .

This is always possible with a gauge transformation.

Originally, Sauter[START_REF] Sauter | Zum kleinschen paradoxon[END_REF] studied a space-dependent electrical field E(x) = E/ cosh 2 (x/a) in the context of the Klein paradox. Both the space-dependent and the time-dependent Sauter electrical fields are amenable to an explicit analytic solutions. The comparison of these two situations is in fact very interesting to understand the essential differences between temporal inhomogeneities (that increase the particle yield) and spatial inhomogeneities (that reduce the particle yield). We will further comment on this difference in the subsection 6.5, when we discuss the worldline instanton approximation.

In case of scalar QED, this factor would be replaced by a Bose enhancement factor 1 + 2 fp(t ).

Note that only the negative energy mode functions appear in this vacuum expectation value. The positive energy mode functions are necessary to compute, for instance, the expectation value of the symmetrized charged current operator, J µ sym (x) ≡ 1 2 eγ µ [ ψ(x), ψ(x)]. However, if the system is charge neutral, the expectation value of the symmetrized charged current is the same as the expectation value of the unsymmetrized operator, J µ (x) ≡ e ψ(x)γ µ ψ(x).

This condition is satisfied for the energy-momentum tensor and the charge current, but the chiral charge is an exception.

If the background field is independent of one or more spatial coordinates, then the mode functions are plane waves in these coordinates, and the sum over the corresponding momenta can be done very efficiently by a Fast Fourier Transform, thereby trading a factor N 2 into N ln(N ) for each spatial direction for which such a simplification happens.

The formula written here is for scalar QED. In the case of an SU(Nc) gauge theory with scalars in the adjoint representations, the Wilson loop is an element of the adjoint representation of the SU(Nc) algebra, that must be traced in order to obtain the effective action.

In the simple case of a constant and uniform background electrical field in the direction 3, the Wilson loop depends only on the area of the loop projected in the plane 34. One can then use some known results on the moments of the distribution of these areas[START_REF] Epelbaum | From lattice Quantum Electrodynamics to the distribution of the algebraic areas enclosed by random walks on Z 2[END_REF][START_REF] Mingo | On the distribution of the area enclosed by a random walk onz2[END_REF] over the set Γn(x, x) in order to recover the standard result on the Schwinger mechanism in a constant electrical field.

Despite technical differences, this approach bears a lot of resemblance with the WKB semi-classical approximation, see e.g.[START_REF] Brezin | Pair production in vacuum by an alternating field[END_REF][START_REF] Popov | Pair production in a variable external field (quasiclassical approximation)[END_REF][START_REF] Kleinert | Electron-Positron Pair Production in Spaceor Time-Dependent Electric Fields[END_REF][START_REF] Strobel | Semiclassical pair production rate for time-dependent electrical fields with more than one component: WKB-approach and world-line instantons[END_REF][START_REF] Strobel | Semiclassical pair production rate for rotating electric fields[END_REF].

In refs.[START_REF] Brezin | Pair production in vacuum by an alternating field[END_REF][START_REF] Taya | Finite pulse effects on e + epair creation from strong electric fields[END_REF], the inverse of this parameter is defined as γ. We follow the original definition by Keldysh.

When taking the limit, one needs to keep the canonical momentum at the final time pz -eEτ constant, in order to obtain the result (241), which is independent of the longitudinal momentum.

Besides the well known analytical solutions with constant and Sauter-type fields, solutions have been obtained for more complicated setups in 1 + 1-dim massless QED, where the fermion fields can be elimnated by a bosonization procedure[START_REF] Chu | Capacitor Discharge and Vacuum Resistance in Massless QED 2[END_REF].

Acknowledgements

FG was supported by the Agence Nationale de la Recherche project 11-BS04-015-01.

whose known 1-loop analytical results serve as a reference for checking the accuracy of the numerical approach. The parameter k can be varied in order to test the algorithm at various spatial scales. The results of this test are shown in the figure 14 203) and the result of a numerical worldline computation based on the Ansatz of eq. ( 200). From ref. [START_REF] Gies | Pair production in inhomogeneous fields[END_REF]. imaginary part of the effective action normalized to the approximation of a locally constant field (i.e. where one uses the result for a constant field with the local value E(X)), as a function of the wavenumber k. The dashed curve is the analytical result by Nikishov. One sees that, while the locally constant field approximation quickly breaks down (the ratio starts deviating from one at very small values of k/m), the worldline numerical approach remains in good quantitative agreement with the exact result up to rather large values of k/m. Note that the limit k/m → 1, where the wavelength of the external field equals the Compton wavelength of the charged particles, is hard to cope with because it requires the e + e -pair to become extremely delocalized in order to reach an energy sufficient to become on-shell (when k = m, the particle yield is known to vanish with the Sauter background field). Consequently, this limit is dominated by large worldlines, whose sampling is difficult. Although imperfect, the ansatz of eq. ( 200) manages to capture the main behavior in this limit.

Lattice worldline formalism

The algorithms described in the previous section do not require any discretization of spacetime, which is arguably the most natural setup if the background field is known analytically.

In this subsection, we present an alternate formulation of the worldline formalism where the (Euclidean) space-time is first discretized on a cubic lattice [START_REF] Schmidt | Determinant calculations with random walk worldline loops[END_REF][START_REF] Schmidt | Determinant calculations using random walk worldline loops[END_REF][START_REF] Schmidt | Matter determinants in background fields using random walk world line loops on the lattice[END_REF][START_REF] Epelbaum | From lattice Quantum Electrodynamics to the distribution of the algebraic areas enclosed by random walks on Z 2[END_REF]. Although it has not yet been used to evaluate tunneling phenomena such as the Schwinger mechanism, such a lattice version of the worldline formalism seems promising when the background itself is obtained as the result of a lattice computation. Let us call a the lattice spacing (for simplicity, we assume that the lattice has identical lattice spacings in all directions, but it is easy to depart from this restriction).

On the lattice, it is more convenient to replace the representation of eq. ( 183) for the