Electric Field Decay Without Pair Production: Lattice, Bosonization and Novel Worldline Instantons

Electric fields can spontaneously decay via the Schwinger effect, the nucleation of a charged particle-anti particle pair separated by a critical distance $d$. What happens if the available distance is smaller than $d$? Previous work on this question has produced contradictory results. Here, we study the quantum evolution of electric fields when the field points in a compact direction with circumference $L<d$ using the massive Schwinger model, quantum electrodynamics in one space dimension with massive charged fermions. We uncover a new and previously unknown set of instantons that result in novel physics that disagrees with all previous estimates. In parameter regimes where the field value can be well-defined in the quantum theory, generic initial fields $E$ are in fact stable and do not decay, while initial values that are quantized in half-integer units of the charge $E = (k/2) g$ with $k\in \mathbb Z$ oscillate in time from $+(k/2) g$ to $-(k/2) g$, with exponentially small probability of ever taking any other value. We verify our results with four distinct techniques: numerically by measuring the decay directly in Lorentzian time on the lattice, numerically using the spectrum of the Hamiltonian, numerically and semi-analytically using the bosonized description of the Schwinger model, and analytically via our instanton estimate.


Introduction
It is said that in our quantum world, the totalitarian principle applies: everything not forbidden is compulsory. In quantum electrodynamics, non-zero electric fields can decay via the quantum nucleation of charged particles. This process has been very well-studied since the seminal work of Schwinger [1]. When the electric field is below the critical field E c ∼ m 2 /g (where m, g are the mass and charge of the particles) this process is nonperturbative and exponentially slow, with a rate that scales as e −πm 2 /gE . The particleanti-particle pairs responsible for this effect nucleate a distance d = 2m/gE apart and are separated along the direction the field is pointing. The instanton responsible for this process is simply an oriented circle with radius d (Figure 1a).
Many generalizations of [1] have been considered over the years, incorporating inhomogeneities or finite extent in the electric field. Without introducing inhomogeneities, another generalization is to consider the decay of an initial field when the direction along which it points is periodically identified. 1 If the circumference L of this compact dimension is larger than d one expects the standard instanton in Figure 1a to govern the decay and the pair-production rate to be approximately unchanged. However, when L < d the instanton no longer "fits" in the compact space.
In the regime L < d there is a stark disagreement in the recent literature over the decay rate and how to calculate it. In [2], Brown argues that the instanton responsible for the decay is similar to the standard circular instanton, but squeezed into a "lemon" so as to fit in the compact dimension ( Figure 1b). By contrast in [3], Medina and Ogilvie argue for a configuration that includes a different set of arcs of the circular instanton ( Figure 1d). The action for these configurations scales very differently from Brown's and is discontinuous as a function of L. The fact that these two estimates differ parametrically in the exponent is remarkable for such a simple and classic problem.
In this work we study this question in a simple case, namely the massive Schwinger model, quantum electrodynamics with massive charged fermions in one space and one time dimension. Neither of the "instantons" proposed in [2,3] actually solve the equations of motion of the theory (due to the discontinuous first derivative at certain points in the trajectory) and so neither should necessarily be expected correctly describe the decay. Here, we identify a novel set of instantons that do solve the equations of motion. The resulting prediction for the rate differs (in the exponent) from both [2,3] (although it is closer to [2]). We check our prediction using multiple numerical techniques, including directly time evolving the initial state in a lattice version of the theory and identifying the physical properties of the eigenstates of the lattice Hamiltonian. The massive Schwinger model has several features that distinguish it from QED in higher dimensions. Most importantly, there is no magnetic field and the gauge field is non-propagating, so the dynamics of the theory is entirely due to the charged particles (that interact through the electric fields they produce). Classically, this means that the electric field is almost entirely determined by the charge configuration (almost, because a constant background field is allowed). In the quantum theory there is a constraint that, in a charge eigenstate, determines the field (see (2.3)).
In particular, the possible field values are quantized in integer multiples of the charge g plus the constant background. One implication is that unlike for larger values, an initial field in the range −g/2 < E < g/2 cannot decay, because the energy density E 2 /2 would be larger were the field to increase or decrease by g. Indeed, one can think of a field in this range as a parameter defining the specific theory under study. This parameter is the famous θ angle, so-called because as the field varies the spectrum of the Hamiltonian undergoes a spectral flow that is periodic under E → E + g.
A novel and surprising feature of our results is that when the theory is compactified on a small circle there is a very sharp distinction between generic values of the initial field and those equal to integer or half-integer multiples of the charge, E = (k/2)g with k ∈ Z. In the regime of would-be slow decay and when E is not equal to one of these quantized values, the field never decays, in contrast to the non-compact or large-field case. This can be understood in our analysis from the fact that we do not find any instanton solution for such cases. We perform multiple distinct numerical analyses that confirm that both the quantum expectation value of the field and its (small) variance remain very close to their initial values for arbitrarily long times. By contrast, when E = (k/2)g the field oscillates coherently and sinusoidally in time between its initial value (k/2)g and −(k/2)g, with exponentially suppressed probability of taking any other value no matter how long one waits.
A feature of the massive Schwinger model is the existence of a bosonized version of the theory, the massive Sine-Gordon model. This constitutes a type of strong-weak duality, because the fermionic field theory is weakly coupled when g/m 1, while the bosonized description approaches a free field theory in the opposite limit g/m → ∞. The bosonized description gives an alternative view of electric field decay that can help account for the surprising features mentioned above. Meta-stable values of the field correspond to states in which the boson is localized in one of the local minima of the cosine potential. In the non-compact theory the field can decay via 1+1 dimensional bubble nucleation, where the walls of the bubble are a kink and an anti-kink (corresponding to a fermion and an anti-fermion). However, if the theory is compactified on a sufficiently small circle there is not enough room for a critical bubble to form. When the Kaluza-Klein modes can be ignored, the theory reduces to the quantum mechanics of the bosonic zero mode. The distinction between half integer and non-half integer E can then be understood as follows. For half-integer E the quantum mechanical potential has a Z 2 reflection symmetry, so that for each local minimum at E = (k/2)g there is a symmetric minimum at E = −(k/2)g. In the semi-classical regime the energy eigenstates are symmetric and anti-symmetric states localized in these two wells. Hence (much like the classic symmetric double-well problem in quantum mechanics) the field oscillates sinusoidally with an exponentially long period between the values E = ±(k/2)g. By contrast, for non-half integer values of the field the potential is not symmetric and the energy eigenstates are localized in a single well, so the field remains constant and localized for all time.
In summary, we arrive at our conclusions with the following set of distinct and complementary techniques: • We directly measure the time evolution of the electric field using a real-time lattice code for the fermionic massive Schwinger model • We compute the Hamiltonian eigenvalues from the lattice code, identify the states relevant to the electric field evolution, and estimate the rate of oscillation (when it exists) via energy differences • We identify the bosonized counterpart of the massive Schwinger model and compute the time-dependence of the electric field numerically and semi-analytically in the bosonic theory (using the approximation where the circle is small enough to ignore Kaluza-Klein modes).
• We estimate the rate using a novel set of field-theory instantons. 2 Relation to previous work: As already mentioned above, Brown [2] considered QED with a compact spatial or Euclidean time dimension and argued that the instanton illustrated in Figure 1b computes the rate of decay of the electric field. In contrast, Medina and Ogilvie [3] proposed a different set of instantons with a different action and hence an exponentially different rate (Figure 1d). Korwar and Thalapillil [4] considered the case of non-zero E and B fields; their results reduce to those of [3] when B = 0. Draper [5] also computed the rate, relying on the instanton illustrated in Figure 1c. None of these analyses coincide with our main results, although their focus is on higher-dimensional versions of QED. Qiu and Sorbo in [6] considered scalar QED in 1+1 dimensions in the linearized approximation around a background field. Their results also do not coincide with ours, although this may be due to the fact that they considered scalar QED and focused on the short-time behavior, whereas we consider the fermionic model and are mostly concerned with the evolution at long times. Electric field decay in the Schwinger model was studied numerically in e.g. [7,8]. Finally, [9] investigated the phenomenon of "flux unwinding" in the compact Schwinger model where a pair of charges can discharge many units of flux by repeatedly traversing the circle. Our present analysis differs from that of [9] because we focus on the regime where the circle is too small to accommodate pair production. The structure of this work is as follows. In Section 2, we describe the massive Schwinger model in the continuum and on the lattice. In Section 3, we outline the numerical techniques used on the lattice and present our results. In Section 4, we study the bosonized description of the theory in the small circle regime and numerically compute the transition rates. In Section 5 we use the worldline formalism of path integrals to identify a novel set of worldline instantons and compute the associated rates. We conclude in Section 6. In Appendix A, we review the lattice formulation of the massive Schwinger model. In Appendix B, we report technical details about the Hilbert space cutoff on lattice and the continuum extrapolation protocol for the energy spectrum in our numerics. In Appendix C, we discuss the continuum extrapolation of the prefactor c appearing in the bosonized action in Section 4. Appendix D contains calculations that complement Section 5.
(a) Schwinger's worldline instanton for pair production in [1,10], on a non-compact or compact spatial dimension of length L larger than the nu- WKB instanton on a compact spatial dimension in [2], which can be adapted to the finite temperature case.
x t E L (c) Worldline instanton on a compact spatial dimension in [5].

