Spectral representation of the particle production out of equilibrium—Schwinger mechanism in pulsed electric fields

Quantum field theory has been quite successful in describing non-trivial contents of the vacuum and the S-matrix elements from the vacuum-to-vacuum amplitude using the Lehmann–Symanzik–Zimmermann (LSZ) reduction formula. For static quantities the numerical Monte-Carlo simulation of the lattice discretized theory is so powerful that one could perform the first-principle calculation in a non-perturbative way. In contrast to tremendous achievements for static quantities, the numerical machinery for solving the real-time dynamics (or the initial value problem) has not been established beyond the scope of the linear response theory. The point is that one should compute not the vacuum-to-vacuum amplitude but an expectation value at time t, which involves time-evolution operators from $-\infty$ to t and also its Hermite conjugate. It is known as the closed-time path (CTP) formalism [1, 2] how to deal with two time-evolution operators with generalized Greenʼs functions. The microscopic derivation of the Boltzmann equation was pioneered by Kadanoff and Baym [3] (see also [4, 5] for recent studies), in which the Wigner transform of correlation functions translates to the distribution function and the spectral function.

Generally speaking, the spectral functions provide us with detailed information on physical contents in the system. Even in the case of equilibrated matter the spectral function represents in-medium dispersion relations of collective excitations such as the plasmon, the zero sound, etc (see [6, 7] for classical textbooks). One can even infer the real-time properties near equilibrium by the analytical continuation of Greenʼs functions once a spectral function is available. In this way, using the spectral function (or the imaginary part of the retarded self-energy), one can evaluate the thermal emission rate of a pair of particle and anti-particle (or hole in condensed matter systems) [8, 9]. Such thermal processes are allowed in a medium where thermally excited particles are brought in. In this kind of calculation in equilibrated matter, the translational invariance in time needs not be violated and the ordinary field-theory techniques are useful.

A more non-trivial example of the particle production is the process induced by the presence of time-dependent external field. The pioneering work by Heisenberg and Euler [10] has revealed that the one-loop effective action on top of electromagnetic background fields has an imaginary part. This indicates that the vacuum becomes unstable; in other words, the particle production can occur from the vacuum. The vacuum permittivity has also been formulated in the field-theoretical manner by Schwinger [11]. Named after his seminal work, the pair production of particle and anti-particle from the vacuum under electric field is commonly referred to as the Schwinger mechanism (see [12] for a comprehensive review). This could be regarded as a special example of the Landau–Zener effect (see [13] for example).

The essence for the particle production from the vacuum is concisely represented by the Bogoliubov transformation of the creation/annihilation operators, with which the positive- and the negative-energy states are mixed together. We note that the celebrated Hawking radiation, that is the particle production under gravitational effects, belongs to the same class of physics. In short, the vacuum defined in the 'in-state' is filled with particles and anti-particles if seen in the 'out-state' where the observer stands, and the Bogoliubov transformation connects the in- and the out-states by a unitary rotation. It should be mentioned that this unitary rotation would be a trivial phase if there is no external field, and so, no particle production occurs then.

The Schwinger mechanism and the Hawking radiation are quantum (tunnelling) phenomena and have been intensively studied in the semi-classical method like the Wentzel–Kramers–Brillouin (WKB) approximation (see [14] for a review, and also [15] for the WKB formulation of the Hawking radiation). Although the mixing via the Bogoliubov transformation is straightforward and the semi-classical methods appeal to our intuition, it would be more desirable to develop a systematic formulation in terms of the field correlators. We believe that this reformulation is indispensable for future progresses; someday one might be able to execute real-time numerical simulations, and then, the canonical quantization with creation/annihilation operators is not quite compatible with numerical algorithms. Ideally, if we can express the Schwinger mechanism using some spectral functions in analogy with the thermal emission rate, we could attain a unified view of the particle production near and out of equilibrium. We also comment that the Schwinger particle production is not much different from the thermal emission in the formalism; the latter arises from perturbative processes in media, while the former is also possible perturbatively as a result of interaction with time-dependent gauge potential from which an energy is supplied. The difference in physics lies in the fact that the non-perturbative resummation in the former case of the Schwinger problem would usually lead to an essential singularity with respect to the coupling constant associated with a tunneling phenomenon. So, in this sense, the non-perturbative resummation by solving the equation of motion always involves a mixed contribution from perturbative processes too.

Some time ago, the present author formulated the Schwinger mechanism in a form similar to the LSZ reduction formula in [16], which is based on preceding works [17, 18]. A variant of this formula is also used in a recent attempt to utilize the classical statistical approximation to simulate the Schwinger mechanism numerically [19]. As we discuss later, though the LSZ-type formula looks reasonable, the treatment of the in-state has some subtlety. If we consider the inclusive spectrum only, in fact, we can easily derive another formula given in terms of the spectral functions without any ambiguity. In this case the translational invariance in time is lost and we should handle the dynamical (time-dependent) spectral functions. It should be thus a natural idea to look into the temporal change of the spectral functions in accord with the quasi-particle contents affected by the time-dependent background fields. We note that our approach is along the similar lines to the Wigner formalism [20] from which the quantum kinetic theory is derived [21–23].

We stress that reformulating the problem of the particle production makes an important building block in a timely subject; real-time dynamics is the key issue in various fields of physics. In the research of the quark-gluon plasma formation for instance, the thermalization process is under active dispute (see [24] and references therein). Large laser facilities are aiming to detect the production of a pair of electron and positron and it has been discovered that the dynamically assisted Schwinger mechanism significantly reduces the critical strength of the electric field [25, 26]. For precise theoretical predictions it is strongly demanded to invent a new scheme for the full quantum real-time simulation. Probably, to achieve this goal, the stochastic quantization is one of the most promising approaches [27]. However, the conventional description of the Schwinger mechanism or the Hawking radiation does not fit in with the functional language with which the stochastic quantization is written. This highly motivates us to think of the spectral representation of the particle production out of equilibrium, as is the main topic of this work.

In this paper we will first give a detailed account of the derivation of our formula with the spectral functions. Then, we will investigate the general properties of the spectral functions associated with the in- and the out-states involving the Bogoliubov transformation. We can understand that the Schwinger mechanism accesses only a small portion of the whole spectral functions. This means, hence, that the spectral functions contain much more information than the Schwinger mechanism and new possibilities for a better detection might be still buried in them. The dynamical spectral functions thus deserve serious investigations and we will construct them concretely for a special case of homogeneous pulsed electric fields to dig non-trivial features out.

Let us consider a general setup to formulate the particle production out of equilibrium in quantum field theory. In this paper we focus only on a single-component complex scalar field (i.e. scalar QED [28]) to simplify the expressions, but the generalization to other fields such as fermions and multi-component fields is not difficult [29].

