Path Integral Quantization of Generalized Stueckelberg Electrodynamics: A Faddeev-Jackiw Approach

We build a setup for path integral quantization through the Faddeev-Jackiw approach, extending it to include Grassmannian degrees of freedom, to be later implemented in a model of generalized electrodynamics that involves fourth-order derivatives in the components of a massive vector field being endowed with gauge freedom, due to an additional scalar field. Namely, the generalized Stueckelberg electrodynamics. In the first instance, we work on the free case to gain some familiarity with the program and, subsequently, we add the interaction with fermionic matter fields to complete our aim. In addition to deriving the correct classical brackets for such a model, we get the full expression for the associated generating functional and its associated integration measure.


I. INTRODUCTORY ASPECTS
Considering the context in which the quantum equations of motion were further developed from classical dynamics, Dirac [1] described them with a formal language using the correspondence principle.He discovered that there is a correspondence between classical dynamics in phase space (described by dynamic functions and Poisson brackets (PB)) and quantum dynamics in a Hilbert space (described by self-adjoint operators and commutators).Some time later, Dirac made the first thorough analysis of constrained systems in a Hamiltonian framework, probably one of his main achievements being the classification of constraints into first and second class, since the former were recognized as generators of gauge transformations.Besides that, Dirac's work gave rise to the so-called Dirac brackets, which are a generalized version of the PB and the correct ones for that kind of system.This Dirac's approach was further extended to connect classical and quantum dynamics of constrained systems via canonical quantization, by the Dirac-Bergmann (DB) algorithm [2,3].
The need for a relativistic description of the interactions that take place in nature led to the development of a covariant language with gauge freedom [4,5].However, this language involves more degrees of freedom than necessary, which implies the presence of constraints.Faddeev was the first to formulate the connections between classical and quantum descriptions of physical systems with constraints in a functional formalism, and Senjanovic later extended these ideas [6].
The quantization of a gauge theory is only possible in the physical degrees of freedom, resulting in a loss of explicit covariance in the equations.Nevertheless, in order to maintain explicit covariance at the quantum level, Faddeev and others developed a method that introduces non-physical ghost particles [7].
Following the historical development of quantization procedures of constrained systems, the Faddeev-Jackiw (FJ) method revitalized the subject in the 1980s [8].This formalism takes a classical geometric approach based on the symplectic structure of the phase space, and it is only applied to Lagrangians with linear dependence in the velocities.The method was originally based on the generalized Darboux's theorem for pre-symplectic forms, which are characteristic of phase space of constrained systems, to identify the constraints of the theory as well as the coordinates that describe the associated true degrees of freedom.In these coordinates, a reduced phase space is defined, which is equipped with a symplectic form whose inverse allows us to determine the modified Poisson brackets.In the majority of the cases of interest, it coincides with the Dirac brackets.In fact, many studies have explored the equivalence between the FJ and DB formalisms, with the equivalence being proved in some but not all cases [9][10][11][12][13].It is worth mentioning that in this approach one does not have to go through the steps of Dirac's method [14].
Certainly, the FJ approach offers several advantages, such as not having to distinguish between types of constraints and bypassing the Dirac conjecture.This led to an increase in interest in the formalism, so that a few years later Barcelos-Neto and Wotzasek made a contribution to the method by proposing an alternative path that circumvents the use of Darboux's theorem and instead uses the emerging constraints to extend the phase space iteratively inserting Lagrange multipliers as new coordinates until obtaining a new set that allows defining a symplectic structure on it [15,16].
On the other hand, by studying the fundamental role of the Lagrangian in quantum mechanics, Dirac introduced the very first idea of what is known today as transition amplitude and its relation with the classical action [17].This idea was later taken by Feynman to develop its path integral formulation, in which the classical action plays a key role [18].Subsequently, Schiwnger postulates its (variational) quantum action principle in the Heisenberg description [19,20].As we will see below, the functional integration formalism will trigger several lines of our research in the study of the quantization procedure, since, despite its primarily canonical nature, the FJ method has also been used in path integral quantization [21,22].Now, taking into account that the wide study of models that involve not only gauge freedom, but also massive vector fields, in the context of higher-order theories [23][24][25][26][27][28], as well as the different treatments of them (such as the DB algorithm or the FJ symplectic analysis) have led to several investigations on the quantization procedure [29][30][31][32][33], it is plausible to complete these investigations and explore the subtleties of the FJ quantization approach in a new environment [34].
Namely, the Generalized Stueckelberg (GS) electrodynamics.Finally, it is worth mentioning the fact that the Ostrogadsky instabilities can be overcome by quantum radiative corrections [35,36].This article is organized as follows: In Sec.II we present a general background of the free GS model, subsequently, the FJ framework is detailed in Sec.III to be later applied to the free GS model in Sec.IV.Thereupon, the path integral quantization in the FJ perspective and its application in the free GS model are presented in Sec.V. Thereafter, looking to incorporate interactions to the previous model, in Sec.VI we extend the FJ theory to the case in which Grassmann-valued fields are present.With all these tools ready, we make in Sec.VII the FJ analysis and the corresponding path integral quantization of the GS electrodynamics.Finally, Sec.
VIII contains the epilogue of the investigations, with our conclusions and perspectives of future work.