Massive Schwinger model on a compact spatial dimension
The massive Schwinger model [1,11] is quantum electrodynamics with a massive Dirac fermion in one space and one time dimension. The Lagrangian is where ψ is a two-component Dirac fermion, m is the fermion mass, D µ ≡ ∂ µ + igA µ is the covariant derivative with g being the electric charge and A µ being a U (1) gauge field, and F µν = ∂ µ A ν − ∂ ν A µ is the field strength. Throughout this work, we work in natural units = c = 1 with Minkowski metric η µν = diag(+1, −1) where the indices µ, ν runs from 0 to 1. In the temporal gauge A 0 = 0 such that the electric field E = F 10 = −Ȧ 1 (t), the Hamiltonian reads We are interested in this model on a compact spatial dimension of size L, i.e. with identification x ∼ x + L. We impose periodic boundary condition on both the electric field and the fermion, E(x + L) = E(x) and ψ(x + L) = +ψ(x).
Gauss' law ∂ 1 E = gj 0 implies that the electric field is given by where j 0 =ψγ 0 ψ and F is a constant background field which is physically significant in 1 + 1 dimensions [11]. Physical states satisfying Gauss' law have zero net charge on the circle. When the space is non-compact, the constant F is naturally fixed by the electric field value at infinity. On a circle of circumference L, we pick an arbitrary point to be the origin, and take it as the lower limit of the integral in (2.3); then F = E(x = 0) = E(x = L).

Massive Schwinger model on a lattice
Unlike its massless counterpart, the massive Schwinger model is not exactly solvable (although as we will see, its bosonized form is weakly coupled when the dimensionless coupling g/m is large). For this reason it is useful to simulate the model on a lattice and compare numerical results against analytical predictions. In this work, we perform numerical studies on a lattice using the Kogut-Susskind formulation of lattice gauge theory [12][13][14], which we review in Appendix A. Here we summarize the key elements of the construction. We put the model defined by the Hamiltonian (2.2) on a spatial staggered lattice, on which electrons and positrons respectively occupy odd and even sites, while the gauge field is implemented by Wilson lines (also called link fields) between adjacent sites. Further performing a Jordan-Wigner transformation which maps the fermions to spins [15], the lattice Hamiltonian takes the form where a is the spatial lattice spacing, n is the lattice site index, σ 1,2,3 are the Pauli matrices with σ ± (n) ≡ 1 2 (σ 1 (n) ± iσ 2 (n)), and α = F/g is the dimensionless background field. L(n) ≡ (E(n) − F )/g is the dimensionless lattice electric field (with the background field subtracted) operator defined on the n-th link between n-th and (n + 1)-th sites, and the link phase on the n-th link θ(n) is related to the gauge field by θ(n) ≡ agA 1 (n). 3 The total number of lattice sites N has to be even in this formulation.
We introduce dimensionless lattice observables by defining the dimensionless parameters in terms of which the dimensionless Hamiltonian is aH lat ≡ W/2x, where W given by The dimensionless lattice electric field L(n) obeys the lattice Gauss law As in the continuum case, we impose periodic boundary conditions (PBC) on the lattice, so the integer n in (2.4) and (2.6) runs from 0 to N − 1.
Note that l can be any integer, after we canonically quantize the gauge field sector-this is explained in Appendix A after (A.9). As such, the lattice Hilbert space is actually infinite dimensional, even with a finite N . We perform a truncation on the electric field, such that |l(n)| ≤ l max on all links n. This results in a finite dimensional lattice Hilbert space H lat . For the low energy physics of interest in this work, the error of lattice simulation is negligible as long as the cutoff l max is sufficiently large. In practice, how large the cutoff l max needs to be depends on the coupling m/g and the spatial size mL. We summarize in Appendix B the explicit values of l max used for various parameter choices in our lattice simulation. The finite dimensional lattice Hilbert space H lat is spanned by the states in (2.9), |{l}, {σ 3 } j , j = 1, · · · , dim H lat , with the truncation on l so that |l(n)| ≤ l max for all links n. For the j-th basis state |{l}, {σ 3 } j , we denote the eigenvalues of the electric field on the n-th link and the spin on the n-th site respectively by l(n) j and σ 3 (n) j . Then, we have matrix elements {l},

Observables and probability
Given a state |ψ expanded in the lattice basis, we compute the main lattice observables in this work.
• The expectation value of the spatially averaged electric fieldĒ is • The standard deviation of the electric field is • The probability of measuring the electric field value E m /g on the n-th link can be computed by first constructing the projection operator on the latticê where the sum j is over all the basis states with l(n) j = Em g − α. The probability is then given by the expectation value of the projection, (2.14) The quantum state on all odd (even) links is identical due to the periodic boundary conditions. There is a slight difference between odd and even links due to the staggered nature of the lattice, but this brings no important modifications to our analysis. Hence, without loss of generality we will focus on the probability of measuring the electric field on the central link (i.e. ( N 2 − 1)-th link) throughout this work.
• To calculate the number of pairs in a state, we recall that in the staggered fermion formalism, (−1) n σ 3 (n) = +1 when an electron (positron) appears at the n-th site where n is odd (even), and (−1) n σ 3 (n) = −1 when the n-th site is empty. The basis states |{l}, {σ 3 } j are eigenstates of electron (positron) number operatorN e (N p ): The number of pairs operator is defined asN pairs = 1 2 (N e +N p ) with eigenvalues N pairs,j = (N e,j + N p,j )/2. In the state |ψ as given by the expansion (2.10), the number of pairs is given by (2.17)

Exact Diagonalization and time evolution
We study both static properties and dynamics of the compact massive Schwinger model on the lattice, using the simple technique of exact diagonalization: after constructing the Hamiltonian matrix, we numerically diagonalize it to determine the energy levels and eigenstates. Time evolution is done by decomposing the initial state into Hamiltonian eigenstates and multiplying them by the phase factors exp{−iE i t} with E i being the i-th Hamiltonian eigenvalue.

Bounds on lattice parameters
In order to simulate the continuum theory on a lattice with spacing a, we need to characterize the range of lattice parameters that ensure the lattice results are physical in the continuum limit. The two bounds we impose on the lattice simulation are as follows.
• Lattice artifacts are suppressed if: • Semiclassical instanton calculations that we will introduce in the later sections are valid if: Combining these constraints, valid lattice simulation results are obtained within the range for x: In addition, an extrapolation procedure is required for the lattice results to approach the continuum limit. Unlike most previous work [16,17] which explored the properties of the noncompact Schwinger model, we are interested in the finite size effects brought by the compact spatial dimension. The size of the spatial circle L = N a is now a relevant quantity. Instead of using the two-step extrapolation, N → ∞ followed by x → ∞ as proposed in [16,17], we keep N √ x fixed and take N → ∞ to approach the continuum limit. This extrapolation approach keeps fixed the size of the space with respect to the coupling gL, while taking the spacing ga → 0. The continuum limit is thus approached with the finite size effects intact. For more details about the continuum extrapolation, see Appendix B and Appendix C.

Lattice dynamics and decay of an initial state electric field
In this section, we proceed to investigate the dynamics of the massive Schwinger model on a circle making use of the lattice simulation code. We focus on the time evolution of the electric field and measure the timescales for certain types of quantum quench. In addition, the energy splittings relevant to those typical quenches are identified. This provides a computationally cheaper way to study the dynamics of the electric field in the weak field regime.