We require the existence of well-defined asymptotic states, namely, the in-state at $t=-\infty$ and the out-state at $t=\infty$, where the interactions should be turned off. In our convention we put 'in' and 'out' in the subscript to refer to quantities that belong to the in-state and the out-state, respectively. With increasing time, thus, the energy dispersion relation should evolve from ${{p}_{0}}={{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$ to ${{p}_{0}}={{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)$ driven by interactions with external fields, and we would like to compute the particle number associated with this change. For this purpose the expression for the inclusive spectrum is our starting point, which is given by the number operator as

$$\begin{eqnarray}&&\frac{{\rm d}N}{{{{\rm d}}^{3}}{\boldsymbol{p}} }=\frac{1}{{{\left( 2\pi \right)}^{3}}}\left\langle {{{\hat{\rho }}}_{{\rm in}}}\;{{{\hat{n}}}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right\rangle :=\frac{1}{{{\left( 2\pi \right)}^{3}}}\left\langle {{{\hat{\rho }}}_{{\rm in}}}\;\hat{a}_{{\rm out}}^{\dagger }\left( {\boldsymbol{p}} \right){{{\hat{a}}}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right\rangle .\end{eqnarray} \tag{ 1 }$$

Here ${{\hat{\rho }}_{{\rm in}}}$ represents the density matrix that characterizes the in-state. If we choose it to be a pure state of the initial vacuum, ${{\hat{\rho }}_{{\rm in}}}=|{{0}_{{\rm in}}}\left. {} \right\rangle \left\langle {} \right.{{0}_{{\rm in}}}|$, there is no contribution to (1) from the initial state. Then (1) counts the number of produced particles only. We make a remark that, if the initial state contains particles, we may utilize (1) to address the problem of particle absorption as well as particle production.

Our goal at the moment is to find an alternative expression of (1) in terms of field variables instead of creation/annihilation operators. To this end we need a prescription to identify creation/annihilation operators under background fields. These operators are related to the field operator $\hat{\phi }\left( x \right)$ via the expansion on complete basis, which is a clean procedure in the asymptotic states. In the out-state the annihilation operator is extracted through

$$\begin{eqnarray}\begin{array}{rcl} \sqrt{2{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\;{{{\hat{a}}}_{{\rm out}}}\left( {\boldsymbol{p}} \right) & = & \mathop{{\rm lim} }\limits_{t\to \infty }\;{\rm i}\int {{{\rm d}}^{3}}{\boldsymbol{x}} \;{{{\rm e}}^{{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)t-{\rm i}{\boldsymbol{p}} \cdot {\boldsymbol{x}} }}\;\left[ {{\partial }_{t}}-{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right]\hat{\phi }\left( t,x \right) \\ {} & = & \mathop{{\rm lim} }\limits_{t\to \infty }\;{\rm i}\;{{{\rm e}}^{{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)t}}\;\left[ {{\partial }_{t}}-{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right]\;\hat{\phi }\left( t,{\boldsymbol{p}} \right). \\ \end{array}\end{eqnarray} \tag{ 2 }$$

We use the same notation $\hat{\phi }$ also for the Fourier transformed field as long as no confusion arises. In our convention the normalization above is consistent with the commutation relation, $\left[ {{{\hat{a}}}_{{\rm out}}}\left( {\boldsymbol{p}} \right),\hat{a}_{{\rm out}}^{\dagger }\left( {\boldsymbol{p}} ^{\prime} \right) \right]={{\left( 2\pi \right)}^{3}}\delta \left( {\boldsymbol{p}} -{\boldsymbol{p}} ^{\prime} \right)$. This (2) is a basic relation frequently used in the derivation of the LSZ reduction formula in many textbooks. Because the number operator involves the creation/annihilation operators at the equal time, we can drop the exponential part and simplify the formula as

$$\begin{eqnarray}&&\hat{a}_{{\rm out}}^{\dagger }\left( {\boldsymbol{p}} \right){{\hat{a}}_{{\rm out}}}\left( {\boldsymbol{p}} \right)=\frac{1}{2{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\mathop{{\rm lim} }\limits_{{{t}_{1}}={{t}_{2}}=t\to \infty }\left[ {{\partial }_{{{t}_{1}}}}+{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right]\;\left[ {{\partial }_{{{t}_{2}}}}-{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right]\;{{\hat{\phi }}^{\dagger }}\left( {{t}_{1}},{\boldsymbol{p}} \right)\hat{\phi }\left( {{t}_{2}},{\boldsymbol{p}} \right).\end{eqnarray} \tag{ 3 }$$

If we are interested in the production of anti-particles, we can find a similar formula replacing $\hat{a}_{{\rm out}}^{\dagger }\left( {\boldsymbol{p}} \right){{\hat{a}}_{{\rm out}}}\left( {\boldsymbol{p}} \right)$ with $\hat{b}_{{\rm out}}^{\dagger }\left( -{\boldsymbol{p}} \right){{\hat{b}}_{{\rm out}}}\left( -{\boldsymbol{p}} \right)$. Owing to the conservation of U(1) charge (electric charge), the number of produced anti-particles should be anyway identical with that of produced particles, so we do not calculate it explicitly here.

In view of this form it is already clear that we can translate (1) into a representation by means of the Wightman function [30]. If we define a Wigner transform of the Wightman function, we can interpret it as a distribution function and we can identify the corresponding quantum kinetic equation [20–23]. The particle production is then to be characterized in terms of the solution of the kinetic equation. In this work we will pursue an alternative expression using the spectral functions instead of the Wightman function, for the spectral functions provide us with more intuition about physical contents of the system. Let us define the spectral functions or the Wigner transformed Pauli–Jordanʼs (denoted by $\mathcal{A}$) and Hadamardʼs (denoted by $\mathcal{D}$) functions as follows;

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{A}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right) & := & \frac{1}{V}\int _{-\Delta }^{\Delta }{\rm d}\delta t\;{{{\rm e}}^{{\rm i}{{p}_{0}}\delta t}}\;\left\langle {{{\hat{\rho }}}_{{\rm in}}}\;\left[ \hat{\phi }\left( t+\frac{1}{2}\delta t,{\boldsymbol{p}} \right),{{{\hat{\phi }}}^{\dagger }}\left( t-\frac{1}{2}\delta t,{\boldsymbol{p}} \right) \right] \right\rangle , \\ {{\mathcal{D}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right) & := & \frac{1}{V}\int _{-\Delta }^{\Delta }{\rm d}\delta t\;{{{\rm e}}^{{\rm i}{{p}_{0}}\delta t}}\;\left\langle {{{\hat{\rho }}}_{{\rm in}}}\;\left\{ \hat{\phi }\left( t+\frac{1}{2}\delta t,{\boldsymbol{p}} \right),{{{\hat{\phi }}}^{\dagger }}\left( t-\frac{1}{2}\delta t,{\boldsymbol{p}} \right) \right\} \right\rangle . \\ \end{array}\end{eqnarray} \tag{ 4 }$$

Conventionally only ${{\mathcal{A}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right)$ is called the spectral function, while ${{\mathcal{D}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right)$ is called the symmetry or statistical two-point function. In this work, since we treat them on equal footing, we generalize the meaning of the spectral function to refer to both ${{\mathcal{A}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right)$ and ${{\mathcal{D}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right)$ collectively. Here, in (4), we put a volume factor V because we look at the same momenta ${\boldsymbol{p}}$ and trivially there arises $2\pi \delta \left( 0 \right)=V$. One could define the spectral functions with two momentum arguments, which would be useful in the presence of spatially modulated background fields as used in [20]. In this work, however, we consider only the spatially homogeneous case, so that the above definition (4) suffices for our goal. It should be mentioned that our definition of (4) explicitly depends on an extra parameter $\Delta$. In the Wigner transformation, usually, one formally takes $\Delta \to \infty$. For a practical application to the numerical analysis, a finite $\Delta$ as above would be legitimate, as we will discuss later. To extract information on the in- or the out-state, as a matter of fact, one should keep the ordering, $|t|\gg \Delta$, when we formally take $\Delta \to \infty$; otherwise the spectral functions are affected by the interaction even for t that is far outside of the interacting region. Roughly speaking, $\Delta$ should be interpreted as an 'observation time' with which the quasi-particle oscillation is resolved.

We can change the variables from ${{t}_{1}}$ and ${{t}_{2}}$ to $t=\frac{1}{2}\left( {{t}_{1}}+{{t}_{2}} \right)$ and $\delta t={{t}_{1}}-{{t}_{2}}$, so that we can finally arrive at the following formula,

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}N}{{{{\rm d}}^{3}}{\boldsymbol{p}} } & = & \mathop{{\rm lim} }\limits_{t\to \infty }\frac{V}{{{\left( 2\pi \right)}^{3}}}\int \frac{{\rm d}{{p}_{0}}}{2\pi }\;\frac{1}{4{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\left[ \frac{1}{4}\partial _{t}^{2}+{{\left( {{p}_{0}}+{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right)}^{2}} \right] \\ {} & {} & \quad \times \left[ {{\mathcal{D}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right)-{{\mathcal{A}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right) \right]. \\ \end{array}\end{eqnarray} \tag{ 5 }$$

It is important to stress that (5) does not rely on a choice of the integration range $\Delta$ in the definition of (4). This is because the ${{p}_{0}}$-integration picks $\delta \left( \delta t \right)$ up to realize ${{t}_{1}}={{t}_{2}}$ after taking each derivative on ${{t}_{1}}$ and ${{t}_{2}}$. Although the results should be the same regardless of $\Delta$, the physical picture becomes more vivid if we choose an appropriate value of $\Delta$. Using that ${{\mathcal{D}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right)$ and ${{\mathcal{A}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right)$ are even and odd functions of ${{p}_{0}}$, respectively, we can readily confirm that the contribution from anti-particles amounts to just the same answer, which should be guaranteed by the charge conservation.

Only for completeness let us make a remark on another expression with use of Greenʼs functions. We can perform the Wigner transform for the retarded and advanced propagators to define $D_{{\rm R}}^{\Delta }\left( t,{{p}_{0}},{\boldsymbol{p}} \right)$ and $D_{{\rm A}}^{\Delta }\left( t,{{p}_{0}},{\boldsymbol{p}} \right)$ as well as the Feynman (time-ordered) propagator, $D_{{\rm F}}^{\Delta }\left( t,{{p}_{0}},{\boldsymbol{p}} \right)$. Then, (5) is just equivalent with

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}N}{{{{\rm d}}^{3}}{\boldsymbol{p}} } & = & \mathop{{\rm lim} }\limits_{t\to \infty }\frac{V}{{{\left( 2\pi \right)}^{3}}}\int \frac{{\rm d}{{p}_{0}}}{2\pi }\;\frac{1}{2{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\left[ \frac{1}{4}\partial _{t}^{2}+{{\left( {{p}_{0}}+{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right)}^{2}} \right] \\ {} & {} & \quad \times \left[ D_{{\rm F}}^{\Delta }\left( t,{{p}_{0}},{\boldsymbol{p}} \right)-D_{{\rm R}}^{\Delta }\left( t,{{p}_{0}},{\boldsymbol{p}} \right) \right]. \\ \end{array}\end{eqnarray} \tag{ 6 }$$

This expression might be more tractable if one wants to apply the Schwinger–Keldysh formalism to compute the correlation functions.

Here we note that our formulas (5) and (6) look significantly different from the LSZ-type expression as used in [16, 19]. Because of the subtlety in the treatment of the time-dependent dispersion relation, it is quite non-trivial how to rewrite them into a covariant form that would lead to the LSZ-type formula and it would be theoretically cleaner to utilize our formula directly derived from the expectation value of the number operator. We give detailed discussions on this point in appendix A.

To gain some feeling about how our formula works, we will take a quick look at the typical behaviour of these spectral functions in the asymptotic states where we can expand the field in terms of plane waves.

Because the dynamical spectral functions are less known objects than more conventional ones in equilibrated matter, we will devote this section to the exploration of how they look like in the asymptotic in- and out-states. To reduce unnecessary complication, we shall limit our discussion to the choice of ${{\hat{\rho }}_{{\rm in}}}=\left| {{0}_{{\rm in}}} \right\rangle \left\langle {{0}_{{\rm in}}} \right\rangle$ for the moment. We will address an extension to the finite temperature/density environment in the later section. We denote the annihilation operators, ${{\hat{a}}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$ and ${{\hat{b}}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$, respectively, for particles and anti-particles, with which the vacuum $\left| {{0}_{{\rm in}}} \right\rangle$ is defined. For a practical purpose we take the range of t from $-T$ to T with a sufficiently large T.

The field operator in the in-state around $t=-T$ is then a superposition of the plane waves with ${{\hat{a}}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$ for the positive-energy oscillation and $\hat{b}_{{\rm in}}^{\dagger }\left( {\boldsymbol{p}} \right)$ for the negative-energy oscillation, i.e.

$$\begin{eqnarray}&&\hat{\phi }\left( t\sim -T,{\boldsymbol{p}} \right)=\frac{1}{\sqrt{2{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)}}\left[ {{{\hat{a}}}_{{\rm in}}}\left( {\boldsymbol{p}} \right){{{\rm e}}^{-{\rm i}{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)t}}+\hat{b}_{{\rm in}}^{\dagger }\left( -{\boldsymbol{p}} \right){{{\rm e}}^{{\rm i}{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)t}} \right].\end{eqnarray} \tag{ 7 }$$

The plane waves should be the solutions of the equation of motion around $t=-T$, and they evolve to a mixture of the positive- and the negative-energy states as t elapses toward the interacting region. We can then parametrize this mixing effect as

$$\begin{eqnarray}&&\frac{{{{\rm e}}^{-{\rm i}{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)t}}}{\sqrt{2{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)}}\ \ \left( t\sim -T \right)\ \longrightarrow \ \frac{{{\alpha }_{{\boldsymbol{p}} }}\;{{{\rm e}}^{-{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)t}}+\beta _{{\boldsymbol{p}} }^{*}\;{{{\rm e}}^{{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)t}}}{\sqrt{2{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}}\ \ \left( t\sim T \right)\ ,\end{eqnarray} \tag{ 8 }$$

where the Bogoliubov coefficients, ${{\alpha }_{{\boldsymbol{p}} }}$ and ${{\beta }_{{\boldsymbol{p}} }}$, are determined according to the equation of motion, and a similar relation should hold for another branch of solution starting with $\propto {{{\rm e}}^{{\rm i}{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)t}}$. In fact, if the Hamiltonian is Hermite, the complex conjugate of the above relation is true, so the field operator in the out-state at t = T then reads,

$$\begin{eqnarray}\begin{array}{rcl} \hat{\phi }\left( t\sim T,{\boldsymbol{p}} \right) & = & \frac{1}{\sqrt{2{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}}\left\{ \left[ {{\alpha }_{{\boldsymbol{p}} }}{{{\hat{a}}}_{{\rm in}}}\left( {\boldsymbol{p}} \right)+{{\beta }_{{\boldsymbol{p}} }}\hat{b}_{{\rm in}}^{\dagger }\left( -{\boldsymbol{p}} \right) \right]{{{\rm e}}^{-{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)t}} \right. \\ {} & {} & \left. \qquad +\left[ \alpha _{{\boldsymbol{p}} }^{*}\hat{b}_{{\rm in}}^{\dagger }\left( -{\boldsymbol{p}} \right)+\beta _{{\boldsymbol{p}} }^{*}{{{\hat{a}}}_{{\rm in}}}\left( {\boldsymbol{p}} \right) \right]{{{\rm e}}^{{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)t}} \right\}. \\ \end{array}\end{eqnarray} \tag{ 9 }$$

This defines the creation/annihilation operators in the out-state, and by requiring the canonical commutation relation for them, we can find the normalization condition, $|{{\alpha }_{{\boldsymbol{p}} }}{{|}^{2}}-|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}}=1$.

At this point we can immediately recover the known result for the Schwinger mechanism directly from (3). Applying the operator $\left[ {{\partial }_{t}}-{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right]$ on $\hat{\phi }$ we project the positive-energy part out, and then we can plug the number operator of (3) into (1), which yields an estimate of produced particles as

$$\begin{eqnarray}&&\frac{{\rm d}N}{{{{\rm d}}^{3}}{\boldsymbol{p}} }=\frac{1}{{{\left( 2\pi \right)}^{3}}}\left\langle {{0}_{{\rm in}}}|\beta _{{\boldsymbol{p}} }^{*}{{{\hat{b}}}_{{\rm in}}}\left( -{\boldsymbol{p}} \right)\;{{\beta }_{{\boldsymbol{p}} }}\hat{b}_{{\rm in}}^{\dagger }\left( -{\boldsymbol{p}} \right)|{{0}_{{\rm in}}} \right\rangle =\frac{V|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}}}{{{\left( 2\pi \right)}^{3}}}.\end{eqnarray} \tag{ 10 }$$

This is a standard formula for the particle production obtained via the Bogoliubov transformation [31]. Now it is intriguing to check how our formula (5) gives rise to the same answer.

We can immediately compute the spectral functions from the asymptotic forms (7) and (9) if we take $T\gg \Delta$. Then, the spectral functions at $t=-T$ have no access to the region with non-vanishing background fields, so they take a familiar expression just for non-interacting particles;

$$\begin{eqnarray}&&{{\mathcal{A}}^{\Delta }}\left( t\sim -T,{{p}_{0}},{\boldsymbol{p}} \right)=\frac{\pi }{{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)}\left[ \delta \left( {{p}_{0}}-{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right) \right)-\delta \left( {{p}_{0}}+{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right) \right) \right],\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\mathcal{D}}^{\Delta }}\left( t\sim -T,{{p}_{0}},{\boldsymbol{p}} \right)=\frac{\pi }{{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)}\left[ \delta \left( {{p}_{0}}-{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right) \right)+\delta \left( {{p}_{0}}+{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right) \right) \right].\end{eqnarray} \tag{ 12 }$$

In this case ${{\mathcal{D}}^{\Delta }}-{{\mathcal{A}}^{\Delta }}$ has only a term that is proportional to $\delta \left( {{p}_{0}}+{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right) \right)$, and thus the produced particle is vanishing as is obvious from (5). Now let us go into later time when these functions should change their shape. Once (9) eventually follows, it is just a simple arithmetic procedure to reach,

$$\begin{eqnarray}&&\quad {{\mathcal{A}}^{\Delta }}\left( t\sim T,{{p}_{0}},{\boldsymbol{p}} \right)=\frac{\pi }{{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\left[ \delta \left( {{p}_{0}}-{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right)-\delta \left( {{p}_{0}}+{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right) \right],\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}\begin{array}{rcl} \quad {{\mathcal{D}}^{\Delta }}\left( t\sim T,{{p}_{0}},{\boldsymbol{p}} \right) & = & \left[ |{{\alpha }_{{\boldsymbol{p}} }}{{|}^{2}}+|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}} \right]\;\frac{\pi }{{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\;\left[ \delta \left( {{p}_{0}}-{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right)+\delta \left( {{p}_{0}}+{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right) \right] \\ {} & {} & \qquad +\frac{2}{{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\operatorname{Re}\left[ {{\alpha }_{{\boldsymbol{p}} }}{{\beta }_{{\boldsymbol{p}} }}\;{{{\rm e}}^{-2{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)t}} \right]\;2\pi \delta \left( {{p}_{0}} \right). \\ \end{array}\end{eqnarray} \tag{ 14 }$$

Here, again, we required $T\gg \Delta$. There are two interesting observations as perceived from the above: (1) Pauli–Jordanʼs function ${{\mathcal{A}}^{\Delta }}$ is insensitive to the Bogoliubov transformation and the overall factor is $|{{\alpha }_{{\boldsymbol{p}} }}{{|}^{2}}-|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}}=1$. (2) Hadamardʼs function ${{\mathcal{D}}^{\Delta }}$ is affected by the mixing effect by $|{{\alpha }_{{\boldsymbol{p}} }}{{|}^{2}}+|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}}\ne 1$ and, besides, an interference term $\propto {{\alpha }_{{\boldsymbol{p}} }}{{\beta }_{{\boldsymbol{p}} }}$ appears. We emphasize that such an interference term is usually absent and is quite peculiar to the Bogoliubov mixing effect.

Then, the difference between these two spectral functions consists of three terms as follows,

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{D}}^{\Delta }}-{{\mathcal{A}}^{\Delta }} & = & \left| {} \right.{{\beta }_{{\boldsymbol{p}} }}{{\left| {} \right.}^{2}}\frac{2\pi }{{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\delta \left( {{p}_{0}}-{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right)+\left| {} \right.{{\alpha }_{{\boldsymbol{p}} }}{{\left| {} \right.}^{2}}\frac{2\pi }{{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\delta \left( {{p}_{0}}+{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right) \\ {} & {} & \qquad +\frac{2}{{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\operatorname{Re}\left[ {{\alpha }_{{\boldsymbol{p}} }}{{\beta }_{{\boldsymbol{p}} }}{{{\rm e}}^{-2{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)t}} \right]2\pi \delta \left( {{p}_{0}} \right). \\ \end{array}\end{eqnarray} \tag{ 15 }$$

We can make it sure by the explicit calculation that the second ($\propto \;|{{\alpha }_{{\boldsymbol{p}} }}{{|}^{2}}$) and the third ($\propto \;{{\alpha }_{{\boldsymbol{p}} }}{{\beta }_{{\boldsymbol{p}} }}$) terms have no finite contribution if applied to (5), and only the first ($\propto \;|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}}$) term is relevant to the particle production, which yields exactly the same answer as (10).

Although the calculations are very easy, the expression of (15) in the out-state has profound implications. In many situations we typically have ${{\alpha }_{{\boldsymbol{p}} }}\approx 1$ and $|{{\beta }_{{\boldsymbol{p}} }}|\ll 1$, for which the first term is much smaller than the third interference term. In the next section, indeed, we will numerically compute the spectral functions and confirm that this is the case. It would be an interesting future problem to think of a way to make use of the interference term in order to probe the Bogoliubov mixing effect experimentally.

Before closing this section, it would be instructive to understand how the standard propagators are modified by the Bogoliubov transformation. Surprisingly, we find that the retarded propagator is intact under the mixing effect and only the Feynman propagator depends on the Bogoliubov coefficients. That is,

$$\begin{eqnarray}&&D_{{\rm R}}^{\Delta }\left( t\sim T,{{p}_{0}},{\boldsymbol{p}} \right)=P\frac{{\rm i}}{p_{0}^{2}-E_{{\rm out}}^{2}\left( {\boldsymbol{p}} \right)}+\frac{\pi }{2{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\left[ \delta \left( {{p}_{0}}-{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right)-\delta \left( {{p}_{0}}+{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right) \right],\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}\begin{array}{rcl} D_{{\rm F}}^{\Delta }\left( t\sim T,{{p}_{0}},{\boldsymbol{p}} \right) & = & P\frac{{\rm i}}{p_{0}^{2}-E_{{\rm out}}^{2}\left( {\boldsymbol{p}} \right)}+\left( |{{\alpha }_{{\boldsymbol{p}} }}{{|}^{2}}+|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}} \right)\frac{\pi }{2{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\left[ \delta \left( {{p}_{0}}-{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right) \right. \\ {} & {} & \left. \quad +\delta \left( {{p}_{0}}+{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right) \right) \right]+\frac{1}{{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\operatorname{Re}\left( {{\alpha }_{{\boldsymbol{p}} }}{{\beta }_{{\boldsymbol{p}} }}{{{\rm e}}^{-2{\rm i}{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)t}} \right)2\pi \delta \left( {{p}_{0}} \right)\ , \\ \end{array}\end{eqnarray} \tag{ 17 }$$

where P stands for taking the principal value. It is clear at glance that ${{\mathcal{D}}^{\Delta }}-{{\mathcal{A}}^{\Delta }}=2\left( D_{{\rm F}}^{\Delta }-D_{{\rm R}}^{\Delta } \right)$ holds as it should.

Now it is time to take one step forward to understand how the spectral functions should evolve continuously from the in-state to the out-state as a function of t. In the aim of visualizing the behaviour with increasing t, we need to perform numerical calculations. In the next section we present our numerical results.

We here solve the equation of motion for given electric fields. In principle, we can numerically deal with arbitrary electric fields within our present approximation to neglect the back-reaction. Although the analytical solution is not demanded here, we shall choose one of the most well-investigated profile known as the Sauter potential [32], which is solvable and identifiable with a pulsed electric field,

$$\begin{eqnarray}&&E\left( t \right)=E\;{\rm sec}{{{\rm h}}^{2}}\left( \omega t \right).\end{eqnarray} \tag{ 18 }$$

The frequency parameter ω characterizes the life time of the applied electric field. Let us choose the z axis along the direction of the electric field, and then the associated vector potential reads,

$$\begin{eqnarray}&&{{A}_{z}}\left( x \right)=\frac{E}{\omega }\left[ {\rm tanh} \left( \omega t \right)-1 \right].\end{eqnarray} \tag{ 19 }$$

Then, we can find two independent solutions, $\psi _{p}^{\left( \pm \right)}\left( t \right)$, by solving the equation of motion under this vector potential,

$$\begin{eqnarray}&&\left[ \partial _{t}^{2}+{{\left( {{p}_{z}}+\frac{eE}{\omega }\left[ {\rm tanh} \left( \omega t \right)-1 \right] \right)}^{2}}+m_{\bot }^{2} \right]\psi _{p}^{\left( \pm \right)}\left( t \right)=0\ ,\end{eqnarray} \tag{ 20 }$$

where ${{m}_{\bot }}$ represents the transverse mass, $m_{\bot }^{2}:=p_{x}^{2}+p_{y}^{2}+{{m}^{2}}$. We should impose the following boundary conditions;

$$\begin{eqnarray}&&\psi _{p}^{\left( \pm \right)}\left( t\sim -T \right)=\frac{1}{\sqrt{2{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)}}\;{{{\rm e}}^{\mp {\rm i}{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)t}}\ ,\end{eqnarray} \tag{ 21 }$$

for large enough T, with the dispersion relations,

$$\begin{eqnarray}&&{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)=\sqrt{{{\left( {{p}_{z}}-2eE/\omega \right)}^{2}}+m_{\bot }^{2}}\ ,\qquad {{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)=\sqrt{p_{z}^{2}+m_{\bot }^{2}}.\end{eqnarray} \tag{ 22 }$$

We note that the condition (21) specifies not only the initial value but the initial derivative too, which are enough to characterize the unique solution of the second-order differential equation.

One can write the analytical expressions of $\psi _{p}^{\left( \pm \right)}\left( t \right)$ using the hyper-geometric functions. Therefore, the number of produced particle or $|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}}$ is analytically known. Hereafter we shall refer to all quantities with mass dimensions in unit of the electric field $\sqrt{eE}$ and present our results with dimensionless numbers. In this work we work with a specific choice of

$$\begin{eqnarray}&&\frac{\omega }{\sqrt{eE}}=1,\end{eqnarray} \tag{ 23 }$$

to investigate the effect of pulsed electric fields. A different choice of ω makes no qualitative change in our resulting spectral functions. We make a plot in figure 1 to show the analytical structure of $|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}}$ as a function of the longitudinal momentum ${{p}_{z}}$ and the transverse mass ${{m}_{\bot }}$.

Zoom In Zoom Out Reset image size

Obviously from figure 1, the particle number becomes greater for smaller ${{m}_{\bot }}$. The produced particles are accelerated to the positive z direction by the electric field, and as understood from (20), ${{p}_{z}}$ is shifted by $0\sim 2eE/\omega$ during the time evolution. This means that the momentum distribution of the produced particles should spread over ${{p}_{z}}=0\sim 2eE/\omega$. We can confirm this expectation explicitly on figure 1.

To discuss the effect of the particle production in a reasonably visible manner, we will look at the point of following momenta,

$$\begin{eqnarray}&&\frac{{{m}_{\bot }}}{\sqrt{eE}}=0\ ,\qquad \frac{{{p}_{z}}}{\sqrt{eE}}=\frac{\sqrt{eE}}{\omega }=1.5\ ,\end{eqnarray} \tag{ 24 }$$

which deviates from the pronounced peak seen in figure 1. From a theoretical point of view one might have thought that the exact peak position would be a better choice, but if we choose ${{m}_{\bot }}=0$ and ${{p}_{z}}=2$, for instance, the numerical calculations result in singularity out of control. From a pragmatic point of view, on the other hand, this choice might have looked too idealized; ${{m}_{\bot }}=0$ means ${{p}_{x}}={{p}_{y}}=0$ and m = 0. In reality, however, $\sqrt{eE}$ created by experiment is much smaller than the lightest particle (i.e. electron) mass and it is impossible for any present experiment to probe the regime of (24). As written shortly above (24), we made such a parameter choice intentionally for the illustration purpose here. We could have taken ${{m}_{\bot }}/\sqrt{eE}\gg 1$, but the particle production would be exponentially suppressed then, and so, it should be impossible that we perceive a secondary peak in the spectral functions corresponding to the Bogoliubov mixing. We should therefore keep in mind that the presented results with (23) and (24) are exaggerated for the ease of perception; however, qualitative physics is correctly captured.

With these parameters we solve the equation of motion (20) numerically to find $\psi _{p}^{\left( \pm \right)}\left( t \right)$, the results of which are shown in figure 2. We used the fourth-order Runge–Kutta (RK4) method and took 20 000 points to discretize along the time direction. We imposed the boundary conditions (21) at $t=-T$ with T = 100. Because the equation of motion is real, $\psi _{p}^{\left( - \right)}\left( t \right)=\psi _{p}^{\left( + \right)*}\left( t \right)$ follows immediately, and this means that the real part of them should be identical. This is why we present the imaginary part in figure 2, and indeed, we can make it sure that our numerical calculations go correctly to respect ${\rm Im}\psi _{p}^{\left( - \right)}\left( t \right)=-{\rm Im}\psi _{p}^{\left( + \right)}\left( t \right)$.

Zoom In Zoom Out Reset image size

The in-state around $t\sim -T$ has the field amplitude of the correct normalization $1/\sqrt{2{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)}$, while at later time, as seen in figure 2, the amplitude deviates from $1/\sqrt{2{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}$. This discrepancy is attributed to the mixing between the positive- and the negative-energy states and thus signals for the Bogoliubov transformation.

Once we have the wave-functions, we can construct the spectral functions for any t, i.e. a simple calculation leads to

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{A}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right) & = & \int _{-\Delta }^{\Delta }{\rm d}\delta t\;{{{\rm e}}^{{\rm i}{{p}_{0}}\delta t}}\left[ \psi _{{\boldsymbol{p}} }^{\left( + \right)}\left( t+\frac{1}{2}\delta t \right)\psi _{{\boldsymbol{p}} }^{\left( + \right)*}\left( t-\frac{1}{2}\delta t \right)-\psi _{{\boldsymbol{p}} }^{\left( - \right)}\left( t+\frac{1}{2}\delta t \right)\psi _{{\boldsymbol{p}} }^{\left( - \right)*} \right. \\ {} & {} & \quad \left. \times \left( t-\frac{1}{2}\delta t \right) \right], \\ \end{array}\end{eqnarray} \tag{ 25 }$$

for Pauli–Jordanʼs function and we can find a similar expression for Hadamardʼs function. The time evolution of the spectral functions may have dependence on the choice of $\Delta$. Intuitively, $\Delta$ corresponds to the observation time, as we already mentioned, to detect the quasi-particle behaviour in the oscillation pattern. For the concrete demonstration, let us take a look at figure 2 again; the temporal oscillation shows a constant pattern except near the origin where the system is disturbed by pulsed electric fields. So, around t = 30 for example, if $\Delta$ is less than 30, the quasi-particle behaviour is well separated from the interaction region at the origin and the spectral functions should be close to (15) then. If $\Delta$ is greater then 30, however, the integration region covers the pulsed electric fields, which should alter the spectral shape. Indeed, as we can see in figure 3, we can confirm this anticipation by comparing the results at $\Delta =25<t=30$ and $\Delta =50$ for ${{\mathcal{D}}^{\Delta }}-{{\mathcal{A}}^{\Delta }}$.

Zoom In Zoom Out Reset image size

Figure 3 already indicates the Schwinger process of the particle production. We can see a peak at ${{p}_{0}}\sim {{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)=1.5$ and its height corresponds to $|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}}$ according to (5). Precisely speaking, if we take $\Delta \to \infty$, the peak becomes a Diracʼs delta function and the ${{p}_{0}}$-integration in (5) has a contribution from a point ${{p}_{0}}={{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)$ only. Now that we implement the Wigner transformation with a finite $\Delta$, the peak is broadened and we should keep the ${{p}_{0}}$-integration over the range of the order of $1/\Delta$.

From (15) we can understand two more peaks are expected at ${{p}_{0}}\sim -{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)=-1.5$ and ${{p}_{0}}\sim 0$ which are evident in figure 3. For $\Delta =50$, we are also aware of some enhancement around ${{p}_{0}}\sim -1$ and 0.5, which is quite non-trivial. They arise from the effect of the background fields when $\Delta$ is comparable to or greater than t.

We are interested in the full temporal profile of these enhanced regions, so we make a density plot for ${{\mathcal{D}}^{\Delta }}$ and ${{\mathcal{A}}^{\Delta }}$ individually, as is shown in figure 4. This figure provides us with useful messages about the flow of the spectral peaks. First of all, the oscillatory peak at ${{p}_{0}}=0$ appears only in ${{\mathcal{D}}^{\Delta }}$ at late t, as is the case in the out-state (13). Second, we can notice that the intermediate enhancement around ${{p}_{0}}\sim -1$ and 0.5 emerges in both ${{\mathcal{D}}^{\Delta }}$ and ${{\mathcal{A}}^{\Delta }}$. Because the enhancement originates from time-dependent background fields, it would be conceivable that no simple pattern but complicated time dependence may well occur. This is not the case, however, and the enhancement goes rather straight in time. We can observe this in a clearer way in the form of not the density plot but the three-dimensional plot as presented in figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 5 we plot ${{\mathcal{D}}^{\Delta }}-{{\mathcal{A}}^{\Delta }}$ that is the difference between two in figure 4. Near the in-state, as seen in figure 5, there stands only one peak in the vicinity of ${{p}_{0}}=-{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$, which agrees perfectly with the asymptotic analysis (11). Then, this peak diminishes with increasing time, and meanwhile, the intermediate oscillatory modes grow up at ${{p}_{0}}\sim -1$ and 0.5. Eventually, these modes fade away, and at the same time, the spectral function approaches the asymptotic form of (15). At late time we can recognize a small peak around ${{p}_{0}}\sim 1.5$ and this peak amounts to the Schwinger mechanism. In other words, the Schwinger mechanism is a phenomenon that takes account of such a small portion of the whole spectral shape. It is certainly worth considering other tools to diagnose a wider region of the spectral functions, which is an interesting future problem beyond the present scope.

Finally we shall extend our analysis to a more general situation of the initial state. A more non-trivial but still controllable example is the finite temperature/density calculation of the Schwinger mechanism [33]. Let us assume that the initial state is a mixed state characterized by the following density matrix,

$$\begin{eqnarray}&&{{\hat{\rho }}_{\infty }}=\frac{{\rm exp} \left[ -\beta {{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)\hat{a}_{{\rm in}}^{\dagger }\left( {\boldsymbol{p}} \right){{{\hat{a}}}_{{\rm in}}}\left( {\boldsymbol{p}} \right) \right]}{{\rm tr}\left\{ {\rm exp} \left[ -\beta {{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)\hat{a}_{{\rm in}}^{\dagger }\left( {\boldsymbol{p}} \right){{{\hat{a}}}_{{\rm in}}}\left( {\boldsymbol{p}} \right) \right] \right\}},\end{eqnarray} \tag{ 26 }$$

with β being the inverse temperature. Then the straightforward calculation immediately leads to the produced particle number given as

$$\begin{eqnarray}&&\frac{{\rm d}N}{{{{\rm d}}^{3}}{\boldsymbol{p}} }-\frac{V{{f}_{{\boldsymbol{p}} }}}{{{\left( 2\pi \right)}^{3}}}=\frac{V|{{\beta }_{{\boldsymbol{p}} }}{{|}^{2}}}{{{\left( 2\pi \right)}^{3}}}\left( 1+{{f}_{{\boldsymbol{p}} }}+{{{\bar{f}}}_{-{\boldsymbol{p}} }} \right)\ ,\end{eqnarray} \tag{ 27 }$$

where ${{f}_{{\boldsymbol{p}} }}$ and ${{\bar{f}}_{{\boldsymbol{p}} }}$ represent the Bose–Einstein distribution function, respectively, for the particle and the anti-particle with the in-state energy ${{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$, namely,

$$\begin{eqnarray}&&{{f}_{{\boldsymbol{p}} }}:=\frac{1}{{{{\rm e}}^{\beta \left[ {{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)-\mu \right]}}-1}\ ,\qquad {{\bar{f}}_{{\boldsymbol{p}} }}:=\frac{1}{{{{\rm e}}^{\beta \left[ {{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)+\mu \right]}}-1}\end{eqnarray} \tag{ 28 }$$

with a chemical potential μ introduced. We note that, in the left-hand side of (27), the number of thermal particles is subtracted since they are irrelevant to the particle production. This result is understandable also from the spectral functions. We can find that Hadamardʼs function picks the Bose–Einstein distribution function up as

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{D}}^{\Delta }}\left( t,{{p}_{0}},{\boldsymbol{p}} \right) & = & \int _{-\Delta }^{\Delta }{\rm d}\delta t\;{{{\rm e}}^{{\rm i}{{p}_{0}}\delta t}}\left[ \psi _{{\boldsymbol{p}} }^{\left( + \right)}\left( t+\frac{1}{2}\delta t \right)\psi _{{\boldsymbol{p}} }^{\left( + \right)*}\left( t-\frac{1}{2}\delta t \right)\left( 1+2{{f}_{{\boldsymbol{p}} }} \right) \right. \\ {} & {} & \left. \qquad \qquad \qquad +\psi _{{\boldsymbol{p}} }^{\left( - \right)}\left( t+\frac{1}{2}\delta t \right)\psi _{{\boldsymbol{p}} }^{\left( - \right)*}\left( t-\frac{1}{2}\delta t \right)\left( 1+2{{{\bar{f}}}_{-{\boldsymbol{p}} }} \right) \right], \\ \end{array}\end{eqnarray} \tag{ 29 }$$

but Pauli–Jordanʼs function ${{\mathcal{A}}^{\Delta }}$ remains independent of the medium effect and is not changed from the vacuum expression (25). Such a qualitative difference between ${{\mathcal{D}}^{\Delta }}$ and ${{\mathcal{A}}^{\Delta }}$ makes a sharp contract and is not quite comprehensible on the intuitive level.

The particle production is increased by the Bose–Einstein distribution function. Therefore, if ${{f}_{{\boldsymbol{p}} }}$ takes a macroscopic value, this increase must be a sizable effect. In the high-T limit, in fact, the distribution function approaches, ${{f}_{{\boldsymbol{p}} }}\sim {{\bar{f}}_{{\boldsymbol{p}} }}\to T/{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$, and the particle production is significantly enhanced by $2\;T/{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$, which is a substantial factor if T is large or ${{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$ is small enough, as is the case in the quark-gluon plasma. However, the in-state already contains as many particles as ${{f}_{{\boldsymbol{p}} }}$ and so the particle production is not practically enhanced if measured relative to the number of particles in the in-state.

Another interesting limit lies in a finite chemical potential that makes ${{\bar{f}}_{{\boldsymbol{p}} }}$ (or ${{f}_{{\boldsymbol{p}} }}$) be much bigger than ${{f}_{{\boldsymbol{p}} }}$ (or ${{\bar{f}}_{{\boldsymbol{p}} }}$). This may well opens a new possibility for the experimental detection of the Schwinger process. For example, if we have a macroscopic occupation number like the Bose–Einstein condensate of scalar particles (that is actually a superconductor), the anti-particle (hole) that did not exist in the in-state is produced with a gigantic enhancement factor by the macroscopic occupation number. It may be worth pursuing this possibility further in the future research.

We found a useful formula that relates the particle production to the dynamical (time-dependent) spectral functions. We then clarified the basic properties of these spectral functions and proceeded to the numerical calculation of the spectral functions using the solutions of the equation of motion for a complex scalar field theory under pulsed electric fields. We closely studied the time evolution of the spectral functions. Wigner transformed Hadamardʼs function turned out to exhibit an oscillatory mode at ${{p}_{0}}=0$ as a result of the Bogoliubov mixing. This peak is larger by one power of the Bogoliubov coefficient as compared to the other peak corresponding to the Schwinger mechanism. This structure hints a new possibility of measurement that verifies the Bogoliubov mixing. Another non-trivial finding is the appearance of transient enhancement in the intermediate time region. In spite of time-dependent background fields, the enhancement occurs somehow in an organized manner, which implies that some unknown mechanism underlies the real-time dynamics. We also extended our discussions to the finite temperature/density case to identify the medium enhancement factor.

Apart from the Schwinger problem, from a more general perspective of theoretical physics, the dynamical properties of the spectral functions are quite non-trivial and are still less unknown than the equilibrated matter. We emphasize that these spectral functions are essential ingredients to think of real-time physics, and the particle production is actually one of the possible applications. In this sense we should continuously invest our efforts to deepen the theoretical understanding of the dynamical spectral functions and the present work should contribute to the first step along this direction.

This work was supported by JSPS KAKENHI grant number 24740169.

We first recall the standard argument of the LSZ formula without external fields, which connects the matrix elements of the S-matirx and the vacuum expectation values. The key equation is (2) but there is no distiction between the in- and the out-states without electric fields. Then, one can promote this expression at $t=\infty$ to a definition of $\hat{a}\left( t,{\boldsymbol{p}} \right)$ for any t, that is,

$$\begin{eqnarray}&&\sqrt{2E\left( {\boldsymbol{p}} \right)}\;\hat{a}\left( t,{\boldsymbol{p}} \right):={\rm i}\;{{{\rm e}}^{{\rm i}E\left( {\boldsymbol{p}} \right)t}}\left[ {{\partial }_{t}}-{\rm i}E\left( {\boldsymbol{p}} \right) \right]\;\hat{\phi }\left( t,{\boldsymbol{p}} \right).\end{eqnarray} \tag{ A.1 }$$

In principle, of course, the time-dependence in operator can be introduced by the Hamiltonian, and from this point of view, the above definition (A.1) corresponds to the interaction picture in which the free Hamiltonian evolves the system. In any case, what is important is only the asymptotic relation of (2) and, as long as (2) holds, the definition of $\hat{a}\left( t,{\boldsymbol{p}} \right)$ is irrelevant to the final answer.

From the surface terms as a result of the integration of the total derivative in t, we can easily express the asymptotic operator values as

$$\begin{eqnarray}\begin{array}{rcl} \sqrt{2E\left( {\boldsymbol{p}} \right)}\;\left[ \hat{a}\left( \infty ,{\boldsymbol{p}} \right)-\hat{a}\left( -\infty ,{\boldsymbol{p}} \right) \right] & = & \int _{-\infty }^{\infty }{\rm d}t\;{{\partial }_{t}}\;{\rm i}\;{{{\rm e}}^{{\rm i}E\left( {\boldsymbol{p}} \right)t}}\;\left[ {{\partial }_{t}}-{\rm i}E\left( {\boldsymbol{p}} \right) \right]\hat{\phi }\left( t,{\boldsymbol{p}} \right) \\ {} & = & {\rm i}\int _{-\infty }^{\infty }{\rm d}t\;{{{\rm e}}^{{\rm i}E\left( {\boldsymbol{p}} \right)t}}\;\left[ {{\partial }_{t}}+{\rm i}E\left( {\boldsymbol{p}} \right) \right]\left[ {{\partial }_{t}}-{\rm i}E\left( {\boldsymbol{p}} \right) \right]\hat{\phi }\left( t,{\boldsymbol{p}} \right) \\ {} & = & {\rm i}\int {{{\rm d}}^{4}}x\;{{{\rm e}}^{{\rm i}{{p}_{0}}t-{\rm i}{\boldsymbol{p}} \cdot x}}\;\left( \square +{{m}^{2}} \right)\hat{\phi }\left( t,x \right){{\left| {} \right.}_{{{p}_{0}}=E\left( {\boldsymbol{p}} \right)}}. \\ \end{array}\end{eqnarray} \tag{ A.2 }$$

In the computation of the S-matrix the time-ordered product arranges the ordering of $\hat{a}\left( \infty ,{\boldsymbol{p}} \right)$ and $\hat{a}\left( -\infty ,{\boldsymbol{p}} \right)$ properly, which is not required for the present problem of the expectation value of the number operator.

At this point it should be clear why such a manipulation needs some modifications when the time-dependent background fields are imposed. We should alter the definition (A.1), so that we can treat different ${{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$ and ${{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)$. One simplest prescription to do this is the following; in the adiabatic limit with slowly changing vector potential $A\left( t \right)$, we may be able to approximate the dispersion relation as $E\left( t,{\boldsymbol{p}} \right)=\sqrt{{{\left[ {\boldsymbol{p}} +eA\left( t \right) \right]}^{2}}+{{m}^{2}}}$ as assumed in [34, 35] (see also [36] for more discussions on the physical interpretation). Alternatively, one can interpret this $E\left( t,{\boldsymbol{p}} \right)$ as a definition rather than an approximation. In this way we can arrive at an expression similar to (A.2) as

$$\begin{eqnarray}\begin{array}{ll} {} & {} & \sqrt{2{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\;{{{\hat{a}}}_{{\rm out}}}\left( {\boldsymbol{p}} \right)-\sqrt{2{{E}_{{\rm in}}}\left( {\boldsymbol{p}} \right)}\;{{{\hat{a}}}_{{\rm in}}}\left( {\boldsymbol{p}} \right) \\ {} & {} & \qquad ={\rm i}\int _{-\infty }^{\infty }{\rm d}t\;{{{\rm e}}^{{\rm i}{{\int }^{t}}{\rm d}t^{\prime} \,E\left( {{t}^{^{\prime} }},{\boldsymbol{p}} \right)}}\;\left[ {{\partial }_{t}}+{\rm i}E\left( t,{\boldsymbol{p}} \right) \right]\left[ {{\partial }_{t}}-{\rm i}E\left( t,{\boldsymbol{p}} \right) \right]\;\hat{\phi }\left( t,{\boldsymbol{p}} \right). \\ \end{array}\end{eqnarray} \tag{ A.3 }$$

This expression does not reduce to the LSZ-type formula, however, as used in [16, 19]. Then, how is it possible to have $\left( \square +{{m}^{2}} \right)$ in front of $\hat{\phi }\left( t,x \right)$ as is the case in (A.2) but not in (A.3)?

The solution to this problem is quite simple. We should define ${{\hat{a}}_{{\rm out}}}\left( t,{\boldsymbol{p}} \right)$ from (A.1) with $E\left( t,{\boldsymbol{p}} \right)$ replaced with ${{E}_{{\rm out}}}\left( t,{\boldsymbol{p}} \right)$. If the boundary condition is chosen for $A\left( t \right)$ being zero at $t\to \infty$ (which is always possible in a certain gauge), we can find

$$\begin{eqnarray}&&\quad \sqrt{2{{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}\;\left[ {{{\hat{a}}}_{{\rm out}}}\left( {\boldsymbol{p}} \right)-{{{\hat{a}}}_{{\rm out}}}\left( -\infty ,{\boldsymbol{p}} \right) \right]={\rm i}\int {{{\rm d}}^{4}}x\;{{{\rm e}}^{{\rm i}{{p}_{0}}t-{\rm i}{\boldsymbol{p}} \cdot x}}\;\left( \square +{{m}^{2}} \right)\hat{\phi }\left( t,{\boldsymbol{p}} \right){{\left| {} \right.}_{{{p}_{0}}={{E}_{{\rm out}}}\left( {\boldsymbol{p}} \right)}}.\end{eqnarray} \tag{ A.4 }$$

With this definition we can confirm that the standard LSZ-type formula works; however, the price to pay is that ${{\hat{a}}_{{\rm out}}}\left( -\infty ,{\boldsymbol{p}} \right)$ is not the same as ${{\hat{a}}_{{\rm in}}}\left( {\boldsymbol{p}} \right)$. Therefore, we can correctly evaluate the spectrum, ${\rm d}N/{{{\rm d}}^{3}}{\boldsymbol{p}}$, at $t=\infty$, but we should adjust it to throw unphysical contributions away from $t=-\infty$. This is the reason why a subtraction was necessary in [19] to realize the empty vacuum that contains no particle at the initial time.

We emphasize here that the final results do not depend on this artificial definition of $\hat{a}\left( t,{\boldsymbol{p}} \right)$ because only the surface terms are physical. So, (A.3) as well as (A.4) work fine to yield the correct value of ${\rm d}N/{{{\rm d}}^{3}}{\boldsymbol{p}}$ at $t=\infty$. If we adopt our spectral formulas such as (5) and (6) to begin with, they refer to the surface terms only and there is no need to introduce $\hat{a}\left( t,{\boldsymbol{p}} \right)$. Yet, the non-uniqueness of the intermediate annihilation operator is attributed to the undetermined parameter $\Delta$ in the definition of spectral functions. Although the final answer is intact, a different $\Delta$ would lead to a different picture of quasi-particles near the initial and the final states.