II. THE FREE GS MODEL
It is well known that the Podolsky model describes an ensemble of massive and massless generalized photons by means of an Ostrogadskian structure of fourth order in derivatives.This implies an ultraviolet improvement for the case of an interacting theory [37,38].However, even with this refinement, the model still presents infrared divergences.Consequently, one of the goals of a suitable extension is to provide a low-energy regularization for the theory, which can be achieved by a usual Podolsky theory plus a mass term written with a Stueckelberg combination, responsible for keeping the gauge invariance.This combination of U.V. and I.R. mass terms can be obtained by a generalization of the Higgs mechanism [34].The analysis of the interacting version of this model for T = 0 in the context of the Debye screening discussion was performed in [39].
The free case is described by the following Lagrangian density It is important to make a comment on the different contributions in the full Lagrangian density L above.The first term corresponds to the well-known Maxwell Lagrangian density (describing the dynamics of the free electromagnetic field).On the other hand, the two first terms (i.e., L M + L P ) correspond to the (free) Podolsky electrodynamics, which is a second-order theory of a real vector field.Finally, L M + L S leads to the free Stueckelberg model, which describes a coupling between a massive real vector field and a real scalar field B, in order to maintain the gauge freedom.
Considering the limit m → +∞, we get the Stueckelberg-Proca theory improving the I.R. sector; on the other hand, for m → +∞ and M → 0, we recover the free Maxwell electrodynamics.Finally, for M → 0, we get the free Podolsky model leading to U.V completion.Therefore, the complete model ensures improvements in both the U.V. and the I.R. regime.
It is easy to check that the total Lagrangian density L remains invariant under the following local gauge transformations with Λ = Λ(x) completely arbitrary.Interestingly, in the next section, it is going to be proved that the pure gauge Stueckelberg field B does not contribute to the physical degrees of freedom, as it should be.
It is worth mentioning that if we use a suitable covariant gauge 1 , it is possible to eliminate the harmonic Stueckelberg field and that way get the following equation of motion after the introduction of an external source J µ in the Lagrangian in (1) by adding a term −J µ A µ .
For the particular case of a static charge at the origin, i. e. J µ (x) = δ µ0 q δ 3 ( x − 0), the induced vacuum average for the inter-particle potential becomes [39] A with r representing the modulus of spatial position x, and the pole masses are defined as being It is interesting to note that, depending on the values of m and M, a complex pole can be reached leading to a confining phase [40].Additionally, one of the main properties of this model is to present three distinct phases for the potential A 0 .A mixture of Yukawa-like terms for M < m 2 , an exponential decaying one for M = m 2 (resembling the magnetic field behavior inside superconductors), and an oscillatory phase associated with the regime M > m 2 .Just for comparison, if a Proca-like pole M is generated through thermal corrections for the usual interacting Podolsky model, it would be2 M ∼ 4πα 3 T .Therefore, the non-trivial regime is reached at an incredibly large temperature of order 10 23 K.
Although investigating several well-defined possibilities for extensions of physical theories is an interesting subject by itself, the changes must be constrained by experimental observations.For instance, the Podolsky mass is restricted to the range m 370 GeV , to be in accordance with electron-positron interaction phenomenology [41].Besides, there is a very stringent bound on the Proca mass M ≤ 10 −12 eV according to the analysis of samples from extragalactic radio pulsars [42].It is worth noting that the investigation into the origin of the Podolsky mass remains an ongoing subject of current research [43]. 1 The Lorenz-Podolsky-Stueckelberg gauge condition is given by 1 The variety of interesting properties convinced us that such a model can furnish useful insights on various aspects of field theory, if investigated by several phase space approaches such as DB formalism [29], Hamilton-Jacobi analysis [31] and the FJ framework [30,33].This program was implemented for several field theories, revealing important complementary achievements.