Initial state and time evolution
To study the decay or time evolution of an initial state electric field we need to specify an initial quantum state that we will then evolve in time using the lattice Hamiltonian. Throughout this paper we will use the ground state of the α = 0 Hamiltonian as the initial state. Being the ground state, this should be a "minimum uncertainty state" in the sense that it minimizes some product of the variance in the field and its time derivative.
Due to symmetry this α = 0 ground state has Ē g = 0. We then evolve this initial state using the (numerically diagonalized) Hamiltonian with α = k/2 for some k ∈ R (this is a quantum quench). The spectrum of the Hamiltonian is periodic under α → α + 1 via spectral flow, but the initial state defined this way has Ē g = 0 + α with respect to this Hamiltonian. Hence the time evolution from this state reflects the behavior of an initial state electric field satisfying Ē = αg. 4 The observables calculated on the lattice include the expectation value and the standard deviation of the spatially averaged electric field, the probability to measure a dimensionless electric field value on the central link (C-L) (i.e. ( N 2 − 1)-th link) and the number of electron-positron pairs.
The dynamics is particularly interesting when the background fields are equal to halfinteger multiples of the coupling, α = k/2 with k ∈ Z. We consider the following three types of quenches: (1) |k| = 1 quench, (2) |k| = 2 quench, and (3) |k| = 3 quench. Without loss of generality, we will discuss quenches with k > 0 only; a quench with −k can be obtained immediately from the one with k by reversing the electric field direction. Figures 2 -4 respectively illustrate k = 1, 2, 3 quenches, the evolution behaviors of which vary with the spatial circle size mL.  The expectation values of the spatially averaged electric field Ē /g, the standard deviation of the spatially averaged electric field σ E/g , the probability to measure field E m on the central link P (E C-L = E m ), and the number of electron-positron pairs N pairs as functions of time in the k = 1 quench for m g = 50, with three different spatial circle sizes mL = 8 (top panel), mL = 4 (middle panel) and mL = 1 (bottom panel) on an N = 10 lattice. On the lattice, we set l max = 10 for mL = 8, 4 simulations and l max = 30 for mL = 1 simulation to achieve error O(10 −10 ) for the lowest four energy levels E/g.
By fitting the evolution data to the Ansatz  we extract the timescales τ /a, which are related to the periods T of the oscillations by T /a = 2πτ /a. As shown in Table 1, the timescales extracted respectively from the expectation values and measuring probabilities agree excellently.
It is striking that Ē /g oscillates roughly sinusoidally for all values of mL (left column, blue line of Figures 2, 3 and 4). However, there are significant qualitative and quantitative differences as mL varies which reveal the different origins of these oscillations. For mL = 8, the initial field is well-localized in the sense that its standard deviation satisfies σ E/g 1 and the period of oscillation is very long. From the probabilities for specific E-field values (middle column) and the time dependence of σ E/g (left column, orange line) one can see that the field is oscillating from Ē g = k/2 to  This behavior can be understood in terms of the novel instantons we discover, which have an action proportional to L (Section 5). When mL 1 the instanton action is large and the initial field is well-localized to a specific value, with exponentially suppressed probability to take any other value. As we will see, these instantons lack a negative mode. Rather than mediating a traditional decay, they predict an oscillation from E to −E, precisely as is seen in the top row of these figures. When mL is smaller the instanton action is not large and the field is not localized to a specific initial value (in units of g).
In this case the timescale is perturbative (the oscillation frequency for mL = 1 is within 10% of g/ √ π). One might wonder whether a different initial state with σ E/g 1 might be more pertinent to the question of electric field decay in this regime. This is not the caseany such state would include higher energy eigenstates and would very quickly evolve into a state with σ E/g 1. Indeed, this regime is in some ways similar to the massless or light regime of the Schwinger model in non-compact space, where the instanton action is small and the field decays by perturbative nucleation of light pairs of particles (except that in this small-circle regime the decay does not occur by pair production, as can be seen in the right column from the fact that the number of pairs is constant in time).
These behaviors can also be precisely understood in the bosonized theory (Section 4). The bosonized potential is symmetric when k is integer and has multiple local minima with a quadratic envelope. The behavior in the top row is similar to that of a symmetric double well potential when the barriers are high. The behavior in the lower row is that of a quadratic potential with small "wiggles" superimposed. The frequency and amplitude of these oscillations found on the lattice agree quantitatively with those predicted by the bosonized version of the theory.

k /
∈ Z quenches Let us turn to consider k / ∈ Z quench. As before the initial state is evolved using the final Hamiltonian with α = k/2, but now for k / ∈ Z. In Figures 5 and 6 we show the lattice simulation results of k / ∈ Z quenches in both large and small mL regimes.
In the large mL regime (cf. the top row of Figures 5 and 6), the electric field Ē /g undergoes small amplitude and relatively rapid oscillations with probability 90% to coincide with its initial value. This behavior is in sharp contrast to the case k ∈ Z. It is also instructive to compare to the behavior of the massive (m/g 1) Schwinger model in the non-compact case. There, initial fields with −1 < k < 1 (−0.5 <Ē/g < 0.5) are stable because the energy density in the field between a pair of nucleated particles would exceed that in the background field. By contrast, k = 3.4 would decay by Schwinger pair production. Instead, our numerical results show that in the small-circle regime this field is stable (as we will see, there is indeed no instanton that can mediate its decay).
For small mL (cf. the middle and bottom rows of Figures 5 and 6) the behaviors are qualitatively similar to k integer. Again, these behaviors are simple to explain both in terms of our novel instantons and the bosonized description. When mL 1 the would-be instanton action is large and the field is well-localized. However, for non-integer k the instanton ceases to exist and there is no decay channel for the field. In the bosonized description this corresponds to an asymmetric potential with local minima; an initial state localized in one minimum with high barriers has no symmetric state to oscillate to, or any other decay channel. The non-compact massless Schwinger model would behave in a qualitatively similar fashion.
In Figure 7 we plot the minimum probability over a long time period of measuring the central link electric field to be k/2 as a function of the background field at t = 0, i.e. F/(g/2) = Ē t=0 /(g/2) = k, which is again equal to the initial spatially averaged electric field. If this probability remains close to 1 it means the initial field is stable, while if it  is close to zero it means the field decays or oscillates from its initial value to some other values.
In the large mL regime, e.g. mL = 8 for m/g = 50 as shown by the orange points in Figure 7, when k is not (close to) an integer value the field is very likely to remain at its initial value for an extremely long time (e.g. (t/a) max = 10 8 in our lattice simulation). The sharp dips in P min (E C-L /g = k/2) appear at integer k, corresponding to the long-timescale quenched tunnelings shown in the top panels of Figures 2, 3 and 4. For k / ∈ Z, the timescale associated with the dominating oscillation mode is much shorter than the one for k ∈ Z. One can also see (insets) that the width of this resonance decreases exponentially with increasing k. This is expected because the instanton action scales linearly with k, and (as we will see) in the bosonized description the number of barriers separating E and −E also  In the small mL regime, e.g. mL = 4 for m/g = 50 as denoted in Figure 7 by the green points, the minimum probability for the field on the central link to be equal to the initial average electric field is noticeably away from unity even for the trivial quench with k = 0. As k increases, the central link field gradually loses the memory of the initial value of Ē /g in the sense that the probability P E C-L = Ē t=0 almost vanishes at some later times. In contrast to the large mL case, there are no dips in P min (E C-L /g = k/2) at integer k.

Relation to energy splittings
For large mL, the long-period sinusoidal oscillations in time suggest that only two Hamiltonian eigenstates with small energy splitting are involved, as is the case for tunneling between two symmetric wells in quantum mechanics (this is indeed the case in the bosonized description, as we will see in Section 4). Consider two states Ē /g = ±k/2 ≡ Ē ± in which the electric field is localized aroundĒ/g = ±k/2. Define |S and |A respectively to be the symmetric and antisymmetric superpositions of Ē ± . If |S and |A are Hamiltonian eigenstates, the transition amplitudes between the wells are given by where E S and E A are Hamiltonian eigenvalues for |S and |A respectively, withĒ ≡ (E S + E A )/2 and ∆E ≡ E A − E S . It follows that the transition probability from one well to the other at time t is In Table 1 we compare the oscillation timescale τ from a cosine fit to the lattice real-time data to the inverse of the energy difference ∆E between the k-th eigenstate and the (k − 1)th eigenstate. We find very precise agreement. (This numerical agreement is robust for lattices with varying numbers of sites N , in the large mL regime with 1 m/g < 50) τ a Ē /g and τ a P are extracted from lattice data of the electric field Ē /g and the quantum mechanical probability using fits (3.1) and (3.2) respectively. ∆E 1,0 , ∆E 2,1 and ∆E 3,2 are respectively the energy differences between the first excited and the ground states, between the second and the first excited states, and between the third and the second excited states of the Hamiltonian with α = k/2. The column of a∆E lists the lattice data of the relevant energy differences. The type of a quench is indicated by k.
The massive Schwinger model spectrum is periodic in α with period one, so one only needs to look at α = 0 or α = 1/2 spectrum for k ∈ Z quenches. The energy splittings ∆E/g converge as we perform the fixed N/ √ x continuum (N → ∞) extrapolation. This implies that the energy differences and the quenched tunneling effects are physical rather than lattice artifacts. We also verified numerically that the symmetric and anti-symmetric combinations of these Hamiltonian eigenstates Ē ± indeed correspond to states with expectation value for the spatially averaged field Ē /g = ±k/2 and standard deviation σ E/g 1. For small mL, the oscillations in time deviate from pure sinusoidal functions. It is instructive to perform a discrete Fourier analysis of the evolution data to extract the timescale of the quench and determine which energy differences dominantly govern the dynamics. Consider a sequence of time evolution data {q nt } where n t = 1, · · · N t with N t being the total number of time steps. In our case, q nt represents either Ē /g or the relevant transition probability at time t nt /a. The discrete Fourier transform from The power spectrum for the quantity q is then defined by S q (k) = |Q k | 2 . After a standard rescaling, we obtain the power spectrum for observable q as a function of dimensionless frequency, S q (f /a −1 ). The significant peaks in the power spectra of the electric field expectation value S Ē /g (f /a −1 ) and the relevant transition probability S P (f /a −1 ) correspond to the characteristic frequencies of the quench. Each peak (f /a −1 ) peak is related to a character- Comparing energy differences ∆E c /g with those in the Hamiltonian spectrum, we find that the ∆E c /g derived from the most predominant peak in the power spectra is ∆E 1,0 /g for all three types of quenches, in the small mL regime (e.g. the middle panels of Figures 2, 3 and 4). The energy differences of the higher excited states such as ∆E 2,1 and ∆E 3,2 play a minor role in these quenches, serving to modulate the amplitude of oscillations.

