{'fraction_non_alphanumeric': 0.03479490806223479, 'fraction_numerical': 0.010264699939381693, 'mean_word_length': 4.40896174863388, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variational Quantum Electrodynamics

A variational method is discussed, based on the principle of minimal variance. The method seems to be suited for gauge interacting fermions, and the simple case of quantum electrodynamics is discussed in detail. The issue of renormalization is addressed, and the renormalized propagators are shown to be the solution of a set of finite integral equations. The method is proven to be viable and, by a spectral representation, the multi-dimensional integral equations are recast in one-dimensional equations for the spectral weights. The UV divergences are subtracted exactly, yielding a set of coupled Volterra integral equations that can be solved iteratively, and are known to have a unique solution.


I. INTRODUCTION
In the last years there has been a renewed interest on variational methods for gauge theories [1][2][3], because of the relevance of non-Abelian gauge theories, that are known to be asymptotically free. The high energy asymptotic behavior of these theories is known exactly, which is one of the most important requirements for a viable variational approach to quantum field theory. On the other hand, important issues like quark confinement and the low energy phase diagram of QCD still lack a consistent analytical description because of the strong coupling that rules out the use of perturbation theory.
Recently, a new higher order extension of the GEP has been proposed [21], based on the method of minimal variance [22,23], and has been shown to predict nontrivial results even for fermions. The variance of the interaction contains second order terms and seems to be suited for dealing with the minimal coupling of gauge theories.
In this paper we explore the potentiality of the method of minimal variance by a study of the simple U (1) gauge theory with an interacting fermion, i.e. quantum electrodynamics (QED). An important merit of the method, shared with other techniques like the PGEP, is the paradox that the standard formalism of perturbation theory is used, while retaining a genuine variational nature, without the need of any small coupling. In fact, we do not assume that the coupling is small, but at any stage we check that expanding the results in powers of the coupling, the standard known properties of QED are recovered in the phenomenological weak coupling limit.
The issue of renormalization is addressed, and the standard renormalization scheme of perturbation theory is modified in order to obtain finite stationary conditions for the optimized propagators, which emerge as solutions of a set of coupled integral equations. Their numerical solution would be a first step toward the study of more complex non-Abelian theories in the strong coupling limit. However, even in this simple case, the numerical solution might not be so straightforward and seems to need some more effort. In that respect we discuss a method that is based on the spectral representation of the propagators. Under some assumptions, the multi-dimensional integral equations are recast in one-dimensional equations for the spectral weights, and the UV divergences are subtracted exactly, yielding a set of coupled Volterra integral equations, which can be solved iteratively, and are known to have a unique solution.
The paper is organized as follows: in Section II the method of minimal variance is described in detail for the simple case of QED; in Section III the problem of renormalization is addressed, yielding a set of finite stationary equations; in Section III the method of spectral representation is discussed, and the stationary equations are recast in a set of Volterra integral equations.

II. QED BY A GENERALIZED VARIATIONAL METHOD
The method of minimal variance [22,23] is based on a second order variational criterion that is suited to describe gauge theories with a minimal coupling like QED [21], where first order approximations like the GEP do not add anything to the standard treatment of perturbation theory [19]. The method has been discussed in some detail in Ref. [21]. Let us consider the basic U (1) gauge theory of a single massive fermion interacting through an Abelian gauge field where the last term is the gauge fixing term in Feynman gauge, and the electromagnetic tensor is F µν = ∂ µ A ν − ∂ ν A µ . We do not assume that the coupling e 2 is small, unless we would like to compare the results with the phenomenological QED. The quantum effective action Γ[a] can be evaluated by a shift a µ for the gauge field and is given by the sum of connected vacuum 1PI graphs [24] for the action S. Here the action can be split as S = S 0 + S I , and we define the trial action S 0 as where D µν (x, y) and G(x, y) are unknown trial matrix functions. By comparison with the definition of L in Eq.(1), the interaction can be written as the sum of three terms where ∆ µν (x, y) and g m (x, y) are free-particle propagators. Their Fourier transform can be expressed as where η µν is the metric tensor, andm = m − e a is a shifted mass matrix term. An implicit dependence on a is assumed in G, D, S 0 and S I . If the U (1) symmetry is not broken, in the physical vacuum a µ = 0 and the mass term becomesm = m. Of course, the trial functions G −1 , D −1 cancel in the total action S which is exact and cannot depend on them. Thus this formal decomposition holds for any arbitrary choice of the trial functions, provided that the integrals converge. The effective action Γ[a] can be evaluated by perturbation theory order by order as a sum of Feynman diagrams according to the general path integral representation of Eq.(2) By our decomposition of the action functional, we must associate the trial propagators G(x, y), D(x, y) to the free-particle lines of the diagrams, while the vertices are read from the interaction terms in S I . The three vertices that come out from the three interaction terms in Eq.(4) are reported in the first line of Fig.1. At any finite order, the approximate effective action does depend on the trial functions G, D, which must be fixed by a variational criterion. Several variational strategies have been discussed [21]: the variations δG, δD affect both S 0 and S I , and the optimal choice is the one that minimizes the effects of the interaction S I in the vacuum of S 0 , ensuring that the expansion makes sense even without any small parameter in the Lagrangian [25].
Denoting by X the quantum average the effective action can be written as where the zeroth order contribution can be evaluated exactly, since S 0 is quadratic and the remaining terms can be written by expansion of the logarithm in moments of S I which is equivalent to the sum of all connected 1PI vacuum diagrams arising from the interaction S I , as emerges from a direct evaluation of the averages by Wick's theorem. In our notation Γ n , V n , Σ n are single nth order contributions, while their sum up to nth order is written as Γ (n) , V (n) , Σ (n) , so that The effective potential follows as V (a) = −Γ[a]/Ω where Ω is a total space-time volume. We fix the trial functions by the method of minimal variance, requiring that the functional derivatives of the second order term V 2 are zero [21][22][23] In fact, by inspection of Eq.(10), the second order term can be written as where σ I is the variance of the Euclidean action S E I , That is obvious by Wick rotating, as the operator (iS) becomes the Euclidean action (iS) → −S E and the quantum action iΓ → −V /Ω. The method is based on the physical idea that in the exact eigenstates of an operator O the variance must be zero because OO = O 2 . For any Hermitian operator, the variance is a positive quantity, bounded from below, and the variational parameters can be tuned by requiring that the variance is minimal. In quantum mechanics the method is not very popular because the accuracy of the standard variational approximation can be easily improved by a better trial wave function with more parameters. In field theory, calculability does not leave too much freedom in the choice of the wave functional, which must be Gaussian. When the simple first order stationary condition fails, a second order extension can be achieved by the method of minimal variance [23] as discussed in Ref. [21]. Among the other variational strategies, we cite the method of minimal sensitivity [25] that would be equivalent to a search for the stationary point of the total second order effective potential V (2) instead of the single term V 2 . Actually, for the simple theory of a self-interacting scalar field, the total effective potential V (2) is unbounded and has no stationary points [20], while the stationary conditions Eq. (12) have been shown to have a solution [22], since the variance is always perfectly bounded.
The stationary conditions Eq.(12) are readily evaluated in terms of self-energy and polarization graphs, without the need to write the effective potential. In fact a general connection has been proven in Ref. [21] between the functional derivatives of the effective potential and the two-point functions, where the polarization function Π µν and the self-energy Σ ab are the sum of all connected two-point graphs without tadpoles. Explicit spinor indices have been inserted in the trial function G ab . First and second order twopoint graphs are shown in Fig.1.
Making use of Eq. (15) and Eq.(16) the stationary conditions Eq.(12) can be written as The first order two-point functions are given by a single tree graph each, as shown in Fig.1. Making use of the explicit form of the vertices in the interaction Eq.(4) we can write The proper self-energy and polarization contain one second order term each, the one-loop graphs of Fig.1 Σ These would be the usual proper two-point functions of QED if the functions D and G were replaced by the bare propagators ∆ and g m . The total second order contributions to the two-point functions follow by the sum of all second order graphs in Fig.1 where matrix products have been introduced in the notation. The stationary conditions Eq.(17) then read and are a set of coupled integral equations for the trial functions G, D. Their solution is equivalent to the optimization of an infinite set of variational parameters.
Before proceeding further, it is instructive to examine the first order approximation. The first order GEP is obtained by imposing that the first order effective potential V (1) is stationary. The general relations Eqs. (15), (16) for n = 1 give 1 where V (1) = V 0 + V 1 and the stationary conditions are equivalent to the vanishing of first order self-energy and polarization. Inserting the explicit expressions of Eq. (18), the stationary conditions of the GEP yield the trivial result D = ∆ and G = g m . Thus the GEP is equivalent to the free theory, and any meaningful variational approximation requires the inclusion of second order terms at least. By insertion of the explicit expressions for the first order functions Eq.(18), the second order coupled integral equations Eq.(21) can be recast as where the proper functions Π ⋆ 2 , Σ ⋆ 2 are given by Eq. (19). While this result resembles the simple lowest order approximation for the propagators in perturbation theory, it differs from it in two important ways: the presence of a minus sign in front of the second order term, and the functional dependence on the unknown propagators D, G in the proper functions in Eq. (19). Because of this dependence, the stationary conditions are a set of coupled integral equations, and their self-consistent solution is equivalent to the sum of an infinite set of Feynman graphs. In fact, despite the appearance, the stationary conditions are not a second order approximation of an expansion in powers of the coupling e 2 , but they make sense even when the coupling is large as they derive from a variational constraint on the variance.
Once the best trial functions are determined, as solutions of the coupled integral equations Eq.(23), perturbation theory can be used for determining higher order corrections with the optimized interaction S I and zeroth order propagators given by the solutions G, D. For instance, the second order propagator G (2) can be obtained by standard Feynman rules. We assume that the U (1) symmetry is not broken, and a = 0 in the physical vacuum. In terms of the proper self-energy and by inserting the explicit expressions for the first order self-energy Σ 1 = G −1 − g −1 m and the bare propagator g m , we find 1 As discussed in Ref. [21] the general relations Eqs.(15), (16) hold even for n = 0, 1 provided that we define Σ 0 = iδV 0 /δG, Π 0 = −2iδV 0 /δD, and take Σ −1 = 0, Π −1 = 0.
which looks like the standard one-loop result of QED, but differs for the functions G and D that must be inserted in the one-loop Σ ⋆ 2 in Eq. (19) instead of the bare propagators g m , ∆. If we expand the stationary conditions Eqs. (23) in powers of the coupling e 2 , take the lowest order approximation G ≈ g m , D ≈ ∆, and substitute back in the one-loop proper self energy Σ ⋆ 2 , then Eq.(25) becomes exactly equal to the one-loop propagator of QED. In fact, we can state that the variational method agrees with the standard results of perturbation theory when the equations are expanded in powers of the coupling. Thus, in the phenomenological limit of weak coupling the method of minimal variance would predict the standard results of QED. On the other hand, a numerical solution of the stationary conditions Eq.(23) would allow a study of the strong coupling limit.

III. RENORMALIZATION
Any numerical solution of the stationary equations Eq.(23) requires a regularization of the integrals and renormalization of the bare parameters in the Lagrangian. One of the main advantages of the present formalism is its Lagrangian approach that allows for a formal use of standard perturbation theory, while retaining a genuine variational nature of the approximation that is non-perturbative and valid even in the strong coupling limit. Thus, the problem of regularization and renormalization can be addressed by the standard techniques of perturbation theory, at any order in the optimized interaction S I , assuming convergence as a byproduct of the variational method. We use the standard dimensional regularization scheme of QED and define renormalized fields and couplings where µ is an arbitrary energy scale, and the space dimension is d = 4 − ǫ.
Gauge invariance requires that Z e = 1/ √ Z A at any order. We can also define renormalized trial functions and write the action as Everything goes as before with the substitution so that defining new renormalized proper functions in ddimensional space the stationary conditions Eq.(23) now read As usual, we expand the differences (Z − 1) in powers of the interaction S I , and denote by δZ the lowest order non-vanishing contribution. In the optimized theory δZ must be small, and we may regard it as a small parameter in the expansion. While the first order approximation does not require any renormalization, we find a non-vanishing δZ in the second order approximation, and assume that Moreover, at the same order of approximation, we may neglect higher powers of δZ in the stationary equations. For instance, we may completely neglect δZ in the second order two-point functions Eq. (20), while retaining a first power of δZ in the first order two-point functions Eq. (18). That is equivalent to dropping the factor (Z Ψ /Z A ) 2 in the last term of Eqs.(32), which can be written, in a compact notation, as having inserted a renormalized g R that satisfies, up to first order in δZ, In the minimal subtraction scheme (MS), the constants δZ Ψ , δZ m , δZ A are defined by the requirement that the quantities inside the square brackets of Eqs.(34) are finite, and are given by the polar diverging parts of the one-loop proper functions. For instance, assuming that the U (1) symmetry is not broken, and a µ = 0 in the vacuum, by Lorentz and gauge invariance we can write and defining by A ∞ , B ∞ , Π ∞ the polar diverging parts of A, B, and Π, respectively, in the limit ǫ → 0, the renormalization constants follow and setting Σ ∞ = A ∞ + kB ∞ , the renormalized stationary conditions can be written in the simple shape which are UV finite and can be solved for the functions D R , G R . Notice that the renormalization constants in Eq.(38) are the opposite of the standard definitions in QED. That is perfectly reasonable, as the aim of the present renormalization scheme is a finite integral equation for the functions D R , G R that play the role of zeroth order propagators in the perturbation expansion. The equivalent of the one-loop propagator is the second order function G (2) in Eq. (25), obtained by perturbation theory as the sum of all Feynman graphs up to second order, with the free lines given by the optimized renormalized propagators D R , G R . As a result, if these propagators are finite, the function G (2) is not, while if we want to make the second order function G (2) finite we must renormalize backward, and the zeroth order functions D R , G R would acquire diverging renormalization factors as for the bare propagators in QED. That seems more evident if we evaluate the second order function G (2) in the following two steps. Suppose we obtained finite functions D R , G R as solution of the integral equations Eq.(39), and want to write the first order function G (1) by perturbation theory. We need the first order proper self-energy, which is given by Eq. (18), and can be written in our renormalization scheme, according to Eq.(30), as The first order function follows which contains the diverging term δZ Ψ . If we would like to make the first order function finite, we must add a wave function renormalization term δZ ′ Ψ = −δZ Ψ . This is a backward renormalization that cancels the previous renormalization, since the first order approximation just gives back the bare propagator. Next, for evaluating the second order function G (2) we need the second order proper self-energy, which is given by Eq. (20). Neglecting higher order powers of δZ A comparison with Eq.(34) shows that the renormalization constants must be the opposite of Eq.(38) in order to get a finite second order propagator. After having canceled the renormalization in the first step, an opposite renormalization is required in this second step, going from a first to second order approximation. This opposite renormalization agrees exactly with the standard renormalization of QED. A similar analysis can be done for the polarization function and the renormalization constant δZ A . Thus the apparent wrong sign of the renormalization constants in Eq.(38) is just a consequence of the different aim of the present renormalization scheme that renormalizes backward with respect to the standard scheme, in order to get finite zeroth order propagators.

IV. SPECTRAL REPRESENTATION
A numerical solution of the coupled integral equations Eqs.(39) would give a variational estimate for the optimized propagators D R , G R . These functions are just the zeroth order approximation in the optimized expansion, but nevertheless they are expected to contain important physical insight. As the total action S does not depend on the trial functions D, G, they could be freely chosen as arbitrary variational parameters, and are not required to satisfy any physical condition, apart from convergence of the integrals. Of course, we expect that even if the trial functions were unphysical in some respect, the optimized functions G R , G (1) , G (2) would progressively acquire a physical nature if the expansion makes sense. However, as for any variational problem, physical constraints might be imposed on the trial functions in order to make the problem more tractable. If we impose that the functions D, G must be the propagators of some physical theory, then their spectral representation can be used in the integral equations Eqs.(39). That would be a way to cancel the divergences exactly, before dealing with the numerical problem. Moreover the multi-dimensional integral equations would give rise to one-dimensional integral equations for the spectral weights.
We illustrate the method by a weaker approximation and restrict the gauge field propagator to its free-particle value in Feynman gauge D = ∆, while assuming for G the Källén-Lehmann spectral representation [26] Hereafter, we drop the subscript R everywhere as we are dealing with renormalized quantities. We assume that the U (1) symmetry is not broken and a µ = 0 in the vacuum, so that the renormalized free propagator in Eq.(35) reads g −1 = k − m where m is the renormalized mass. Basically,in Eqs.(39), we ignore the second equation, as D is not varied, and optimize the choice of the spectral weights f 0 (ω), ρ 0 (ω) and of m 0 by the first equation. A full numerical calculation would require the inclusion of a finite mass for the photon, to be sent to zero at the end of the calculation, once the IR singularity has canceled. That is not a major problem, and we ignore it at the moment for brevity.
As we prefer to maintain the pole at the renormalized mass m, we modify the MS renormalization scheme of Eqs.(38) a little. The first of Eqs.(34) can be written as which is finite if we take where, as before, the superscript ∞ indicates the polar diverging part in the limit ǫ → 0. The first of Eqs.(39) still holds with and can be written as which is the stationary integral equation to be solved.
With this choice the function G has a first order pole at k = m, with a finite residue Thus, the spectral representation can be written as where the lower bound in the integral has been set at the renormalized mass m because of the vanishing of the photon mass. Here the spectral weight functions ρ(ω), f (ω) must be regular in ω = m, having taken the pole apart in the first term. Before inserting the spectral representation in the stationary equation Eq.(48), we find useful to introduce the regularized function where Σ ⋆ ω and Σ ∞ ω are evaluated by insertion of the function G ω (k) instead of G(k) in the definition of Σ ⋆ , Eq.(31).
With this notation, the subtracted proper function can be written as where Σ m is the regularized function Σ ω of Eq.(51) evaluated for ρ = 1 and ω = m. With the same notation the stationary equation Eq.(48) becomes Taking now the imaginary part, we obtain an integral equation for the regular weight functions ρ, f By insertion of the propagators Eqs.(43),(52) in the first of Eqs.(31), the function Σ ⋆ ω is given by the well-known QED proper self energy [26,27] with the mass replaced by ω, and odd powers of the mass multiplied by ρ The external integration can be evaluated yielding where µ has been rescaled as µ 2 → µ 2 e γ /(4π) and α is the standard QED coupling constant α = e 2 /(4π). We immediately extract the diverging terms Σ ⋆ ω (m) and and subtracting according to Eqs.(47),(51), the regularized function can be written as with the functions a(k), b(k) that follow, in terms of the complex varaiableω 2 = ω 2 − iη, from the integral representation By an elementary integration, real and imaginary parts follow in terms of the Heaviside step function Θ(x), and of and the residue Z 0 that by Eq.(67) now reads The Volterra integral equations are known to admit a solution, which is unique, and can be numerically evaluated by iteration. Of course, a full numerical analysis would require some extra care for the regularization of the IR divergence. In fact the zeroth order functions θ 0 , φ 0 have a pole at k = m, which is the lower integration limit. The insertion of a finite mass for the photon would raise the lower limit to a higher value m + > m, and would remove the divergence.
In the weak coupling limit, up to first order in α, we obtain for Z 0 the standard result of QED It would be interesting to study the behavior of Z 0 in the strong coupling limit, by a numerical solution, as the vanishing of Z 0 would be the sign of the onset of a new vacuum without single particle excitations. We do not expect it to occur in the present case, as we are keeping D = ∆ fixed, and we are neglecting the pair excitations that would contribute to the polarization function. However, the technique can be extended to the study of the full set of coupled stationary equations that come out from the method of minimal variance. Including a spectral representation for the trial photon propagator D, the paired stationary equations Eqs.(39) could be studied numerically by the same technique, yielding non-linear coupled integral equations for the weight functions. The present analysis shows that, at least in the weaker approximation of a fixed D = ∆, the method of minimal variance yields a non-trivial solution.