III. STANDARD FJ FORMULATION
Since our aim is to implement the Faddeev-Jackiw approach on the free GS model, we present its formal framework in the context of a field theory, following Barcelos-Neto and Wotzasek [15].
The starting point is a Lagrangian written in its canonical form, i.e., linear in the velocities in which ξ I denotes the independent 3 dynamical phase space fields and its coefficients θ I are the components of the so-called canonical 1−form (which, in general, depends on the fields) and H is the Hamiltonian.For computational purposes, it is useful to work with the Lagrangian density 4 with θ I = θ I (ξ ; ∇ξ ) and H = H (ξ ; ∇ξ ).The Euler-Lagrange equations obtained from the above Lagrangian are given by [30,33] ω IJ (x; y) ξ J (y)d 3 y = In the non-singular case, ω is the symplectic form on the canonical phase space and, since is not degenerate, possesses an inverse, with components ω IJ , obeying ω IK (x; z)ω KJ (z; y) this relation can be used to express the equation of motion (8) in terms of ω IJ ξ I (x) = ω IJ (x; y) δ H δ ξ J (y) On the other hand, the time evolution of an arbitrary dynamical (space) functional F can be obtained as follows Then, considering the Hamiltonian as the generator (in the classical sense) of time evolution, according to Ḟ .= {F; H}, leads us to the following definition of Poisson brackets from the above definition, one readily obtains the important result which relates the elements of the matrix representation of ω −1 with the PB of the canonical phase space fields.This is, in fact, one of the main advantages of the FJ formulation.
In the more general case, the exterior derivative of the canonical 1−form yields a presymplectic form, which implies that its matrix representation is singular.The fact that any degenerate bilinear form has a non-trivial kernel can be understood in terms of the eigenvectors corresponding to the null eigenvalue6 of the associated matrix; these eigenvectors are often called the zero modes of such a matrix.We will denote them as in which r is a label indicating the (possible) degree degeneracy of the null eigenvalue.Making use of the equations of motion it is straightforward to check that applying v r on the Hamiltonian lead to the (dynamical) constraints in which we have introduced the constraint densities G r , since G r actually are also spatial functionals.These constraints determine a surface on the initial phase space and we must express all the quantities only on such a surface.Therefore, we write In many situations, some of the fields ξ B are present in the Hamiltonian density by multiplying some of the constraints and consequently, they disappear when the constraints are imposed.Taking this into account, we have denoted the surviving fields as ζ I ′ .Due to the dynamical character of the constraints G r , we need to require their stability along the time evolution of the system; this is getting by introducing Lagrange multipliers fields The presence of the λ r implies an extension of the old constraint surface, with field-coordinates , with field-coordinates ({ζ I ′ }; {λ r }) which we denote by ξ (1)I .In matrix notation (which will be often used) ξ (1) = ζ I ′ λ r t .Then, considering the objects θ (1) [ξ (1) ] the above Lagrangian can be written as L[ξ (1) , ξ (1) ] = θ (1) Note that we have arrived at an identical expression in eq. ( 6) and so, we are able to repeat all the above procedures.The quantities with superscript (1) correspond to the so-called first iteration of the Faddeev-Jackiw algorithm.The program stops when we end up with a (non-degenerate) symplectic form ω (n) , since its inverse determines a Poisson structure on the (extended) phase space M (n) .
If at some stage (iteration), say m, of the algorithm we obtain zero-modes of the pre-symplectic form that do not give rise to new constraints (in the sense that they result in identically null tautologies), then the system has gauge freedom.To see that, let us take an arbitrary variation of eq. ( 6) IJ defined as in eq. ( 8).If, in particular, we choose δ as a gauge variation, we know that the LHS of the equation becomes zero; On the other hand, it can be shown that the m−iterated Hamiltonian H (m) , which takes into account all the constraints that arise among the algorithm, is gauge invariant, and therefore thus, we find that the vector field X G = δ G ξ (m)I (x) δ δ ξ (m)I d 3 x cancels out the gauge invariant Hamiltonian H (m) .To continue with the program, appropriate gauge conditions must be imposed by hand, replacing each identically null constraint.
Once we get a symplectic form ω (n) , we can determine the time evolution of all of the variables present in the (extended) phase space M (n) by means of eq. ( 8) The knowledge of the model phase space structure leads to a well-defined path integral measure built by means of the determinant of the matrix representation of ω (n) .The construction of the path integral is based on the Darboux theorem, seen in the discussion [22], for example.

