Finite temperature Schwinger pair production in coexistent electric and magnetic fields

We compute Schwinger pair production rates at finite temperature, in the presence of homogeneous, concurrent electric and magnetic fields. Expressions are obtained using the semiclassical worldline instanton formalism, to leading order, for spin-$0$ and spin-$\frac{1}{2}$ particles. The derived results are valid for weak coupling and fields. We thereby extend previous seminal results in the literature, to coexistent electric and magnetic fields, and fermions.


I. INTRODUCTION
The non-perturbative pair production of electrically and magnetically charged particles in the background of large field strengths has garnered much interest and study over the years. Sauter [1], as well as Heisenberg and Euler [2], had speculated that sufficiently large electric fields could lead to spontaneous pair production of e + -e − . The notion was further sharpened and investigated comprehensively by Schwinger [3]; deriving the imaginary part of the QED one-loop effective action. These results were then further generalised by various authors to diverse cases -for instance, to extended objects such as magnetic monopoles [4], spatial or temporal inhomogeneous fields [5,6] and arbitrary gauge couplings [7,8], to cite a few examples (see for instance [9,10] and references therein for a more complete discussion). Exact analytic expressions are known nevertheless only for a few special cases and extending investigations into hitherto unexplored regimes is an ongoing endeavour.
The worldline path integral formalism has proven to be a potent method for perturbative and non-perturbative quantum field theoretic computations. The origins of the method may be traced to ideas by Fock [11], Nambu [12] and the Feynman worldline representation of one-loop effective actions [13,14]. The formalism was, for instance, leveraged to compute pair-production rates for magnetic monopoles at strong coupling [4,7]. With the development of string theoretic techniques towards understanding gauge theory scattering amplitudes [15][16][17][18][19], the method found further resurgence and applications (see for example [20] and related references); particularly, in our context, conveniently accommodating computations with large external fields [21][22][23][24][25].
Our aim in this work is to extend these results to the case when there are homogeneous (spatially and temporally) electric and magnetic fields simultaneously present. We compute leading order thermal corrections, using worldline path integral techniques, to the nonperturbative vacuum decay rates when there are coexistent electric and magnetic fields. We work in a regime where the coupling constant is small, and the external fields are also relatively weak. As far as we know, these expressions have not been computed before in the literature. We will largely follow techniques developed in [4-6, 20, 41]. In the limit of vanishing temperature (T → 0), one recovers the well-known results in literature [3,[43][44][45][46][47][48]. When the magnetic field vanishes (B → 0), in the case of scalar quantum electrodynamics (SQED), the results are seen to relapse into the known expressions for pure homogeneous electric fields, computed recently [41]. In quantum electrodynamics (QED), with fermions, we also obtain new expressions in the B → 0 limit that complement these recent SQED results.
It is well known that even at zero temperature (T = 0), the presence of a magnetic field parallel to the electric field (E B), leads to interesting modifications to vacuum decay rates, relative to the pure electric field case. The vacuum decay rates, per unit volume, at T = 0 for homogenous E B are given by [3,[43][44][45][46][47][48] Here, m and q are the mass and electric charge of the particle under consideration. Note that in addition to the usual enhancement due to extra degrees of freedom in the spin-1 2 case, the vacuum decay rates in the fermion case may be further enhanced, relative to the scalar case, when B > E.

arXiv:1808.01295v3 [hep-th] 29 Mar 2019
Note also that, for any homogeneous E and B fields, for which the Lorentz invariant E · B = 0, one may go to a frame of reference with boost ( υ) given by [49] where the transformed fields ( E and B) are parallel to each other. This is the so-called centre-of-field frame.
Since the vacuum decay rate per unit volume is a Lorentz invariant, one may conveniently compute it in this centreof-field frame. The formulas for homogeneous E B are therefore potentially of wide applicability. For homogeneous fields with E · B = 0, but the fields not equal in magnitude, a reference frame may be found where the transformed field is purely electric or magnetic [49]. In this latter scenario, the relevant expressions are those of single field Schwinger pair production. Apart from being of significant theoretical interest, scenarios with parallel electric and magnetic fields are also relevant in various astrophysical systems. For instance, it is believed that neutron stars such as pulsars have strong electrical fields parallel to the magnetic field in their polar vacuum gap regions [50]. Neutron star surface temperatures are expected to reach ∼ 10 5 K. Non perturbative production of exotic states such as millicharged particles, which may form a component of dark matter, may occur in these vacuum gap regions and provide hitherto unknown constraints on these states [51] (in the context of constraints from non-perturbative production, in pure E or B fields, also see [52] for millicharged particle bounds from accelerator cavities, and [53,54] for bounds on magnetic monopoles). These settings also, therefore, make the results phenomenologically very relevant.
In Sec. II, we discuss the derivation in the case of SQED. Towards the exposition of necessary techniques and to fix notations, we re-derive the known zero temperature result for the case of E B using worldline instantons, before presenting the main result for finite temperature. Then, in Sec. III, we consider QED. Results are presented for spin-1 2 particles in the zero temperature and finite temperature cases. The finite temperature SQED and QED Schwinger pair production results, for E B, are new and readily generalise earlier seminal results in the literature [8,40,41]. Even for the zero temperature cases, to the best of our knowledge, this is the first time that an explicit and complete derivation is being presented for vacuum decay rates, when E B, using worldline instanton techniques. We summarise our main results, shortcomings of the derivations and future directions in Sec. IV.

II. THERMAL PAIR PRODUCTION FOR E B IN SQED
We would like to calculate decay rates for vacua, made metastable by the presence of large external fields. Let us denote the probability for vacuum to vacuum transitions by 0 out |0 in 2 . In presence of external fields sourced by a potential A, the probability for vacuum to vacuum transitions are given by where W M [A] is the Minkowskian effective action for the theory under consideration. Expressed in terms of Euclidean quantities, specialising now to Scalar electrodynamics (SQED), we have explicitly Here, In the above Euclidean expressions, the covariant derivative D µ := (∂ µ + iqA µ ), external gauge field A µ := (A 1 , A 2 , A 3 , A 4 ) with A 4 := −iA 0 , and field tensor In terms of Euclidean quantities, the Lorentz invariant vacuum decay rate per unit volume is given by We simplify the effective action further, following a standard technique [5,7], and after performing a functional integration, one obtains Using Frullani's integral identity Tr(ln M) [3,55], dropping terms that do not contribute to the imaginary part, and converting the trace to a path integral, leads then to the well-known expression for the SQED one-loop Euclidean effective action [5,7] Here, x denotes differentiation with respect to τ . An implicit assumption that the coupling constant is small (q 2 1) has been made while writing the above result, by dropping non-local interaction terms that are higher order in the coupling constant. Now, making a substitution τ = z u and z → z/m 2 gives x denotes differentiation with respect to u. Evaluating the z integral above, gives where K 0 (z) is the modified Bessel function of the second kind. For x 1, we have the asymptotic formula K 0 (x) ∼ exp(−x) and hence, when m 1 0ẋ 2 du 1, the above expression may be simplified to The assumption m 1 0ẋ 2 1 is equivalent to making a weak field approximation qE/m 2 2π [5,7]. We will therefore also assume that the external electromagnetic fields are relatively weak and satisfy q|F |/m 2 1, for field strengths |F |. Finally, note that Eq. (11) may equivalently be obtained by making a saddle point approximation, to the z integral in Eq. (9).
Considering the terms in the exponent as part of an effective action, the corresponding Euler-Lagrange equations are The antisymmetry of F µν immediately implies thaṫ where ρ is a constant. Specialising now to temporally and spatially homogeneous E B, let us choose E and B in the x 3 direction, without loss of generality. We then have for the Field tensorF , This leads to the equations of motion To clarify ideas and general techniques, that shall be adopted in the finite temperature derivation, we first derive the well-known result in the T = 0 case, using the worldline instanton formalism. Though this result is well-known and has been derived using many other techniques [43][44][45][46][47][48], we believe that a systematic derivation of this has not been presented before in the literature, using worldline path integral methods.
We note from Eq. (16) that the equations of motion for x 1 , x 2 and x 3 , x 4 are decoupled from each other. The set of equations for x 1 , x 2 give rise to hyperbolic solutions, which fail to satisfy the periodic boundary condition x µ (0) = x µ (1), as required by Eq. (9). Thus, the only solutions for x 1 and x 2 are trivial solutions. For x 3 and x 4 one finds solutions satisfying the required periodic boundary conditions. Let us collectively denote these solutions byx. Note that in the above, one must have ρ = 2πkR = m2kπ qE , to satisfy the boundary conditions. These solutions therefore represent a circle in the x 3 − x 4 plane, with radius R = m/qE. This is equivalent to the situation in the pure E case [7]. The effective action, with these solutions (x), is then given by Let us now compute the fluctuation prefactor for this solution(for general techniques, see for instance [56][57][58]). To leading order, the fluctuation prefactor is proportional to ∼ det[δ 2 S eff /δx ν δx µ ] −1/2 , evaluated at the solutions to the equations of motion, with appropriate boundary conditions.
Define the prefactor matrix at zero temperature . The relevant determinant, with zero modes removed, may be expressed using the matrix determinant lemma (see for instance [41,59]) as is to be interpreted as a Green's function. In the E B case, we have Two of the eigenvalues are -2πqE(l 2 /k − l), corresponding to eigenvectors (0, 0, cos(2lπu), sin(2lπu)) and (0, 0, sin(2lπu), − cos(2lπu)), and 2πqE(l 2 /k + l), corresponding to eigenvectors (0, 0, sin(2lπu), cos(2lπu)) and (0, 0, cos(2lπu), − sin(2lπu)). The other two eigenvalues have the form -2πqE − iBl E + l 2 k , corresponding to eigenvectors ( . In all cases l runs from 1 to ∞. With these, one obtains where N 0 = m 2 kπ/qE. The infinite products may be simplified [55], and one obtains the compact expression It is interesting to compare this to the equivalent expression in the case of pure E [6]. The only part remaining to be calculated is the non-local factor that appears in Eq. Here, for the nontrivial solutions, the only relevant part of C 0 is the (3−4) block. The Green's function G 0 (u, u ) can be obtained in the standard way by constructing a spectral representation. Utilising the relation C −1 with V the eigenvector column matrix and C 0,D the diagonal matrix, one gets With the non-trivial solutions for x 3 and x 4 this gives Therefore, due to the decoupling in Eq. (16) leading to trivial solutions for x 1 and x 2 , the non-local part of the prefactor matrix determinant comes out to be unity, in complete analogy to the pure E case [41].
Putting all the factors together, the fluctuation prefactor for fixed k finally comes out to be The relevant part of the SQED Euclidean effective action then becomes (27) From this, using Eq. (6), the T = 0 vacuum decay rate per unit volume, in SQED for homogeneous E B, may be calculated finally as This matches the well-known zero temperature SQED expression in literature [3,[43][44][45][46][47][48], as given in Eq. (1). Note that it also reduces to the pure E case in the limit B → 0, as expected.
With this warm-up derivation in the zero temperature case, clarifying ideas and techniques, we now proceed to thermal Schwinger pair production in SQED when one has homogeneous E B fields. For calculating finite temperature vacuum decay rates, for scalar particles in the presence of a homogeneous electromagnetic field, we need to calculate the imaginary part of the SQED thermal effective action. The supplemental requirement in the thermal case is that the Euclidean time direction must now be compact with endpoints identified and one requires Here, β −1 is the temperature (T ), that is assumed to be much less than the mass of the particle under consideration (T m). The SQED Euclidean effective action at finite temperature is given by (29) One has again assumed weak fields, qF/m 2 1, and small couplings q 2 1. Note that n = 0 coincides with the expression already derived, for zero temperature. We focus on the n = 0 contributions. The terms in the exponent above, are again to be considered as part of some effective action (S eff ).
To find the relevant thermal instantons, we need to find solutions to the equations of motion Eq. (16), that are now additionally compact in x 4 , with period nβ [56][57][58][60][61][62]. Thus, we have to essentially find local sections of the zero temperature instanton solutions Eq. (17), that are additionally periodic by nβ in the x 4 direction. For such viable solutions to exist, we must have 2R ≥ nβ, as is clear from geometry. This implies a bound n max = 2R/β , where x denotes the integer less than or equal to x. This means that there are no one-loop thermal contributions for T < qE/2m ≡ T * , defining a critical temperature T * for a given mass m and charge q. Since there are no thermal corrections below T = T * , it may provide a partial resolution with some earlier studies [27,33,34], where it was argued that there are no thermal corrections at one-loop (also see discussions in [41,42]). Now, for n ∈ Z − , i.e. solutions satisfying the boundary condition there are two solutions (see Fig. 1). For the smaller path (I − ), subtending angle θ n at the center, Θ = 2πk + θ n is the total angle subtended by k windings. The explicit solution (x T,I − ) in this case is given by with the end-points of x 4 identified. There are again no non-trivial solution for x 1 and x 2 satisfying the requisite periodic boundary conditions, similar to the zero temperature case. The corresponding effective action may be computed for this solution, from Eq. (12), and comes out to be where R = m/qE and T * = qE/2m. The relation between angle subtended θ n and temperature T is As we shall see, a calculation of the fluctuation prefactor for the I − solution (x T,I − ) shows that it does not contribute to the imaginary part of the Euclidean effective action. Therefore, the solutionx T,I − may only contribute to the free energy, and there is no contribution to the vacuum decay rate from it.
(34) The corresponding effective action, using Eq. (12), may be calculated and gives where, as before, T * = qE/2m and R = m/qE. The relation between θ n and n for II − , is same as in Eq. (33). This solution, as we shall demonstrate while calculating the fluctuation prefactor, will be one that does contribute to the vacuum decay rate, by giving an imaginary part to the Euclidean effective action.
For the positive integer case, n ∈ Z + case, we have the requirement There are again two solutions (see Fig. 2). For the smaller path (I + ), subtending angle θ n at the center, Θ = 2πk + θ n is the total angle subtended. k is again the number of windings. As is amply clear from Fig. 2 and geometry, the explicit solution (x T,I + ) for this case is These solutions give for the effective action and as in the n ∈ Z − case the computation of prefactor shows that it only contributes to the free energy and not to pair production. Coming now to the longer path (II + ), subtending an angle 2π − θ n at the center (see Fig. 2) we have for kwindings, a total angle subtended Θ = 2π(k+1)−θ n . The II + solution, similar to II − of Eq. (34), will contribute to the imaginary part of the effective action. This solution (x T,II + ) is explicitly x 3 (u) = R cos(Θu+π+θ n /2) , x 4 (u) = R sin(Θu+π+θ n /2) .
FIG. 2. The I + and II + solutions corresponding to n ∈ Z + . The short path I + , as in the earlier case, does not contribute an imaginary part to the Euclidean effective action. The long path II + does contribute an additional negative eignevalue from the non-local part in the prefactor matrix, and hence contributes to an imaginary part for the Euclidean effective action. It is thus thex T,II + solution again that would contribute to vacuum decay rates we are interested in.
The corresponding effective action is calculated to be T * , R, are as defined earlier and the relation between θ n and n is now Note from Eq. (35) and Eq. (40) that the two solutions, x T,II − andx T,II + contributing to the vacuum decay rate, actually give equivalent expressions for the exponential factor. The contribution to pre-exponential factors will also be seen to be similar, for both solutions. Hence, the full sum over n ∈ Z may be replaced just by twice sum over n ∈ Z + . Hence, from now on, we will just consider the solutionx T,II + for presenting the relevant calculations.
Let us now compute the fluctuation prefactor relevant to thex T,II + solution. Again, define a prefactor matrix The relevant determinant, with the zero modes removed, may be written as [59] det Here, C −1 T := G T (u, u ) is again to be interpreted as an appropriate Green's function, without zero modes. The matrix C T is given in this case by Note the presence of additional elements in the 1-2 block, depending on magnetic field strength B, compared to the equivalent matrix in the pure electric field case [41]. For calculating det [C T ], we utilise the result [57,[63][64][65] det where C T is the matrix formed from C T by excluding all non-diagonal terms. ξ α andξ α are the eigenvalues of C T and C T respectively. The matrices ζ Since the eigen spectrum for C T is unknown, we may use the second equality in terms of the ζ andζ matrices to calculate det [C T ] [57,[63][64][65]. For the homegenous E B case we are considering, the coresponding ζ matrix, with appropriate boundary conditions, comes out to be From this, we have at u = 1 the determinant, . This is always positive. Theζ(u) matrix comes out to be -u · 1 4×4 . Hence, the relevant determinant is just det[ζ(1)] = 1. Putting all the above results together and using Eq. (45), det [C T ] may now be readily computed in our case as N T is a normalization factor that may be fixed explicitly by considering C T and the free theory. The factor (−1) k is related to the Morse index [58,64,66] of the corresponding solution.
The non-local part of the prefactor matrix determinant in Eq. (43), is of the form G T µν are Green's functions satisfying with vanishing Dirichlet boundary conditions. Since x 1 and x 2 do not have non-trivial solutions, satisfying required boundary conditions, the combinations containing G T 11 , G T 12 , G T 21 , and G T 22 in the integral all trivially give zero. The remaining terms are those with G T 33 , G T 34 , G T 43 and G T 44 . These are related to the 3-3, 3-4 and 4-4 elements of C T , which only depend on the electric field E. This immediately suggests that the Green's function should match that computed in the pure E, thermal case [41]. Since (C T ) 33 = (C T ) 44 and (C T ) 43 = −(C T ) 34 , it may be shown that G T 33 = G T 44 and G T 43 = −G T 34 . Solving Eq. (50), considering cases u > u and u < u , give These are in agreement with the expressions found in [41], for the pure E case. Putting all the above results together, the contribution of the non-local part, for T = 0 and homegeneous E B, come out to be This is manifestly negative, giving an extra negative mode for longer paths (II ± ), and thus contributing an imaginary part to the Euclidean effective action W E,E B T =0,scalar . Note that this is because the fluctuation prefactor is proportional to ∼ det[δ 2 S eff /δx ν δx µ ] −1/2 , evaluated at the stationary solutions. As alluded to before, the longer path solutions, therefore, contribute to vacuum decay rates. In contrast, substituting Θ corresponding to the shorter paths (I ± ), in place of Θ, would give a nonlocal contribution which is positive. This finally makes the fluctuation prefactor real and hence contributes only to free energy. In the pure electric field case, this was checked by matching the E → 0 limit of the short-path expressions [41], with the exact free energy density of a non-interacting relativistic particle [67], when β → ∞. A derivation of the free energy density using the standard proper time representation of the effective potential [34], in an external electric field, also matches that derived from the short-path expression. There is nevertheless some disagreement in the literature regarding the appropriate choice of path [8,[40][41][42].
Taking all the contributions into account, the thermal SQED fluctuation prefactor, for fixed k, comes out to be Finally, combining all the exponential and preexponential factors, the thermal vacuum decay rate per unit unit volume, to leading order, in the background of homogeneous E B comes out to be with Γ E B T =0,scalar given by Eq. (28) and Γ E B T,scalar by In above, n max = 2R/β , H(x) is the Heaviside step function, and Θ = 2π(k + 1) − θ n = 2π(k + 1) − 2 arcsin( nT * T ). In the limit of B → 0, Γ E B T =0,scalar reduces to the known expression for Γ E T =0,scalar [41]. Also, note that when T < T * ≡ qE/2m, the periodic boundary conditions on x 4 cannot be satisfied and there are no thermal corrections. In this case the result relapses to the zero temperature expression.

III. THERMAL PAIR PRODUCTION FOR E B IN QED
We now proceed to compute the non-perturbative pairproduction rates for spin-1 2 particles in Quantum Electrodynamics (QED). The derivation is analogous to the SQED derivation, with some subtleties coming from the additional Pauli spin term and the necessities of fermionic functional integrations.
For QED, the Euclidean effective action for fermion field Ψ is given by with Here, we define / are Euclidean gamma matrices, which are related to the Minkowskian gamma matrices through the relations in our convention. They satisfy For brevity, henceforth we will remove the subscript (E) from the Euclidean gamma matrices. Performing the fermion functional integral gives where σ µν = − i 2 [γ µ , γ ν ]. Using the Frullani integral identity [3,55], this may be expressed as (61) Note the additional factor of 1/2 compared to the scalar case as well as the additional Pauli spin term. These lead to interesting differences between the SQED and QED results. Introducing fermionic coherent states [20,68,69] and simplifying, the above Euclidean effective action may be re-written as In above, Tr f denotes a fermionic trace and we have assumed q 2 1. In terms of fermionic coherent states (η), this one-loop QED Euclidean effective action explicitly takes the form [20] Let us define and also note that d E is the number of Euclidean dimensions. Following a standard technique [70], let us then re-write the Eu-clidean effective action as As before,F is the electromagnetic field tensor with components F 12 = −F 21 = B and F 34 = −F 43 = iE. Note that det 1/2 (1 − 2iqF ( d dτ ) −1 ) should be a Lorentz scalar. Hence it should depend onF 2 and hence only on the coupling constant as q 2 [70]. Thus, we may relate . From this, we can write Using these definitions, In the case of interest, we havē SinceF 2 is diagonal, the factor Z may be evaluated readily as From this, we find The above determinant may be obtained in the usual way, by solving the eigenvalue problem with anti-periodic boundary condition f (z) = −f (0). The eigenfunctions satisfying these boundary conditions are The corresponding eigenvalues are given by Substituting this in Eq. (71), and taking into account the two-fold degeneracy, we get after a simplification of the infinite products [55], = cosh 2 (qBz) cos 2 (qEz) .
Substituting these results back, one obtains We now make a change of variable τ → zu, z → z/m 2 , as before, and perform the z integral using a saddle point approximation. The additional cosine term, in the fermion case above, gives an imaginary part in the exponential and hence does not modify the saddle point [5]. Also, the hyperbolic cosine term when written in its exponential form contributes a factor ±qBz/m 2 to the integrand's exponent. In the limit of weak fields, q|F |/m 2 1, this does not modify the saddle point either. Hence the saddle point for the z integral turns out to be z 0 = m 2 ( 1 0ẋ 2 du) 1/2 . The relevant one-loop Euclidean effective action in QED is then given by Considering the terms in the exponent above as part of an effective action, the corresponding Euler-Lagrange equations are again given by Let us initially consider the T = 0 case, as before. For E B, E and B assumed to be in the x 3 direction, there are again no non-trivial solutions for x 1 and x 2 , satisfying the periodic boundary conditions. Hence, in complete analogy to SQED, the only non-trivial solutions are This leads to the effective action and the exponential part of the QED vacuum decay rate The additional factors, for z 0 = mρ/2 = m 2 kπ/qE, come out to be cos(qEz 0 /m 2 ) cosh(qBz 0 /m 2 ) → (−1) k cosh(kπB/E) .
Let us now turn to the finite temperature case (T = 0). The computation follows the zero temperature case largely, with few additional complexities introduced by the requirement of the periodicity criteria along x 4 [56][57][58][60][61][62], as in SQED. To compute the fermion pair production at finite temperature, we again must consider solutions that are compact in the x 4 direction, with endpoints identified, and separated by nβ. Based on Eq. (62) and Eq. (76), the one-loop effective action for fermions at finite temperature is For E B, in the qF/m 2 1 regime, the equations of motion do not change compared to the corresponding scalar case. Hence, nor does the value of S eff (x T,II + ), computed earlier in Eq. (40). This leads to an exponent with S eff (x T,II + ), which using z 0 = mRΘ/2 for T = 0 leads to a factor exp(−S eff (x T,II + )) cos qEz 0 m 2 cosh The determinant of the prefactor matrix det [P T,fermion ], which appears in the computation of the QED fluctuation prefactor, also mostly remains the same as in the thermal SQED case. The relevant fluctuation prefactor hence becomes = −2 · (−1) k iV 3 β 4 . Combining all the above results, the leading order QED vacuum decay rate, per unit volume at finite temperature, in the background of coexistent, homogeneous electric and magnetic fields, is given by Here, as before, n max = 2R/β , H(x) is the Heaviside function, T * ≡ qE/2m, and Θ = 2π(k + 1) − θ n = 2π(k + 1) − 2 arcsin( nT * T ). In the limit B → 0, Γ E B T =0,fermion reduces to Γ E T =0,fermion and one obtains the leading order thermal corrections in QED for the pure E case, thereby complementing the known result for scalars [41]. For T < T * , again there are no thermal corrections and the result reverts to the T = 0 expressions.

IV. SUMMARY
The worldline path integral formalism provides a powerful and systematic way to compute nonperturbative vacuum decay rates in various situations. In this work, we computed leading order thermal corrections to vacuum decay rates, in SQED and QED, for the case of homogeneous, coexistent electric and magnetic fields. Apart from its theoretical importance, the results are relevant in astrophysical settings where large electric and magnetic fields may coexist in a thermal environment.
There are a few natural avenues to follow up on that were outside the scope of the present study. The Gaussian approximation to the fluctuation prefactor is inadequate, leading to spurious singularities at thermal thresholds, and one should include higher order terms to potentially mitigate this. This is challenging even in the zero temperature case, but based on the hard thermal loop framework [71,72], it has been argued that such spurious singularities may be softened and the result correctly interpreted [41]. Explicit calculation of these higher order terms beyond the Gaussian approximation would shed more light on the analytic structure of the terms at these thresholds. Another subtle point to note is that, even at zero temperature, the vacuum decay rate is not technically the same as the average, particle pair production rate [73,74]. In the zero temperature case, it may be shown that the physical observable-the mean pair production rate-is just the first term in the series for the vacuum decay rate [73]. Hence, for weak fields, the distinction is mostly pedantic. The thermal vacuum decay rates we compute are therefore expected to closely match the actual particle pair production rates for weak fields, but a more careful calculation is required to make the correspondence clear and rigorous. It would also be appealing to have a better physical understanding of the various results and reach a consensus on the remaining disagreements in the literature [8,[40][41][42]. Doing away with the assumption of relatively weak fields and extending the study to arbitrary coupling strengths would also be pertinent, as well as incorporating modifications due to field inhomogeneities.