The bosonized description
The massive Schwinger model defined by the Lagrangian (2.1) admits an equivalent bosonized description [11,[18][19][20][21] -a theory of a massive real scalar field with a cosine interaction potential. The Euclidean action of the bosonized theory takes the form Here g and m are the gauge field coupling and the mass of the Dirac fermion in the original theory (2.1), and c is a dimensionless prefactor of the cosine potential. The background angle θ = 2πF/g [11]. On a compact spatial dimension with identification x ∼ x + L, the scalar field is periodic, φ(x + L) = φ(x), and admits a mode expansion If the size of the spatial circle is small enough and we are interested in the low energy physics, we can truncate the massive Kaluza-Klein (KK) modes with p = 0. The (1 + 1)dimensional field theory then reduces to a (0 + 1)-dimensional quantum mechanics of the zero mode. As we will see, this quantum mechanics generates accurate predictions for the energy differences relevant to the quenches discussed in Section 3. The dimensional reduction yields whereφ ≡ √ Lφ (0) and '· · · ' stands for terms involving massive KK modes. The effective potential in the quantum mechanical theory For future convenience, we introduce the dimensionless variable ϕ ≡φ/ √ πL and the dimensionless potential The dimensionless potential in terms of ϕ is explicitly given by where U 0 is chosen such that the global minima of U (ϕ) vanishes, i.e. U min = 0.

Numerical bosonized quantum mechanics
We first study the bosonized quantum mechanics by numerically solving the timeindependent Schrödinger equation In the above u(φ) ≡ φ ψ is the wavefunction of the state |ψ in theφ representation, and E is a Hamiltonian eigenvalue. We further write the wavefunction as u(φ) = Cφ ,ϕ u(ϕ) with u(ϕ) ≡ ϕ|ψ and Cφ ,ϕ a dimensionful constant determined by the normalization condition. 5 To make connections to the lattice simulation data, we rewrite (4.5) in terms of lattice parameters and the dimensionless wavefunction u(ϕ) as By solving (4.6), one obtains the dimensionless Hamiltonian eigenvalues and dimensionless eigenstates. From these we immediately know the energy difference ∆E i,j /g = (a∆E i,j ) √ x between the i-th and the j-th eigenstates. 6 Additionally, under bosonization the spatiallyaveraged electric field operator is proportional to the ϕ operator, It follows that the expectation value and the standard deviation of the electric field for the i-th Hamiltonian eigenstates |ψ i are With the Schrödinger equation and the observables to compute at hand, we can use the standard shooting method to solve the quantum mechanics numerically. To do so we need one more input, the value of the prefactor c. This factor hasn't been computed analytically for arbitrary coupling g. Instead, we determine the prefactor c for given m/g and mL numerically by comparing the predictions of a single benchmark quantity from the numerical quantum mechanics and the lattice simulation data. We choose this benchmark quantity to be the energy difference between the first excited and ground states ∆E 1,0 /g for α = 0.5. By numerically solving the quantum mechanics for a set of trial values of c, we find the optimal value of c up to O(0.01) that generates the closest value of the benchmark 5 The constant Cφ ,ϕ is determined such that dφ u(φ) 2 = dϕ |u(ϕ)| 2 = 1. 6 The ground state is referred to as 0-th state.
quantity to the one obtained from lattice data for specified m/g, mL and N . Let us denote the prefactor obtained in this way by c(m/g, mL; N ). Tables 2 and 3 summarize the optimal c(m/g, mL; N ) obtained from numerical quantum mechanics (nQM), the lowest few energy differences and σ E/g,0 and σ E/g,1 computed using the same optimal c, together with the same quantities calculated using lattice simulation (LAT) for the examples of m/g = 50, mL = 8, 4, 1 on N = 10 lattice.
Once the best-fit value c(m/g, mL; N ) is determined, we compute other observables with this c and find that the bosonized theory indeed corresponds to the original fermionic theory. For the examples presented in the tables, the energy differences between higher excited states and σ E/g,0 , σ E/g,1 are close to their lattice counterparts even for α = 0.5. We confirm that the consistency holds true throughout the ranges of m/g and mL we explored. In this sense, the prefactors c(m/g, mL; N ) are quite reliable. The physical prefactor, which is independent of lattice parameters, is attained by a continuum extrapolation of c(m/g, mL; N ). More details about the extrapolation procedure and the prefactor's dependence on m/g and mL are summarized in Appendix C.
mL α optimal c optimal c optimal c ∆E 1,0 /g  Table 2: The lowest three energy gaps and the prefactor c in the bosonized action (4.1) for m g = 50 with various mL and diffrent background field α on N = 10 lattice. The values in the '(nQM)' columns are calculated using numerical bosonized quantum mechanics while those in the '(LAT)' columns are obtained from lattice field simulation. The 'optimal c (nQM)' is the prefactor c(m/g, mL; N ) elaborated in Subsection 4.1. The 'optimal c (Garg (4.14))' and the 'optimal c (WCL (4.19))' are the prefactors c that produce the lattice simulated value of ∆E 1,0 /g according to (4.14) and (4.19), respectively, at α = 0.5. The '×' sign indicates that the formula is invalid for certain mL as expected, and thus the value of c there is listed for comparison purpose only. For the same mL but different α, we use the same optimal c to carry out the numerical quantum mechanics calculation.  Table 3: The standard deviations σ E/g for m g = 50 with various mL and diffrent background field α on N = 10 lattice. The values in the '(nQM)' columns are calculated using numerical bosonized quantum mechanics while those in the '(LAT)' columns are obtained from lattice field simulation. Figure 8 shows the effective quantum mechanical potential and the wavefunctions u(ϕ) of relevant states for the quenches with k = 1, 2, 3 in large mL regime (cf. the top panels of Figures 2, 3 and 4). For the k ∈ Z quench, the relevant states are |ψ k and |ψ k−1 . Due to sufficiently high barriers between the wells, the relevant wavefunctions are locally supported in the two wells of the potential, i.e. around ϕ = ±k/2. These states exactly fit in the role of the symmetric and the antisymmetric states proposed in Section 3.2, confirming the double-well tunneling scenario mentioned there. The corresponding parameter regime of m/g and mL is referred to as the double-well regime. Figure 9 illustrates the potential, and the lowest four wavefunctions u(ϕ) for the k = 1, 2, 3 quenches in small mL regime (cf. the middle panels of Figures 2, 3 and 4). In this intermediate regime the potential barriers are still significant but are no longer high enough to confine the probability to a single well or symmetric pair of wells, and the wavefunctions spread across multiple wells. We refer to this parameter regime as the multiwell regime.
When mL is even smaller, the cosine term in the potential V (φ) becomes negligible in comparison to the quadratic term and potential approximately reduces to that of a simple harmonic oscillator. This scenario is illustrated in Figure 10, and we refer to it as the nearly SHO regime.

Semi-analytical bosonized quantum mechanics
In this subsection, we present semi-analytical calculations of the ground splitting ∆E 1,0 /g in the double-well regime, and of the standard deviations σ E/g,0 and σ E/g,1 in all three regimes for a k = 1 quench.
For a k = 1 quench in the double-well regime, the tunneling is occurring between the lowest two neighboring wells of V (φ). The states relevant to the tunneling are the ground and the first excited states. When it comes to these two states, V (φ) approximately reduces to a symmetric double-well potential. For the symmetric double well potential problem,  one can solve it using either WKB [22] or instanton method [23]. Here we take the instanton approach due to [24].
Let us denote the locations of the two global minima of V (φ) byφ ± . The action of the single instanton involving the quantum tunneling betweenφ − andφ + is where ϕ ± =φ ± / √ πL. Expanding U (ϕ) [i.e. eq.(4.4)] with α = 1/2 around ϕ + gives the frequency of the small oscillation around this well: The ground splitting is then given by [24] ∆E 1,0 = 2ωφ ωφφ 2 + π 1 2 e A e −S inst , (4.12) where ωφ = g √ π ω ϕ and The ground energy splitting in the unit of coupling is then (4.14) In the weak coupling regime m/g 1, ϕ + ≈ 1/2 and then The ground splitting in the weak coupling limit (WCL) reduces to We determine the prefactor c by comparing (4.14) and (4.19) with lattice data. The values of c are summarized in Table 2 for the examples of m/g = 50, mL = 8, 4 on N = 10 lattice. They are consistent with those derived from numerical bosonized quantum mechanics in the double-well regime. For m/g = 50 and mL = 4 in the multiwell regime, the prefactor c derived from (4.14) and (4.19), as distinct from the numerical quantum mechanical calculation, is not valid anymore.
Another set of useful observables are the standard deviations in electric field for the lowest eigenstates. As we will see in Subsection 4.3, they play the role of order parameters characterizing double-well-to-multiwell and multiwell-to-SHO phase transitions. Here we compute the standard deviation for the ground state σ E/g,0 and for the first excited state σ E/g,1 .
In the double well regime for k = 1, the wavefunctions of the ground state u 0 (ϕ) = ϕ|ψ 0 and the first excited states u 1 (ϕ) = ϕ|ψ 1 are locally supported around ϕ = ± 1 2 , symmetric and antisymmetric about ϕ = 0 respectively. Therefore for both states we have Ê /g ≈ 0 and Ê2 /g 2 ≈ 1/4. It then follows from (4.9) that σ E/g,0 ≈ 1/2 and σ E/g,1 ≈ 1/2. Because each of the wells is not perfectly symmetric about the local minimum, and the wavefunctions u(ϕ) have small but non-vanishing supports in the wells beyond the central double wells, these values are not exact.
In the multiwell regime, as the wavefunctions of the ground and first excited states are nonvanishing in higher wells |ϕ| > 1 2 , σ E/g,0 and σ E/g,1 exceeds 1/2 significantly. In the nearly SHO regime, the quadratic term dominates the effective potential V (φ). In the crudest approximation, one may neglect the cosine term and treat the potential as a pure quadratic one. The dimensionless Schrödinger equation (4.6) reduces to This is nothing but the Schrödinger equation for a simple harmonic oscillator of unit mass with angular frequency With facts about a quantum SHO and the operator relation (4.7), the standard deviations σ E/g,0 and σ E/g,1 are readily obtained as In the above |0 and |1 are the ground and the first excited eigenstates associated with the lowest two E ϕ eigenvalues, respectively. They correspond to the physical Hamiltonian eigenstates |ψ 0 and |ψ 1 for a fixed gL. In addition, from the SHO's equally spaced energy levels ∆E ϕ = Ω ϕ , one deduces that ∆E i+1,i /g = 1/ √ π ≈ 0.56419. This is exactly the energy gap ∆E 1,0 /g in the massless (m → 0) Schwinger model where the bosonized potential V (φ) is exactly quadratic. As we argued, in the massive Schwinger model on a circle of sufficiently small size mL (i.e. in the nearly SHO regime), the lowest energy difference ∆E 1,0 /g is approximately 1/ √ π.

Three qualitative regimes
In the preceding subsections, we recognized three regimes that we termed the doublewell (DW), multiwell (MW) and nearly SHO regimes, according to the behavior of electric field evolution and to the configurations of potential and lowest-lying wavefunctions. Qualitatively the distinction is that in the DW regime the barriers are high enough that the electric field can be localized to a specific value and will remain at this value for a long time with high probability, whereas in the MW and nearly SHO regimes it will rapidly evolve to other values.
For a fixed m/g, the three regimes can be reached by tuning the size of the spatial circle mL. The right columns of Figures 8, 9 and 10 illustrate the transitions from DW to MW and from MW to nearly SHO in terms of wavefunctions of the ground and the first excited states. We now proceed to determine the boundaries separating the regimes in terms of mL by using the electric field deviation σ E/g as an indicator.
As discussed in subsection 4.2, σ E/g ≈ 1/2 for both the ground and the first excited state in DW regime while it significantly exceeds 1/2 in MW regime for α = 1/2 (i.e. k = 1). For a given m/g with α = 1/2, we define the DW-MW boundary x bdy on the N -site lattice as the value of x such that where i = 0, 1 corresponds to the ground and the first excited states respectively. Figure  11 plots σ E/g,0 and σ E/g,1 as functions of x for m/g = 50 on N = 10 lattice. Figure 12 shows the DW-MW boundaries x bdy based on both σ E/g,0 and σ E/g,1 for m/g = 50 and N = 6, 8, 10, 12, 14, 16. The data points (x bdy , N ) are fit pretty accurately to with A a fitting constant. Together with N √ x = mL m/g , which is fixed in our continuum extrapolation, we identify the factor A = mL bdy m/g . In this way, we obtained the physical DW-MW boundary mL bdy . The fitting curves and the corresponding mL bdy are plotted and compared against various mL in Figure 12.
The transition between MW and the nearly SHO regimes can be seen in Figure 11 as well. As x gets very large and thus the spatial circle size L shrinks to about a single fermion Compton wavelength, both σ E/g,0 and σ E/g,1 approaches the values for a SHO (cf. (4.21) and (4.22)).

Worldline instantons in the massive Schwinger model
In this section, we use the worldline formalism of path integrals to study tunneling processes in the massive Schwinger model on a circle. We identify the worldline instantons responsible for k ∈ Z quenched tunneling without pair production when 2m 2 gE > mL 1, i.e. the DW regime (e.g. top rows of Figures 2, 3 and 4), and compute the corresponding one-loop transition amplitudes.
We begin with a review on the worldline formalism applied to quantum electrodynamics, following [10,25,26]; readers who are familiar with it may skip directly to the next subsection. Consider quantum electrodynamics in d-dimensional Euclidean spacetime with  Figure 11: The standard deviations in the spatially averaged electric field for the ground state σ E/g,0 and for the first excited state σ E/g,1 as functions of x ≡ 1/(ga) 2 on N = 10 lattice with m/g = 50 and α = 0.5. A smaller spatial circle size corresponds to a larger x according to mL = m g N √ x . The blue and orange points are respectively the lattice data of σ E/g,0 and σ E/g,1 . The green dot-dashed curve denotes σ E/g,0 for a pure quantum SHO given by (4.21) whereas the red dotted curve is σ E/g,1 for the SHO given by (4.22). The purple dashed line denotes the constant σ E/g,0 = σ E/g,1 = 1/2, which corresponds to the perfect DW case.
the gauge field A. The Euclidean action is obtained from the Lorentzian action by a Wick rotation, in which the Euclidean time x d is defined by t = −ix d with t being the Lorentzian time. Accordingly the Euclidean gauge field A µ is related to its Lorentzian counterpart A L µ as A L 0 = iA d , A L j = A j , (j = 1, · · · , d − 1). Let us first consider scalar electrodynamics, with A coupled to a complex scalar φ, and look at the normalized vacuum survival probability amplitude expressed as a path integral: where µ, ν = 1, · · · , d and F µν = ∂ µ A ν − ∂ ν A µ is the Euclidean field strength. Here an appropriate gauge-fixing is to be performed, and we have included the boundary term in S A [A] in order to have a well-defined variational problem with F µν fixed at infinity as boundary condition [27,28]. Integrating out φ, we obtain, at fixed A, the Euclidean effective action Γ E,scalar [A], given by The boundaries x bdy derived from σ E/g,0 and σ E/g,1 on the lattices of N = 6, 8, 10, 12, 14, 16 sites are shown as blue and orange dots respectively. Blue and orange solid curves are the fitting curves (4.24) of σ E/g,0 and σ E/g,1 respectively. The boundary in terms of the physical space size, mL bdy , is derived from the fitting parameter A = mL bdy m/g . The white region above the fitting curves is the double well regime while the orange/brown region below the fitting curves is the multiwell regime. The gray shaded rectangular region (x < (m/g) 2 ) is ruled out due to the violation of eq.(2.21). The gray shaded region in the bottom (x > (N · m/g) 2 , i.e. mL < 1), which is beyond the semiclassical regime, requires much larger l max to explore and is not presented in this work.
whereẋ µ ≡ dx µ /du, This expresses the effective action Γ E,scalar [A], at fixed A, as a path integral of a closed worldline x µ (u) with periodic boundary condition, integrated over the proper time s.

For spinor quantum electrodynamics, we introduce Feynman's spin factor Φ[x(u), A]:
where Tr γ is the trace over the (Euclidean) Gamma matrices, P denotes path-ordering, ψ µ ≡ dψ µ /du, and in the second line, we have expressed it as a fermionic path integral using the coherent state formalism [29][30][31], which is necessary for a consistent semi-classical analysis [32,33]. The effective action is obtained by augmenting (5.8) with Φ[x(u), A], Including the dynamics of the gauge field, the vacuum survival amplitude (5.1) (for both scalar and spinor electrodynamics) becomes The last line requires some explanation. In general, the exponential and the expectation value in the second equality do not commute. And since, from (5.8), Γ E [A] has a natural interpretation as a sum contributed by worldline single-instanton solutions, the non-commutativity is the result of correlations between single-instantons at different locations in the Euclidean spacetime, due to the dynamical gauge field [10]. In the weak-field limit, we can make the dilute instanton gas approximation, wherein these correlations are ignored and one arrives at (5.16). Explicitly, in this approximation, (5.16) becomes [34] At each n, there are n uncorrelated single-instantons at different locations in the Euclidean spacetime, with the factor 1/n! accounting for the exchange of identical instantons. The integral over A is evaluated semi-classically by expanding around a background classical solution A cl,n , and integrating over the fluctuation δA.
From now on, we focus on the spinor case. In this case, it is more convenient to evaluate the effective action in (5.17) by first performing the path integral Z(s) in (5.13) for each worldline, following [26]. 7 Factoring out x(0) = x(s) =x from the path integral, Z spinor (s) reads The classical equations of motion satisfied by the worldline instantons (x µ (u), ψ µ (u)) are which are to be supplemented by that of the gauge field. One obvious class of solutions is with ψ cl = 0, so that at the classical level, the problem reduces to that of scalar electrodynamics. We assume that there are no non-trivial classical fermion solutions contributing to the path integral. In (5.18), we expand the worldline x(u) as a sum of the classical path x cl (u) and fluctuations δx(u), Here the fluctuation δx µ (u) satisfies the Dirichlet boundary condition, since we have factored out the endpoint x(0) = x(s) =x from the path integral in (5.18). This fluctuation gives the following contribution to the path integral: The subscript D in the determinants refers to Dirichlet boundary condition. The factor 1/(4πs) d 2 is a normalization factor; when gF µν (x cl (u)) = 0, the path integral reduces to a free particle path integral with mass m = 1/2. The fermionic path integral around the trivial solution gives the following contribution [35] det where the subscript A refers to anti-periodic boundary condition. We have set the normalization of this ratio of determinants to be unity to be consistent with the fact that, when F = 0 on the worldline, the spin factor defined in (5.10) evaluates to Φ = 1. Including the on-shell contribution exp(−S[x cl ](s)),Z(s,x) in (5.18) becomes (5.23) 7 In Appendix D, we provide alternative calculations for the worldline instantons studied in this section, but for scalar quantum electrodynamics; there, we first integrate over the proper time, and then perform the path integral, as in [10].

Tunneling on a circle
On a small circle of circumference L < 2m/gE, pair production is exponentially suppressed, and a transition of the electric field without pair production dominates. We study the problem with the dynamics of the gauge field turned on, taking into account the backreaction of the worldline instantons to the gauge field. We are interested in finding the worldline instanton solutions that mediate the transition amplitude from one value of the electric field to another.
As in (5.17), we express an (un-normalized) transition amplitude as a sum of multiinstantons, which involves a path integral over the dynamical gauge field, (5.26) and E is the Lorentzian electric field which is related to the Euclidean field strength F µν as E = iF 12 = −iF 21 . In (5.25), we used the dilute instanton gas approximation to commute the gauge field path integral with the sum over multi-instantons. Each Γ (n)

E [A] is contributed by n worldline instantons, satisfying the boundary conditions
For each n, we have n charged point particles x i (u i ) coupled to the gauge field by the Wilson line action along their worldlines C i . The terms in the action involving A read Varying the action wrt. A gives Gauss' law; focusing on the trivial fermion solution ψ cl = 0, which simply says that the on-shell electric field F cl changes by one unit of g as it crosses a charged worldline. Substituting (5.28) back to the action (5.27), the on-shell action is After that, one integrates over the fluctuation δA =δA+δA in (5.25). We have factored out the variationδA due to the variations of the worldlines and the fermions; the corresponding contributions are accounted for in Γ (n) . As it turns out, the quadratic action of the remaining fluctuationδA is that of a free field, decoupled from the worldlines and their fluctuations, as long as the worldlines do not intersect. Therefore, we can absorb its contributions in the normalization of the amplitude.

|k| = 1: Straight line instantons
By calculating the transition amplitude E f = −g/2|E i = +g/2 , we will argue that it is the straight line instantons that account for |k| = 1 quenched tunneling without pair production when mL 1 (cf. the top row of Figure 2). We first consider the survival probability amplitude at E = +g/2, +g/2| + g/2 . The gauge field satisfies the boundary condition E(t E → ±∞) = +g/2. We wish to express the amplitude as a sum of multi-instantons, E is contributed by n worldline instantons.
x 1  The n-worldline instanton solution which leads to (−Γ (n) E ) and satisfies the boundary condition E(x 2 → ±∞) = +g/2 only exists for n even. It is given by n straight lines running in the compact x 1 -direction with alternating orientations, with n/2 of them in the (+x 1 ) direction and n/2 in the (−x 1 ) direction. This is illustrated in Figure 13a. From (5.12) and (5.18), we have We take the electric fieldĒ ≡ E(x icl (u)) on the worldline to be the average (mean) of that at its two sides. On the trivial fermion solution ψ cl = 0, from Gauss' law (5.28) and E(x 2 → ±∞) = +g/2, we find thatĒ ≡ E(x icl (u i )) = 0. The worldlines as solutions to (5.19) are given by Since E(x icl (u i )) = 0, the ratio of determinants due to the fluctuations (for both bosonic and fermionic) is exactly one, so we havẽ Hence, (5.31) becomes where K ν (x) is the modified Bessel function of the second kind, mL 1 is the semi-classical regime. Moreover, from (5.29), the on-shell action of the gauge field is which can be absorbed in the normalization of the transition amplitude. Therefore, the survival amplitude (5.30) is (5.40) From this, it is apparent that we have normalized the amplitude correctly as we threw away the infinite gauge field action. The exponent of the time scale associated to this process is independent of the electric field value E = ±g/2. Thus, this process cannot be mediated by the production of a real pair by taking energy out of the electric field. We can similarly compute −g/2| + g/2 , in which case the multi-instanton solutions are n straight lines, n odd, along x 1 , with (n − 1)/2 running in the (+x 1 ) direction and (n + 1)/2 running in the (−x 1 ) direction. This is illustrated in Figure 13b. In the end, after an appropriate normalization, one gets −g/2| + g/2 = sin (t(m/π)K 1 (mL)) mL 1 The fact that +g/2|+g/2 and −g/2|+g/2 are given by cosine and sine and their square norms sum to one is in agreement with the bosonized description studied in Section 4, in which E = ±g/2 are the two degenerate global minima of the sine-Gordon potential, and the transitions between them are mediated by the Coleman double-well instantons yielding cosine and sine. Comparing (5.41) with (3.3), we extract the ground state energy splitting ∆E 1,0 /g for background E = g/2 from the straight line instantons, In Figure 14, we plot the instanton estimates (5.42) of E 1,0 /g for α = 0.5 in curves and compare them with the continuum extrapolated lattice results. For relatively large m/g and relatively large mL where the instanton semiclassical computation applies, the instanton estimates agree remarkably well with the lattice results.

|k| = 2: Lemon instantons
Next, we wish to identify the worldline instantons that compute the amplitude +g| − g . This corresponds to |k| = 2 quenched tunneling when mL 1 (cf. the top row of Figure 3). From (5.14) and (5.25), we have the Euclidean action Note that we have now taken u i ∈ [−s i /2, s i /2). For the trivial fermion background ψ cl = 0, the equations of motion (5.19) for each worldline (with the worldline index i suppressed) become ẍ 1 cl = +2giF 12ẋ 2 cl As before, we take the electric field E(x cl (u)) on the worldline to be the average (mean) of that at its two sides; Gauss' law (5.28) and E(x 2 → ±∞) = ±g imply that E(x i (u i )) = ±g/2 ≡ ±Ē.
x 1 We identify a single-instanton solution, illustrated in Figure 15, which we dub the "lemon instanton". It consists of two circle segments (1) and (2) running in the positive x 1direction, patched together at u = 0 and u = ±s/2, where the worldline has a continuous zeroth and first derivatives, while the second derivative changes sign. This solution is periodic since x 1 ∼ x 1 + L. The region enclosed by the worldline has E = 0, while outside we have E = ±g, and on the segments, E(1) = −E(2) =Ē ≡ g/2. Explicitly, the solution is given by  The on-shell action of (5.43) for this single-instanton solution takes the form

One-loop calculations
Around the lemon instanton, the bosonic fluctuation operator in (5.23) takes the form Here we have defined u ≡ sv, so that v ∈ [−1/2, 1/2). The eigenmodes of the fluctuation operator must have continuous zeroth and first derivatives. By diagonalizing the operator, the eigenmodes that satisfy Dirichlet boundary condition are easily found to be where U is the same matrix as in (5.52). Since the operator is first-order, the eigenmodes only need to be continuous at v = 0 (in addition to being anti-periodic at v = ±1/2). They are easily found to be where C 1,2 are constants. They have eigenvalues λ ψ,n = iπ(2n + 1)/s, where n ∈ Z. Thus, the eigenvalues are independent of gĒ, and so the determinant is canceled by that of the free operator. Substituting (5.49), (5.56) and (5.23) back to the effective action (5.12), we have Including the pre-factors, we finally get The dominant saddle, p = 0, has the same action as that of Brown's instanton in [2]. It is instructive to consider different limits. When L → (2m/gĒ) where pair production becomes kinematically favored, S p=0 → πm 2 /gĒ, which is the same as the on-shell action of Schwinger pair production. On the other hand, in the small circle limit L (2m/gĒ), the p > 0 saddles are exponentially suppressed. The p = 0 saddle dominates, with action approaching S p=0 → 2mL. The asymptotic form of the action can be understood intuitively: in this limit, the area enclosed by the lemon vanishes, and the action is given by (m times) the length of the worldline. As in the straight line case, the fact the on-shell action does not depend on the electric field value in the small circle limit suggests that the process does not involve the production of a real pair.
Moreover, in the present case, the gauge field fluctuationδA not due to the variation of the worldline again decouples; its contribution can again be absorbed in the normalization of the amplitude. In the end, if the lemon instanton were the only solution contributing to the amplitude, then we would get, in the weak-field and small circle limits, Comparing (5.67) with (3.3), the energy difference between the second and the first excited states ∆E 2,1 /g for background E = g can be obtained from the lemon instantons as However, by comparing the lemon instanton estimates (5.68) with the continuum extrapolated lattice results of ∆E 2,1 g for various m/g and mL, we find that they do not agree.
Instead, if we modify the prefactor by multiplying it with an extra factor Rc 2L = (m/g) 2 mL , where R c ≡ m gĒ = 2m g 2 is the critical radius of the lemon instanton, the modified lemon instanton estimate is then consistent with the lattice result for relatively large m/g and relatively large mL where the semiclassical calculation applies. In Figure 16, the modified lemon instanton estimates (5.69) of E 2,1 /g for α = 0 (the spectrum of which is equivalent to that for α = 1) are plotted in curves and compared with the continuum extrapolated lattice results for various m/g and mL.

Other instantons?
The discrepancy in the pre-factor of the tunneling time scale between the numerical results and that given by the lemon instanton (5.45) alone suggests that, among other possibilities, we may have missed some other contributing single-instanton solutions. One class of candidates is the "chain instantons", illustrated in Figure 17. Each of them is a worldline intersecting itself some odd number N lemon times, generalizing the single-lemon.  Figure 16: The continuum extrapolated (N/ √ x fixed, N → ∞) lattice results (dots) versus the modified lemon instanton estimate (5.69) (curves) of the energy difference ∆E 2,1 /g between the second and first excited states for background field α = 0. Top panel: each color represents a fixed value of m/g; bottom panel: each color represents a fixed value of mL.
One can apply the one-loop calculations of the single-lemon to these "N lemon -chains". In the small circle limit, the area enclosed by each N lemon -chain vanishes, and the action of each chain is given by S chain → 2mL, same as for the single-lemon. The one-loop prefactor is found to be N lemon m/L, up to a numerical factor. Therefore, after including these chains, the effective action would still take the form Γ E ∝ it m/Le −2mL , not to mention that the sum would be divergent. In other words, this still does not agree with our numerical results. We will investigate it further in a future work.

|k| ≥ 3: Higher-order instantons
We have only discussed instantons mediating the transition amplitudes between eigenstates with E = ±g/2 and E = ±g respectively. We propose that, for tunneling between higher electric field values (i.e. |k| ≥ 3 quenched tunneling when mL 1), the instantons can be constructed by combinations of the lemon and the straight line, in a way that obeys Gauss' law. For example, for +3g/2| − 3g/2 (cf. the top row of Figure 4), we expect the "lime wedge" instanton, illustrated in Figure 18, to contribute. In the small circle limit, its action is approximately 3mL.
Following this pattern, we predict that, in the small circle limit, the time scale of the transition amplitude +kg/2| − kg/2 for any positive integer k has an exponential dependence given by exp(kmL). This exponential dependence on k appears to be consistent with the energy gaps computed via the lattice and numerical quantum mechanics at least up to k = 3 (cf. Table 2). When the circle is sufficiently small that Kaluza-Klein modes can be neglected, it also follows from the bosonized quantum mechanical description in that there are k barriers separating the local minima corresponding to E = ±kg/2 (cf. Figure 8).

Summary and outlook
We have analyzed the massive Schwinger model with periodic boundary conditions using four distinct techniques, all of which agree quantitatively. Our analysis reveals several novel features. These results raise a number of questions for future investigations: • Finding the correct method for computing the pre-factor for the novel instantons described in Section 5 in the case k > 1 and including the chain instantons.
• Exploring the relation between non-perturbative effects in the fermionic and bosonized descriptions. Specifically, how are the (ferminonic) instantons described in Section 5 related to quantum mechanical tunneling in the bosonized description?
• The quantization of the electric field in the Schwinger model serves as a toy model for higher-dimensional theories with higher-form fluxes. Such theories can be used to build models of vacuum energy and inflation (for instance [36][37][38]). It would be interesting to investigate the implications of our findings for these theories, and to explore the bosonized version of this quantization.

A Lattice formulation of massive Schwinger model
In this appendix, we review the Kogut-Susskind approach [12][13][14] to the lattice Hamiltonian formulation of the massive Schwinger model. The lattice Hamiltonian (2.4), which plays a central role in our lattice simulation, will be derived. We also present here some facts about the lattice simulation that are omitted in the main text.

A.1 Staggered fermions
Naively putting spinors on a lattice results in the fermion doubling problem [39]. In 1+1 dimensions, this problem is completely resolved by the staggered fermion approach due to Kogut and Susskind [12], wherein the spinor indices are identified with spatial "indices". That is, the two components of a local Dirac spinor are placed on adjacent spatial sites. We put the upper components on even n sites while the lower components on odd n sites, i.e. We introduce dimensionless spinor φ = a 1 2 χ for later use; it is not to be confused with the dual bosonic field in the main text.

A.2 Discretization of the Hamiltonian
Our goal is to obtain the discretized version of Hamiltonian (2.2) in the temporal gauge A 0 = 0.

Fermion kinetic term
The massless Dirac term in terms of staggered fermions reads Here the gamma matrices in Dirac basis are given by with γ 0 = γ 0 and γ 1 = −γ 1 . A staggered fermion field is well-defined locally on a single lattice site whereas its spatial derivative contains information about the neighborhood of the site it sits on. Since the fermion kinetic term H D contains an "interaction" between the local field and its derivative ψ † ∂ 1 ψ, we need to specify where this "interaction" happens. It is convenient to make it happen on the same site ψ † sits on. For instance, ψ † e ∂ 1 ψ o (n) ≡ ψ † e (n)∂ 1 ψ o (n) is defined for even n, with the spatial derivative of ψ o on even n lattice given by 8 For odd n, ψ † o ∂ 1 ψ e (n) ≡ ψ † o (n)∂ 1 ψ e (n) can be defined in the same manner. Diagrammatically illustrated in Figure 19 is how spatial derivatives on staggered lattice are defined. The fermion field kinetic term (A.3) can then be written as "hopping" terms, 2a χ † (n) (χ(n + 1) − χ(n − 1)) + χ † (n + 1) (χ(n + 2) − χ(n)) if n even i 2a χ † (n + 1) (χ(n + 2) − χ(n)) + χ † (n) (χ(n + 1) − χ(n − 1)) if n odd (A.6) 8 We stress that for even n, ψe(n), ψo(n − 1) and ψo(n + 1) are well-defined according to (A.2). In contrast, ψo(n), ψe(n − 1) and ψe(n + 1) are meaningless for even n. On the lattice, the integral is translated into a summation of operators on all sites as 9 In the second line we redefine the dummy indices n ≡ n − 1 and n ≡ n + 1. In the last line the fermion field χ is transferred to the one-component dimensionless fermion field φ.

Fermion mass term
The discretization of the fermion mass term is The summation runs over either all even n sites or all odd n sites, depending on which expression in (A.6) is used. It is easy to see that these two summations are equivalent by manipulating the dummy index.

Gauge field kinetic term
In temporal gauge A 0 = 0, the field equation for gauge field is −Ȧ 1 = E, which indicates that the electric field E is canonically conjugate to A 1 (n) ≡ A(n). On the lattice this commutation relation is [A(n), E(m)] = i 1 a δ n,m . In terms of the dimensionless link phase θ(n) ≡ agA 1 (n) and the dimensionless lattice electric field L(n), the commutation relation reduces to [θ(n), L(m)] = iδ n,m . The most general operator L(n) takes the form of L(n) = F is a constant background field. The canonical commutation relation immediately implies that e ±iθ(n) is a L-shift operator on the n-th link, where |l(n) is the eigenstate of L(n) with eigenvalue l(n). As seen from (4.1), physics is periodic in θ ≡ 2πα with period 2π. Without loss of generality, we take α ∈ [0, 1) such that L(n) is an integer on the n-th link for a certain n. By Gauss' law (2.7), L(n) for all n are automatically integers. The gauge kinetic term is then discretized in terms of L(n) and constant α as

Gauge invariant interaction term
The gauge field defined on a lattice lives on the links connecting adjacent sites and has certain direction on each link. The gauge field "starting" from n-th site and pointing to (n + 1)-th site is defined as A 1 (n) ≡ A(n) ≡ A(n, +) ≡ A(n, n + 1) . (A.11) A consistency condition for the gauge field living between n-th and (n + 1)-th sites is A(n, n + 1) ≡ A(n, +) = −A(n + 1, −) ≡ −A(n + 1, n). Figure 20 illustrates the definition of gauge field between sites. Figure 20: Gauge field A(n, +) ≡ A(n) ≡ A 1 (n) living between sites is specified by its starting site and its spatial direction.
One can further define the Wilson line (or the link field ) from n-th to (n + 1)-th sites on one dimensional spatial lattice is related to the gauge field A 1 as U (n, +) ≡ U (n, n + 1) ≡ e iagA(n,n+1) = e iagA(n,+) = e iagA(n) ≡ e iθ(n) , To incorporate the gauge invariant interaction term H int = dx gψγ 1 A 1 ψ into the full Hamiltonian, one can simply promote the ordinary derivative to the covariant derivative. This turns the fermion kinetic term H D to the covariant kinetic term H D + H int . The covariant derivative of ψ o on even n sites is defined by which is generalized from (A.5). The covariant derivative acting on ψ e is obtained accordingly. Substitution ∂ µ → D µ promotes (A.3) to The covariant kinetic term becomes dressed "hopping" terms +χ † (n + 1)e iθ(n+1) χ(n + 2) − χ † (n + 1)e −iθ(n) χ(n) , (A.17) • for odd n, The discretization of the covariant term is derived as

Full lattice Hamiltonian
Putting everything together, the full Schwinger Hamiltonian (2.2) is discretized as (A.20)

A.3 Jordan-Wigner transformation
In 1+1 dimensions, fermionic degrees of freedom can be mapped onto spin degrees of freedom via Jordan-Wigner transformation which is explicitly given by which is eq.(2.4) in the main text.
B The cutoff l max and the continuum extrapolation of energy levels In our lattice simulation implementation, the dimension of the lattice Hilbert space is determined by the number of lattice sites N and the electric field cutoff l max . An excessive l max can be computationally expensive, while one that is too small could lead to physical effects being missed out from simulation results. Throughout this work, we use "sufficiently large" values of l max that achieve error O(10 −10 ) for the lowest four energy levels E/g. More explicitly, our choices of l max for the entire parameter regime of 1 ≤ m/g ≤ 50 and 1 ≤ mL ≤ 8 are 15, if 10 < m/g ≤ 30, 2 ≤ mL < 4 or 4 < m/g ≤ 10, 1 ≤ mL < 2 10, otherwise (B.1) As briefly mentioned in Section 2, we perform the finite volume continuum extrapolation with N √ x fixed and N → ∞ to extract physical quantities which are independent of lattice parameters. Among the important quantities is the energy gap ∆E i,i−1 /g. We fit the lattice results of the energy difference ∆E i,i−1 /g(m/g, mL; N ) to using the functions FindFit and NonlinearModelFit in Mathematica. The fitting parameters are β and ξ. And ∆E i,i−1 /g is the extrapolated value of the energy difference between the i-th and (i − 1)-th energy eigenstates.

C More about the prefactor c
In this appendix, we present more details about the prefactor c shown up in the bosonized action (4.1). We determine the physical values of c via the continuum extrapolation of the lattice data of c. Its dependence on the mass-coupling ratio m/g and the size of the spatial circle mL is then investigated.
As proposed in Section 4, we find the optimal c up to O(0.01) such that the numerical bosonized quantum mechanics generates the closest value of the ground energy difference ∆E 1,0 /g for α = 0.5 produced by lattice simulation for given m/g, mL and N . The prefactor c determined this way are denoted by c(m/g, mL; N ). In Tables 4 -7, we summarize values of c(m/g, mL; N ) for 1 ≤ m/g ≤ 50, 1 ≤ mL ≤ 8 and various lattice sites N .
In order to get the physical values of c which do not depend on N , we perform the finite volume continuum extrapolation by fitting c(m/g, mL; N ) to  Tables 4 -7.
To study the dependence of c on the mass-coupling ratio m/g and the size of the spatial circle mL, we further fit c [N →∞] (m/g, mL) to the model using the functions FindFit and NonlinearModelFit in Mathematica. The best fit is found to be In Figure 21, we plot the extrapolated prefactor c [N →∞] as a function of m/g and mL.

D Alternative approach to computing worldline instantons
In this appendix, we provide alternative calculations for the worldline instantons studied in Section 5, but for scalar quantum electrodynamics. Specifically, we first integrate over the proper time, and then perform the path integral, as in [3,10]. The calculations in Section 5 suggest that the spin factor does not make any contributions to the effective actions for the class of instantons we investigated (see (5.34) and (5.59)). Thus, the scalar results here can provide a direct check against the spinor results in Section 5.
Recall from (5.9) that Integrating out the proper time s, in the weak-field regime m Around a single-instanton x cl , the quadratic action is In the last line, we separated the operator M µν (u, u ) into the local part with L µν (u) and the non-local part. 10 Since we have not factored the endpoint x(0) = x(1) out of the path integral the fluctuation δx µ (u) satisfies the periodic boundary condition, Periodic Condition: δx µ (0) = δx µ (1). (D.8) Under this, we take M µν as self-adjoint, so that the spectral theorem applies. The quadratic operator M µν (u, u ) around a saddle point x cl typically has zero modes, which contribute to the path integral by some factor N 0 [42]. The non-zero eigenmodes contribute to the path integral (D.2) by the functional determinant with zero modes discarded, det M Here for each inequivalent single-instantons x cl,n , we have denoted the on-shell actions by S n,0 and the zero mode factors by N n,0 . We have approximated ( 1 0 duẋ 2 ) − 1 4 by its on-shell value at each instanton x cl,n .

D.1 The straight line instantons
The effective action (5.31) for the amplitude +g/2| + g/2 is We have adopted the formula for scalar quantum electrodynamics. From (D.9), for m This also means that the regime m 1 0 duẋ 2 1 corresponds to the semi-classical regime mL 1, i.e. the on-shell action S n,0 = nmL . To find the zero mode factor N n,0 , note that for our straight line instantons, the reparametrization mode redundancy u i −→ u i +δu i coincides with the translation of The solution is given by (see Figure 15) With the self-adjointness condition of L µν sorted out, we now find its eigenvalues and eigenmodes explicitly, for either choice of boundary condition. This is conveniently achieved by first diagonalizing L µν in (D.27): Of these, y (1),n and y (2),n satisfy the Dirichlet condition vanishing at u = ±1/2, and moreover vanish at u = 0. We need to check if y (i),n are orthogonal under the inner product (y (i),n , y (j),n ) ≡ 1/2 −1/2 y (i),n · y (j),n . It can be easily shown that the inner products are orthogonal-all save one: the inner product (y (3),n , y (4),n ) = b (n 2 − n 2 )π 2 ((−1) n+n − 1), n = n , (D. 40) cannot be made vanishing. This means that the modes y (3),n and y (4),n are not orthogonal when n and n are not both even or both odd (n = n in particular). By the spectral theorem, this implies that L µν acting on the space spanned by y (3),n and y (4),n is not self-adjoint. Indeed, we can compute the obstruction (D.34) to the self-adjointness of L µν , We thus reach the conclusion that, L µν is not self-adjoint on the space Ω of periodic, first-order differentiable vectors with the inner product (D.29). The corresponding functional determinant is then ill-defined. To avoid this problem, we therefore choose to impose Dirichlet condition on the fluctuations by pulling the endpoint x(+1/2) = x(−1/2) =x out from the path integral and integrate over it in the end. This space is spanned by y (1),n and y (2),n in (D.38). Since they vanish at u = ±1/2, and also at u = 0, they satisfy the self-adjoint condition (D.34). Thus, in the Dirichlet problem, L µν has no zero modes and has eigenvalues We can also reproduce (D.44) by the Gel'fand-Yaglom method [43], similar to what was done in [3].