IV. FREE GS MODEL: A FJ APPROACH
In order to proceed with the FJ analysis of the GS model, we need to rewrite its Lagrangian density in the canonical form.However, because this theory has higher-order derivatives, we must be careful when defining the canonical momenta.We must notice that the complete model is a second-order theory on the fields A µ (due to the sector L P ), whereas it is a first-order theory on the scalar field B. Thus, we must employ the Ostrogadsky approach [44] for the fields A µ , and the regular canonical construction for the field B.

(i).
Vector field A µ • Independent fields (Ostrogradsky recipe): • Momentum fields (Ostrogradsky recipe) In the above (A µ ; Π µ ) and (Γ µ ; Φ µ ) are canonical pairs.Besides that, note that we end up with the kinematic constraint Since we have introduced new fields, we can readily obtain their corresponding gauge transformations, from eq. ( 2) However, since Γ µ are regarded as being independent of A µ , their gauge variations must be unrelated (in the sense that they should not depend on the same parameter), then From the above, we end up with the following gauge transformations According to the Ostrogadsky procedure, and taking into account eq. ( 27), the corresponding Hamiltonian density is given by with After replacing these velocities in eq. ( 32), one gets So, from eq. ( 32) we recover the canonical form of the Lagrangian density Then, the phase space fields are identified as being In terms of these variables, from eq. ( 35) we construct the corresponding canonical 1−form Thus, by taking the exterior derivative of θ we get the functional 2−form whose functional matrix representation is From this expression we can easily get the column zero-mode v t = 0 ν 0 ν 0 j 0 j 1 0 0 , or in In arrangement with eq. ( 17), we must apply this zero-mode to the Hamiltonian to obtain the following constraint from which we identify the constraint density Now, we must impose this constraint on the previous system before starting with the first iteration of the FJ algorithm.To identify the remaining fields after imposing G we must look at the 1−iterated Hamiltonian density H (1) = H G =0 , which, up a surface term, results Then, taking into account eq. ( 20) and eq.( 22), we construct the 1−iterated canonical Lagrangian density as with λ = λ (x), a Lagrange multiplier field.Since Γ 0 disappears from the entire dynamics we have the surviving fields Below we present the other objects involved in the first iteration • Fields defining the functional phase space M (1) • Functional canonical 1−form: • Exterior derivative of θ ( 1) • Matrix representation 7 of ω (1) The above matrix has the following zero-mode with ϑ = ϑ (x) an arbitrary scalar function.As we know, v (1) give rise to a new constraint according to Due to the arbitrariness of ϑ we establish the constraint density as being: Then, we start the second iteration of the Faddeev-Jackiw algorithm • Hamiltonian density • Canonical Lagrangian density We should note that in this iteration none of the fields were suppressed, so ζ (1) = ξ (1) .
• Fields defining the functional phase space M (2) • Functional canonical 1−form • Exterior derivative of θ ( 2) The above matrix possesses the following column zero-modes whose analytic expressions are given by The constraints generated by these zero-modes should be However, after computing explicitly these expressions, tautologies of type 0 = 0 are obtained and, therefore, do not represent new constraints.Note that v (2) 1,2 H (2) = 0 have the same form of the second equation in eq. ( 23), with v (2) 2 both playing the role of X G .This is, clearly, a necessary condition for the equations of motion to be invariant under the following local transformations which are compatible with the gauge transformations presented in eq. ( 31), with Ξ ↔ ϑ and Λ ↔ σ .Since we have obtained two tautologically null constraints, we must impose two additional conditions 8 to replace them.For this task, we use the Coulomb-Podolsky-Stueckelberg gauge condition, which is consistently attainable 9 , which is the more general (non-covariant) condition in our an Ostrogadskian framework which in turn fix the dynamics of A 0 by means of its equation of motion, G Hence, we end up with these two extra conditions that will lead to an (invertible) symplectic form.By repeating the steps of the FJ algorithm, the following symplectic form is obtained 8 Properly, we use an (achievable) gauge condition and a non-dynamical restriction which arises from the consistency of the former. 9To see this, it is sufficient to replace in which all the differential operators just involve derivatives with respect to x and we have called The inverse of the functional above matrix is given by 10 10 It is important to mention that the inverse object of a (functional) symplectic form is called a (functional) Poisson bivector field and is written in the basis δ δ ξ I (x) ∧ δ δ ξ J (y) ; unlike the symplectic form, which is commonly written in the basis δ ξ Note that the rows and columns corresponding to the λ −variables (Lagrange multiplier fields) were separated by dashed lines.In this way, the correct Poisson brackets of the surviving dynamical phase space fields will be given by the entries in the first block.Thus, the non-vanishing equal-time Poisson brackets are The determinant of ω (3) turns out to be proportional to det G 4 , which is indeed the same value obtained by the DB algorithm in the approach of [34].It can be used to build the path integral measure.As well as in the present formalism, there is recent interest in a wide set of alternative phase space studies.For example, one may cite [45], also in the context of higher derivative theories, in which the extended Chern-Simons theory coupled to scalar matter is analyzed, revealing an interesting structure for the Poisson brackets.

V. PATH INTEGRAL QUANTIZATION IN THE FJ FRAMEWORK
Since the last iteration of the Faddeev-Jackiw algorithm provides a symplectic structure for the (extended) phase space M (n) , the corresponding Lagrangian density L (n) becomes non-singular and therefore the path integral approach can be implemented.In general, the generating functional of correlation functions in the absence of sources will be given by an expression of the form in which D µ is a functional measure of integration to be determined and N is just a normalization factor.At this point, it is worth recalling that if Ω is the standard 11 functional symplectic form, then it is written in the so-called Darboux coordinates, say Ξ.Then, if we perform a change in the phase space M (n) between the Ξ−coordinates and the ξ (n) ones, we will have 11 That is, the one whose matrix representation has the form 0 since this equation corresponds to the components of the multiplication of functional matrices, we can write the above relation as from this equation, we recognize the presence of the Jacobian functional matrix corresponding to the transformation Ξ → ξ (n) .Then, taking the determinant in both sides of the above equation On the other hand, we know that in the Darboux coordinates, the generating functional adopts the simple form then, the change Ξ → ξ (n) implies the measure transformation Thus, replacing this into (72) yields the expression for the generating functional with integrand in terms of the ξ (n) −fields Hence, by comparing with (68) we obtain the desired integration measure Therefore, we see that the Faddeev-Jackiw approach provides the correct measure for the functional integral Z.On the other hand, we know that the ξ (n) −fields have the following structure in which we have denoted by ζ M the dynamical fields; i. e. those remaining after having imposed all the dynamical constraints 12 which arose along the algorithm in the Hamiltonian density.However, these fields are not necessarily independent because there could have been non-independent fields in the canonical sector 13 , since this part is not evaluated on the constraints surface (unlike the Hamiltonian density).With this notation, we can write the n−iterated canonical Lagrangian density as Recall that the original proposal in the Faddeev-Jackiw approach was to impose the timeconsistency of the constraints; however, we made the artifice in (19) just for convenience's sake.
Then, we are free to reverse it Replacing this expression for L (n) in (74) we get with It is easy to check that det ω (n) does not depend on the λ −variables and therefore we can isolate the functional integral on these variables as follows Note that each functional delta can be putted in the form δ and since det( ∂ ∂t ) does not depend on the dynamical fields ζ , it can be absorbed by the normalization factor 14 ; i. e., N → N det( ∂ ∂t ).Therefore 13 With this we are referring to the term which contains the velocities in the canonical Lagrangian 14 This absorption is possible since det( ∂ ∂t ) can be expressed as a functional integral over (real) Grassmann-valued fields C and C as: exp − and there is no coupling between these auxiliary fields (commonly known as ghosts) and the dynamical fields after replacing into the expression for Z.
Taking into account that the 3−iterated Hamiltonian density is equal to we get the following expression for the generating functional As a first step, it is convenient to perform the following integral With this, the expression for Z reduces to 15 Since from the definition of ξ (3) one identifies the dynamical fields are given by Note that since det(G 4 ) does not depend on the fields in the integration measure, could be absorbed by the normalization factor; however, we will keep that factor.On the other hand, it results It is easy to check that the corresponding Jacobian is equal to 1 and then The following step consists in performing the following integral On the other hand, consider the following integral With these considerations, the expression for the generating functional becomes Now, we perform the following integral The relabeling σ → A 0 allows us to write ( Ȧ j − ∂ j σ ) → ( Ȧ j − ∂ j A 0 ) ≡ F 0 j , and therefore: On the other hand, the integration over the field π yields With this, we arrive at the following form for By performing the integral over the fields Φ j , we get thus, we arrive at the following expression which in turn can be written in a more suggestive form Although this amplitude is correct, it is not manifestly covariant.Then, in order to turn it into a covariant form, we should take into account the following: (i).We can write the involved gauge condition in a more convenient way: (ii).Considering this authentic gauge condition, we can make use the Faddeev-Popov (FP) procedure and inserting the following identity object into any path integral with functional dependence on the fields A µ and B: with A µ and B being the transformed gauge fields according to (2).
(iii).Using the fact that the gauge invariance of the FP-determinant, the action and the functional measure involving the gauge fields allows us to reverse the gauge transformed field, we end up with an expression of the form: in which the braced first factor on the RSH is a infinite constant corresponding to the socalled gauge group volume 16 .Also note that the argument in the determinant correspond to the differential operator G, previously defined, and due to its independence from the fields, it can be separated from the functional integral.
With this, the final form for the generating functional reads: The above expression reveals two very important aspects, proper of an Abelian gauge theory: First, the extra factors coming from the initial functional measure end out being dropped due the covariant-like fixation of the (non-covariant) gauge.Second, the resulting expression coincides with the naive gauge redundant ones, up to the gauge volume group and the factor det(G).

VI. FJ THEORY INVOLVING GRASSMANN-VALUED FIELDS
In this section, we present an extension to the Faddeev-Jackiw technique for the case in which we are dealing with Grassmann variables [46].As in the previous case, the construction will be done in the context of classical fields.As we saw, the ξ −variables correspond to the ordinary phase space fields (involving, in most cases, both configuration and momentum fields).
In the present case we also will consider anticommutative (Grassmann-valued) fields, which will be called η−variables.Due to the (degree 2) nilpotency of the Grassmann variables, the terms in the Lagrangian density containing velocities of this kind of fields will be always linear in them, and consequently, the conjugate canonical momenta of these variables turn out to be proportional to the anticommutative configuration fields and so, the η−variables only contain configuration fields.In order to distinguish the ξ −variables from the η−variables we may adopt the standard terminology in Grassmann calculus: The commutative objects are said to have even Grassmann-parity (|ξ I | = 0), while the anticommutative objects have odd parity (|η α | = 1); alternately, we can use the physics terminology, in which we refer the ξ −variables as bosonic and the η−variables, as fermionic.
In this approach we define the unified variables z A .= ({ξ I }; {η α }), or in matrix notation: z st = ξ η .The starting point is, as we know, the Lagrangian written in its canonical form17 this Lagrangian always have even Grassmann-parity (|L| = 0) and as a consequence The corresponding equations of motion result From the explicit form of ω AB we get the following properties: In our convention, the functional exterior derivative acts on the left, i.e.: and has odd parity (|δ | = 1), the later corresponds to the Bernstein-Leites convention Under this convention, computing the exterior derivative of the canonical 1−form θ and hence ω = −δ θ .This relation is in apparent contradiction with (10), due to the negative sign; nevertheless, we must take into account two details: • In the present convention, the ordering of the contractions of the component indices with those of the base elements is contrary to the ordinary case.As a consequence, in the above notation, δ z A (x) and δ z B (y) correspond to the A−th row and the B−th column of the (super)matrix representation of ω, respectively.
• Since we have made the labeling A ! = (I, α) and B ! = (J, β ), we have18 Ignoring the odd variables z α and z β in (109), we stay only with the first term, and hence, according to the previous item, we recover the result of the ordinary case.
In the non-singular case ω is an even supersymplectic form and the components ω AB (x; y) of its inverse are such that Note that, by comparing with (106), the ω AB are antisupercommutative, whereas the ω AB are skew- supercommutative.Like in the ordinary case, the ω AB allows us to introduce the super Poisson Brackets of two dynamical superfunctions F and G as being which in turn implies In the singular case, the 2−form ω has zero-modes v r = v A r (x) δ L δ z A (x) d 3 x; which yields the constraints when applied to the Hamiltonian as in (17).A direct consequence of this is that the constraints (and also the constraint densities) always have even Grassmann-parity and therefore also the corresponding Lagrange multipliers λ r that allow introducing the consistency of those in the dynamics.The rest of the algorithm is essentially identical to the ordinary case but with the caveat that, since |λ r | = 0 the must be inserted in the phase superspace M (1) at the end (or beginning) of the even sector and not simply at the end of the old supercoordinates, i. e.
Besides that, it is important to take into account that since this case carries a Z 2 −graded structure, we are going to deal with (even) supermatrices and therefore we have to use the well-established rules for them.With these new issues, we will proceed to determine the canonical structure of the interacting case of the Podolsky-Stueckelberg model with massive Dirac fermions.

VII. GS ELECTRODYNAMICS (INTERACTING CASE)
The Lagrangian density for the interacting case can be obtained from the minimal coupling prescription in which L F is the Lagrangian density presented in (1) and the other terms correspond to Dirac Lagrangian density with the gauge covariant derivatives being e the electric charge of the fermion.Also, m e is the mass of the fermion.Thus, we arrive at the following Lagrangian density with L D being the free Dirac Lagrangian density.In our approach, each component of the Dirac field is a complex Grassmann-valued map, i. e.: ψ α : R 1,3 → C a and so |ψ α | = 1.The Lagrangian (116) is gauge invariant under the transformations in (2) and the additional finite local The corresponding odd canonical momentum fields are given by Due to the dependence between the fermion fields and its momenta, we establish the phase superspace fields as being z st = ξ ψ α ψ † α , with ξ defined in (36) 19 .The corresponding Hamiltonian density is given by Then, the canonical Lagrangian density results From this, we construct the canonical 1−form as being By using (109) we obtain the 2−form: whose supermatrix representation is20 which is singular with zero-mode v = δ δ Γ 0 (x) The constraint obtained by making G = vH coincides with the got in the free case since there is no coupling between Γ 0 and the fermionic fields.Thus, we can start with the first iteration of the Faddeev-Jackiw algorithm: • Hamiltonian density • Fields defining the phase superspace M (1) : • Canonical Lagrangian density • Canonical 1−form θ (1) = θ (1) • Negative of the exterior derivative of θ (1)   ω (1) = ω (1) The supermatrix representation of ω (1) turns out to be singular.Then, we have • Zero-mode of ω (1) • Generated constraint • Constraint density With this, we start the second iteration of the Faddeev-Jackiw algorithm • Hamiltonian density: • Fields defining the phase superspace M (2)   [z (2) ] st = ξ (2) • Canonical Lagrangian density • Canonical 1−form • Negative of the exterior derivative of θ (2)   ω (2) = ω (2) • Supermatrix representation of ω (2) It is easy to check that the supermatrix representation of ω (2) is singular, with zero-modes v (2) Nevertheless, these zero-modes do not generate new constraints, but tautologically null tautologies, and therefore two gauge conditions must be implemented.On the other hand, this fact suggests to us that the infinitesimal variations for the fermionic fields, along with the presented in (61), left invariant the Lagrangian density (116).In fact, we shall notice that (140) corresponds to the infinitesimal versions of (117), with Λ → σ .As in the free case, we choose the Coulomb-Podolsky-Stueckelberg gauge condition and the underlying equation By repeating the steps of the Faddeev-Jackiw algorithm we get the 2−form This yields the following supermatrix representation Following the rules of inverting an even supermatrix, we obtain: In which R and S are non-square blocks with odd entries given by: and with the same definitions G x m 2 previously adopted for the free case discussion.
From this, we recognize the following non-null equal-time super Poisson brackets involving fermion fields These brackets along the presented in (67) constitute the complete set of classical brackets of the Podolsky-Stueckelberg electrodynamics.
In order to make the functional quantization of the interacting case, it is important to recall that the generating functional is given by (81), but in this time considering the extra terms involving the fermionic fields 21 ) 21 It is worth mentioning that from now on we will consider the replacement ψ † α → ψα (which Jacobian does not affect the functional measure), for the sake of convenience.On the other hand, since in this case ζ st = ζ F ψ α ψα , the functional measure now reads: and The procedure is essentially the same as for the non-interacting case with few differences in the (89) and (90) steps, due to the presence of the term involving the ψ−fields.Thus, we get: As in (98), the expression obtained for the generating functional is non-covariant and so we could proceed as we did in the free case to obtain the desired manifestly covariant expression.However, a more elegant covariant expression for Z compatible with the Coulomb-Podoslky-Stueckelberg gauge is a sort of generalization of the so-called R ε −gauges, whose construction involves a family of Lorenz-Podolsky-Stueckelberg gauges-like characterized by a scalar function f , and multiplied by a suitable Gaussian weight parameterized by ε [34,47,48].
(i).Firstly, according to Faddeev and Popov, we introduce the following object (proportional to identity due to a factor depending on ε) into the expression for Z (ii).By performing a reverse gauge transformation, we get22 Reducing the second delta functional to δ Λ+G −1 (K∂ j A j − M 2 m s B) / det(G) and integrating over Λ, we get was made by implementing the correct functional measure of integration taking into account the Darboux theorem.Finally, after extending the FJ theory for Grassmann-valued fields, not only the complete FJ setup could be obtained for the GS electrodynamics (interacting case) but also the correct integration measure, which was implemented in the path integral quantization, following the same steps as the free case.
Regarding the future perspectives, the formal structure provided in this work can be applied to a wide class of systems: the Gross-Neveau self-interacting fermionic model [49], on the superconductivity and some other correlated descriptions of condensed matter.Moreover, it is also suitable to be applied in the context of supersymmetry [50] and bosonization procedures [51].Since our approach is capable to characterize non-linearly interacting systems, the study of non-Abelian BF models [52] can be conveniently performed since it is already a first-order theory.Another perspective is the study of the renormalization and radiative corrections for the interacting generalized Stueckelberg model.The calculation of the anomalous magnetic moment and comparison with the well-established experimental data can be a method to constrain the Proca and Podolsky masses 23 .Moreover, determining its associated partition function possibly leads to complementary constraints for the system's parameters due to thermodynamic phenomenology.

1 + m 2 ∂ µ A µ − M 2 m s B 2 d 4
Therefore, we end up with the final expression forZ Z = N det(K + M 2 ) exp − i 2ε x × × exp i S[A µ , B, ψ α , ψα ] D[A µ , B, ψ α , ψα ](159)In this Abelian theory case, the only remaining operator normalization can be written in terms of a Faddeev-Popov ghost pair with no coupling with the physical fields.Although it is not important for the present zero temperature case, this indeed has implications in the thermodynamics of the theory.Regarding the degrees of freedom, our final phase superspace has, separating bosonic and fermionic excitations, 20 + 2 × 4 dimensions.However, eliminating the four local excitations associated with the auxiliary fields, and considering the four constraints between the phase space fields, there are 12 + 2 × 4 remaining independent quantities.It corresponds to 6 + 4 configuration space local excitations.These are exactly the ones that correspond to two massive spin 1 fields corresponding to the poles and one Dirac fermion with its 4 degrees of freedom.It is worth mentioning that it reinforces the pure gauge character of the auxiliary Stueckelberg field.Despite being just an artifice, this auxiliary field turns the system into a gauge invariant one, opening the possibility of deriving Ward-Takahashi non-perturbative constraints even for the Boson mass renormalization constant.VIII.CONCLUSION AND FINAL REMARKS As it was possible to testify in the text, we fully investigated the quantization of the GS model in the path integral formalism in the Faddeev-Jackiw framework.First, we present some classical, quantum and thermal physics issues associated with the studied model.Continuing on, we present the standard formulation of FJ in the language of differential forms as well as exterior derivatives and instruct how to deal with the constraints that arise throughout the algorithm.Then we apply the FJ program to analyze the free case of the GS model, deriving the entire symplectic structure and obtaining the Poisson brackets of the model.The path integral quantization of the free GS model