Convergence of a discretization of the Maxwell–Klein–Gordon equation based on finite element methods and lattice gauge theory

The Maxwell–Klein–Gordon equations are a set of coupled nonlinear time‐dependent wave equations, used to model the interaction of an electromagnetic field with a particle. The solutions, expressed with a magnetic vector potential, are invariant under gauge transformations. This characteristic implies a constraint on the solution fields that might be broken at the discrete level. In this article, we propose and study a constraint preserving numerical scheme for this set of equations in dimension 2. At the semidiscrete level, we combine conforming Finite Element discretizations with the so‐called Lattice Gauge Theory to design a compatible gauge invariant discretization, leading to preservation of a discrete constraint. Relying on energy techniques and compactness arguments, we establish the convergence of this semidiscrete scheme, without a priori knowledge of the solution. Finally, at the fully discrete level, we propose a compatible explicit time‐integration strategy of leapfrog type. We implement the resulting fully discrete scheme and provide assessment on academic scenarios.

Model the particles are either fermions or bosons. The fermions are described by spinors, and the forces between them by connections on vector bundles. The connections model bosons. In quantum electrodynamics (QED) these are photons. In particular, in scalar QED, one models spin zero bosons that interact with photons. The equations are derived through a variational principle from a Lagrangian function. The form of the Lagrangian is determined by the symmetry (gauge) group of the theory. If, in addition, one demands the theory to be renormalizable, the maximal order of derivatives in the Lagrangian is reduced to one, see [25]. Through Noether's theorem (see [24]) one can show that gauge symmetry implies preservation of constraints on the solutions of the Euler-Lagrange equations. The Maxwell-Klein-Gordon (MKG) equations are an example of such equations modelling the interaction of an electromagnetic wave with a particle. The corresponding Lagrangian is invariant under a gauge transformation and the resulting Euler-Lagrange equations are a set of coupled nonlinear time-dependent wave equations verifying a charge preservation constraint, consequence of gauge symmetry. We focus on this set of equations for our study.
Gauge symmetry makes the theory consistent, one should thus strive to preserve this symmetry when discretizing the model. However, standard discretization techniques can break this symmetry. As an example, for the related Yang-Mills equation it was shown in [10] that standard FE discretizations produce approximate solutions that strongly violate the constraints. The field of structure preserving discretization addresses this issue. It is now successfully applied in various area of numerical analysis and applied physics (see e.g., [4,5,14,16,19,20,29] and references therein). The idea is that the preservation of the underlying geometric structure will reveal itself through both stability and good qualitative properties of the solutions.
Nonlinear wave equations such as, also, the sine-Gordon equations [1,15], provide a particularly interesting test-bed for such approaches.
In the MKG equations, the coupling between the electromagnetic field and scalar complex field describing the evolution of the particle arises through covariant derivatives. Conventional discretization methods, such as standard finite difference methods (FDM) or finite element methods (FEM), approximate the gradient part and the part containing the electromagnetic field separately. Consequently, the resulting approximation of the covariant derivative has no local transformation law, and the gauge symmetry is broken at the discrete level. This implies a violation of the constraint. The key to preserve the gauge symmetry is to approximate the covariant derivative directly, and not as the sum of the gradient and the product of the gauge field with the scalar field/spinor. This was achieved by Kenneth Wilson [32] in 1974 in the field of high energy physics. He was doing calculations on quarks on a lattice, and the theory he developed is now known as lattice gauge theory (LGT).
LGT approximates the covariant derivative in a consistent way and at the same time preserves the local gauge symmetry. The essence in the procedure is to localize the nonlocal terms arising in finite difference methods: the nonlocal terms which are to be compared, are parallel transported to a common reference point with the gauge potential. By doing this, the discrete theory becomes gauge invariant. In this context, one motivation of this work is to exploit this technique in a finite element framework to design gauge invariant schemes. As a result, we propose and analyze a numerical LGT scheme for the Maxwell-Klein-Gordon equations that satisfies the discrete constraint by preserving the gauge symmetry at the discrete level. Some LGT schemes have already been proposed in the literature. In two other articles [6,7], two of the authors already considered LGT discretization techniques for Maxwell-type equations. In [7], the convergence of the LGT scheme applied to pure Maxwell theory was studied. This was done by comparing the LGT scheme with the classical Yee scheme [33] which is known to converge. In [6], a LGT scheme and a standard Finite Difference discretization for the Maxwell-Klein-Gordon equations were compared. A discrete Noether theorem ensures that the constraint (charge) is preserved. However no proof of convergence of the scheme was provided in that paper. With a more applicative perspective, in [29], lattice schemes are also used (with a Yee-type scheme for the electromagnetic part) to investigate some numerical simulations for the massive Maxwell-Klein-Gordon equation in the context of QED. Here again no numerical analysis is provided.
LGT type discretizations have also been studied for other equations, in particular the Ginzburg-Landau equation [13].
Apart from LGT type discretizations, and regarding numerical analysis works in this context (MKG), in [8], two of the authors proposed a finite element semidiscrete energy preserving scheme in the temporal gauge where gauge symmetry is lost but the constraint preservation is recovered using a Lagrange multiplier. The complete semidiscrete numerical analysis is provided. Let us also mention the recent study [17] considering Maxwell-Klein-Gordon equations in the Coulomb gauge. There, a discretization framework is proposed, based on Finite Elements in space and a modified Crank-Nicolson scheme in time that is energy preserving but not constraint preserving. A complete convergence study and some academical test cases are provided. To the best knowledge of the authors, these are the only studies where a complete numerical analysis work is tackled in the precise context of Maxwell-Klein-Gordon equations.
In this article, we propose a discrete numerical framework for the massive and renormalizable MKG equations in the temporal gauge, based on conforming finite elements (FE) and LGT, aiming at preserving gauge symmetry at the discrete level. The Klein-Gordon (KG) part is discretized as in [6] with LGT techniques, while the Maxwell part is discretized using conforming Nédélec Finite Elements. The resulting semidiscrete scheme preserves a discrete constraint, namely the electric charge. We then prove the convergence of this semidiscrete scheme in two space dimensions. This is done with some inspiration from the methodology used in [8] using compactness arguments and energy principles. Discrete constraint preservation plays a central role in the proof to achieve the adequate bounds to obtain convergence. In a second step, we furthermore propose a fully discrete scheme based on a leapfrog time integration that also preserves the constraint at the discrete level. The fully discrete scheme is implemented and some academic numerical results are given for validation.
The paper is organized as follows: In Section 2, the continuous model is introduced from a variational point of view. In Section 3, we set the discretization in space using the lowest order Nédélec elements [23] on rectangles. The discrete gauge invariant Lagrangian is developed through both LGT and FE. Constraint and energy conservation are shown leading to the proof of the convergence of the scheme in Section 4. Finally, Section 5 numerically assess the fully discrete scheme.

THE MAXWELL-KLEIN-GORDON EQUATION
We first set the equations in a quite general setting of differential forms. Let M be a compact Riemannian manifold without boundary. The space of real-valued k-forms will be denoted Ω k (M). We will often identify one-forms and vector-fields. The real valued L 2 -product on differential forms on M is denoted ⟨⋅, ⋅⟩, and the associated L 2 -norm || ⋅ ||. Similar notations will be used for complex valued forms. All adjoints, denoted (⋅) * , will be taken with respect to these L 2 products.

The Klein-Gordon action
The unknown of the Klein-Gordon theory is a complex scalar field, t  → (t) (an element in Ω 0 (M) ⊗ C), and the corresponding action functional is given by where the dot represents time derivative and ∶ Ω k (M) ⊗ C → Ω k+1 (M) ⊗ C is the exterior derivative acting on complex valued forms. The third term is the mass term ( ≥ 0 is the mass) and the fourth term is the self-coupling term with ≥ 0.
Remark 2.1 The case < 0 will not be envisaged in this work and should deserve a specific study, since it leads to possible nonpositive energy.

The Maxwell-Klein-Gordon action
In this section, we formally explain the classical physical steps to achieve the gauge invariant MKG action. The MKG-equation is obtained by imposing a local U(1)-symmetry, that is, by demanding the action (or more precisely the Lagrangian) to be invariant under the transformation where is a real valued function on M, ∈ Ω 0 (M). The Lagrangian given in Equation (1) is clearly not invariant under this transformation. This is resolved by replacing the usual derivatives with covariant derivatives, that is, t ⇝ D ∶= The function is usually called the electric potential while A is called the magnetic potential. They are related to the electric and magnetic fields by the following equations and they transform as simultaneously with (2). This constitutes the gauge transformation of the field ( , , A) given as It is then easy to check that the following action To complete the action we add the U(1)-invariant Maxwell action given by with c the speed of propagation. The full MKG-action is then given by (cf. [18,25,26,28]) One can finally check that this action is invariant under the gauge transformation.

Euler-Lagrange equations
The stationary points of this action with respect to the different fields are given by the following Euler-Lagrange equations In a strong formulation, the Euler-Lagrange equations are also given by where ⋆ is the adjoint to and D ⋆ A = ⋆ − iqA. We see that the Euler-Lagrange-equations consist of two evolution equations given by (7)- (8) and a constraint equation given by Equation (6). Due to the local gauge invariance, Noether's second theorem can be applied to conclude that the constraint is preserved on the solution of the evolution equations (see [6]), which makes the equations consistent. See [24] for the continuous version(s) of Noether's theorem(s).
These equations together with the differential Bianchi identity (cf. [31]) which is satisfied by construction of the electromagnetic field from a gauge potential, constitute the complete set of the MKG-equations.
In the rest of the article, we consider unitary constants q = 1, c = 1. For the mass and self-coupling constant , we will consider that they are either 1 or 0 in the mathematical proofs to allow for variations on the system of equations considered. All the proofs are of course valid for nonunitary cases.

In the temporal gauge
We restrict the equations to the temporal gauge, = 0, allowed by the gauge symmetry for the rest of the paper. We also now focus on domains in R 2 . From this section, we also leave the differential forms notation behind (namely ) and use rather the usual notation grad = ∇, curl = ∇×, div = ∇⋅.

Strong form of the equation
In the temporal gauge, the strong form of the equations is to find (A, ) such thaẗ where D ⋆ A ⋅ = −div ⋅ −iA⋅ and ℑ is the imaginary part of a complex. The constraint is div(Ȧ) = ℑ(̇).

Notations and definition of weak solutions
We let S be a bounded contractible domain in R 2 with  1 boundary.
We use the classical notations for L p (S) spaces and Sobolev spaces W 1,s (S), H 1 (S), H 1 0 (S) (with seminorm | ⋅ | H 1 (S) ), and H(curl, S) is the space of vector potentials in R 2 considered as vector fields or one forms, with square integrable curl; the analogue space for the divergence will be denoted H(div, S). We also denote where A is the tangential component of A on S, and Time dependent spaces are defined as follows.
For closed intervals I ⊆ [0, T], (I; X) is the space of continuous functions from I to X, and (0, T; X) will denote ([0, T]; X).
We also define for 1 ≤ p ≤ +∞, the Bochner spaces L p (0, T; X) for X a Banach space as in [30]. We now give a rigorous sense to a weak solution in the temporal gauge and use a similar notion as in [8].
Definition 1 (E, A, , ) is said to be a weak solution of (80) in the temporal gauge, if We now turn to the discretization of this equation.

Finite element discretization
We discretize the spatial part of the continuous action. Let h > 0. We assume S to be a rectangular domain with a Cartesian mesh  h , and we will assume homogeneous boundary conditions. Furthermore, for (k, l) ∈ N × N,  k,l (C) is defined as the space of polynomials with complex coefficients in two variables (x 1 , x 2 ) ∈ R 2 with maximum degree k with respect to x 1 and l with respect to x 2 .
The discretization is based on three finite dimensional spaces Z 0 h , Z 1 h , and Z 2 h defined as (see [22,23] for properties of these spaces) We denote by h constructed with real valued functions. These spaces are equipped with real basis functions (w h n ), (w h e ) and (w h f ) respectively, which we choose as the tensor product of the one dimensional Whitney forms [22,23]. With these choices of basis functions, scalar fields h have degrees of freedom at the nodes of the mesh, h n , edge-vector-fields/one-forms A h have degrees of freedom at the edges of the mesh, A h e , while face-vector-fields/two-forms B h have degrees of freedom at the faces of the mesh.
Moreover, the gra and curl operators in 2D relate the finite dimensional spaces Z 0 h , Z 1 h , and Z 2 h , so that we have a complex: They also induce matrices G = (G en ) and R = (R fe ) in the chosen bases, such that where N, E, and F are the sets of vertices, edges and faces respectively. Thus, for a node element function h and an edge element function A h we can write Here, h n and (G h ) e are vertex and edge degrees of freedom, while A h e and (RA h ) f are edge and face degrees of freedom. The quantity (G h ) e represents the differential of the scalar field along the edge e, and has the form ( h m − h n ), where the edge e goes from node m to node n. Since we are considering rectangles, a natural orthogonal coordinate system can be associated to the mesh, and in such a coordinate system the expression ( h m − h n ) represents the differential in one of the two directions (see [6] for a more explicit formulation).
The notation e = {m, n} will denote the edge e which goes from node m to node n.

Scalar products and norms
The information about the shape and size of the rectangles is encoded in the mass matrices. For two edges e and e ′ , we define where ⋅ denotes the scalar product in R 2 . We also define in the cases k = 0 and k = 2, where s and s ′ are respectively two nodes or two faces of the mesh. We note that the matrices M h k are square, symmetric and positive definite. We also note that the nonzero entries in where h denotes the maximum diameter of the elements in the mesh.
These matrices M h k are the representative matrices of the L 2 -products in the corresponding bases w h . Thus, for example, for two edge element fields u h and v h written as where u h e and v h e are the edge degrees of freedom (DoF), we have We use (⋅, ⋅) to denote the real valued scalar product of vectors.

Norm estimates
One can establish estimates between norms and degrees of freedom as in the following Lemmas Furthermore, with p > 2, q > 2 and 1 p + 1 q = 1 2 , we have: There exists a constant C, independent of u and h such that Thus This estimate together with Hölder's inequality gives

Gauge symmetry
We can redo the formal derivation of the gauge invariant action at the discrete level. This amounts to understand how the natural analogue of the action at the discrete level would transform under the gauge and correct the terms to obtain a gauge invariant action.
With the notation developed in the previous sections, we can express the discrete version of the Lagrangian from (1) as wherėh and h are degrees-of-freedom vectors. As in the continuous model we want to impose a local U(1) gauge symmetry. We would like the theory to be invariant under the set of transformations (a) First, as in the continuous case, we replace the ordinary time derivative by the covariant time derivative, that is, (b) Next, we replace the mass matrix M h 0 with a mass lumped version H h 0 , which is both diagonal and positive definite. Its entries are given by Let u, v be continuous scalar functions, and u,v their vectors of nodal degrees-of-freedom. We define the associated bilinear form by ⟨., .⟩ 0,h (which gives the scalar product associated to where the points x appearing in the sum are the vertices of K and (⋅, ⋅) denotes the real valued scalar product for nodal DoF vectors.
For real l > 1, and all u ∈ H l (S) and v h ∈ Z 0 h , there exists a constant C depending on l such that Hölder type inequality. Furthermore, for p and p ′ such that 1 where for p > 1, one defines which is uniformly equivalent on Z 0 h to the true L p norm [3], that is, there exists a constant C As a conclusion, we have modified the term ⟨̇h,̇h⟩ as ⟨D h h , D h h ⟩ 0,h . As for the term ⟨̇h,̇h⟩ = (M h 0̇h ,̇h), we get a gauge invariant expression approximating (a) The mass matrix M h 1 defines the L 2 -product for fields with edge degrees-of-freedom. We lump this matrix with the same numerical quadrature as we used for M h 0 in Equation (12), as follows.
Define the scalar product ⟨., .⟩ 1,h on Z 1 h ×Z 1 h by a diagonal matrix H h 1 in the basis w h e indexed over the edges in the mesh. Let u, v be continuous vector fields/one-forms, and u,v their edge degrees-of-freedom. Then Here (⋅, ⋅) denotes the real valued scalar product for edge DoF vectors, (u,v) = ∑ e∈E ℜ(u e v e ). We observe that H h 1 is both symmetric and positive-definite. Consistency. We have the following error estimate (see [9,11]). For real l > 1, and all u ∈ H l (S) and v h ∈ Z 1 h , Henceforth, the mass matrix M h 1 is replaced with the above described mass lumped version H h 1 in the inner products. However, we observe that This can be resolved by inspiration from Lattice Gauge Theory (LGT) [12,27,32], that is, we make the replacement

3.2
Discrete formulation of the Maxwell-Klein-Gordon equation

Gauge invariant discrete Maxwell-Klein-Gordon action
The LGT inspired discretely gauge invariant Klein-Gordon action is therefore given by To complete the construction of the Maxwell-Klein-Gordon action, we add the Maxwell action The discretization of the Maxwell part of the action is well understood (see e.g., [2,[21][22][23]), so the gauge invariant action we are going to use is With the above considerations, where is the edge element uniquely defined by its edge degrees of freedom for e = {m, n}.
The Euler-Lagrange equations are given by the stationarity of the action, that is, By defining the electric field E h as and by a partial integration in time, the Euler-Lagrange equations read and We see that these equations consist of two evolution equations, Equation (22), and one constraint equation, Equation (23)  We suppose that the following initial conditions, Then we consider the following discrete initial conditions Furthermore, we suppose that they are chosen such that In the rest of Section 3.2, for simplicity of notations we drop the indices h, and consider the situation where h is fixed.

Constraint preservation
One important feature of this scheme concerns the constraint equation (23). The discrete MKG action (17) is gauge invariant, since both terms are. One can therefore use a discrete Noether's theorem to prove constraint preservation, in a similar manner as in [6]. We can also show this by a direct calculation. Proof. We start out by a differentiation in time of the left hand side of Equation (23), denoted . This giveṡ By the evolution equation for the electric field, Equation (22), we have Furthermore, and we can apply the nodal interpolator to both sides of this equality. Then (27) can be rewritteṅ The evolution equation for the Klein-Gordon scalar field giveṡ Since our scalar product is real valued, we obtain In a same manner, since one has and This giveṡ= where, This gives that Thuṡ= 0. This concludes the proof. ▪

Energy conservation
We define the energy of the system at any time with We will show through a direct formal calculation that this energy is preserved by the flow. (E, A, , ) solves the evolution equations (21,22). Then the energy is conserved.

Proposition 3.4 Suppose
Proof. The proof is a mere calculation. We havė Using (29), we deduce that So thaṫ Using (22), we finḋ Since the scalar product is real, one has that ⟨ , i ⟩ 0,h and ⟨| | 2 , i ⟩ 0,h vanish. This gives Using the computation done in the proof of Theorem 1, we conclude that ≡ 0, so that the energy is preserved in time. ▪

Choice of gauge and existence
We choose to work in the temporal gauge, that is, ≡ 0 (as in definition 1) and, for the discretization, Let T > 0. Since we are working on a finite dimensional space, we have local existence of solutions of (22), (23), and (24). Conservation of energy assures that the local solution is a global one: the discrete solutions are defined on the whole interval [0, T].

CONVERGENCE OF THE SEMIDISCRETE SCHEME
In the rest of the paper, C will denote a generic constant (independent of t and h). In the proof of convergence, we will need some results concerning the convergence of approximations. We state them here and postpone their proofs to Appendix A.1.

Lemma 4.1 Let Π 1,h be the edge interpolant. For p > 2, the following inequalities hold. There exists C > 0 such that for all
Furthermore there exists C > 0 such that for all t ∈ [0, T], The following result allows, by the constraint (23), to control the weak divergence of A h appearing in (47) and (48).

Lemma 4.2 Let p > 2. There exists C > 0 such that for all t ∈ [0, T]:
With these results at hand, we are ready to prove the convergence of the weak solution of (22).

Boundedness in the energy norm
The initial energy is bounded uniformly in h, as can be seen from (25), (26), and since it is conserved in time, we can immediately conclude that E h and ∇ × A h are bounded in L ∞ (0, T; L 2 (S)), that is, We can also conclude that h ,̇h and G A h h are bounded in time in the following sense.

Convergence of h
This is obtained in three steps. First, one bounds the H 1 (S)-norm of | h |, then one obtains a bound on the L p (S) (p > 2) norm of the gauge potential A h . In a third step, this is used to conclude that h is bounded in H 1 .

Boundedness of the H 1 -norm of | h |.
We would like to deduce that where A h e is the degree of freedom relative to the edge e. Since iA h e is purely imaginary, we have the following estimate so that by positivity of the diagonal matrix H h 1 , one can conclude that We recall that ∇Π 0,h | h | is the edge element vector field whose degrees of freedom are given by the vector G| h |. By previous estimates, L ∞ (0,T,L 2 (S)) ≤ C, which gives that Π 0,h | h | is bounded in L ∞ (0, T, H 1 0 (S)). In order to extract estimates on the H 1 -norm of h rather than on the H 1 -norm of its modulus | h |, one needs a control of the L p -norm of A h .
Boundedness of the gauge potential. Along the same lines as in [8], one can obtain a bound on A h in the L p norm. To this aim, we consider the discrete Helmholtz decomposition of A h We bound each part in L p . Bound on the discrete divergence-free part. The discrete divergence free part Å h is bounded in L p by the L 2 -norm of the curl of the gauge potential, We won't give any details on this estimate as it can be extracted exactly as in proposition 2.5. of [8].
Bound on the gradient part. One has from the constraint Equation (23) We would like to bound p h in L ∞ (0, T; W 1,p (S)) for p > 2. By classical estimates that can be found in [3], we have where 1 p + 1 p ′ = 1, with p > 2 and p ′ < 2. Let us choose q ′ such that 1 Since 1 2 + 1 p = 1 q , from the energy bound (52), we deduce that ||̇h|| L 2 (S) ≤ C and || h || L p (S) ≤ C. This implies that By the Sobolev embeddings W 1,p′ (S) → L q′ (S), we deduce which implies that p h is bounded in W 1,p (S), and Remark 4.3 If > 0, we can also directly prove (without having to use the boundedness in the H 1 norm of | h |) that A h is bounded in L ∞ (0, T, L 4 (S)), using that in this case h is bounded in L ∞ (0, T, L 4 (S)).
Convergence of h . From this, we deduce that h is bounded in H 1 (S). Indeed, In the last section we obtained where 1 q + 1 p = 1 2 . This means that for h sufficiently small, h ≤ 0 with 0 > 0 given, We also make use of the following inequality As a consequence ) are the degrees of freedom of the product A h Π 0,h (| h |), and one can then conclude that and bounds obtained on both A h and | h |, we can conclude that Remark 4.4 Following Remark 4.3, if > 0, we can also directly obtain this result (without having to use the boundedness in the H 1 norm of | h |) using the boundedness of A h in L ∞ (0, T, L 4 (S)).
By following [8] and using the bound on the energy, we arrive the convergence of h in L ∞ (0, T; L p (S)) (up to a subsequence) using the compactness result from [30] and interpolation estimates on L p (S) spaces.

Convergence of the gauge potential A h
The convergence of the gauge potential is obtained by considering the discrete divergence free part and the gradient part of the discrete Helmholtz decomposition (53) separately. Discrete divergence free part. We follow again [8]. By a Kikuchi type result, one obtains the convergence of Å h in L ∞ (0, T, L 2 (S)). Interpolation estimates then give the convergence in L ∞ (0, T, L p (S)), since one has (54). Gradient part. The equation we are considering is (Equation (55)) h ) * , and one can then find, by the Riesz representation theorem, Let q > 2 be given. We choose r > 0 such that 1 2 < 1 r < 1 2 + 1 q . One has where 1 r + 1 r ′ = 1, and v h = P h (v) is the L 2 orthogonal projection on Y 0 h (which is stable in L r′ (S)). As already shown in a previous section, this implies ||f h || L ∞ (0,T;L r (S)) ≤ C||̇h|| L ∞ (0,T;L 2 (S)) || h || L ∞ (0,T;L q (S)) .

Furthermore, ∀t ∈ [0, T],
From this one concludes that there exists w ∈ L ∞ (0, T; W −1,q (S)) such that where we used the compact embedding from L r (S) into W −1,q (S), and the compactness result of [30].
Define p ∈ W 1,q (S) as the unique solution of We would like to prove that In order to prove this, we use a version of the Strang lemma, that is, which is verified by construction of l, and we can conclude that p h → h→0 p in L ∞ (0, T; W 1,q (S)) for all q > 2.

4.3
The limit equation 4 and from energy bound we directly have in L ∞ (0, T; L 2 (S)) weak-*.

Limit equations
We denote by ′h and A ′h the vectors of the degrees of freedom associated to ′ h and A ′ h respectively. In the temporal gauge, the semidiscrete equations (22) read Study of Equation (57).
⟨̇h,̇′⟩ 0,h , and by weak-* convergence oḟh, Furthermore from (13) and uniform L 2 bound oṅh The convergence of the terms ∫ T 0 ⟨ h , ′ ⟩ 0,h and ∫ T 0 ⟨| h | 2 h , ′ ⟩ 0,h directly follows from the convergence of h in L ∞ (0, T; L p (S)) and the convergence for the test functions. We now study the convergence of the term It will be obtained in several steps. First, we decompose the quantity of interest We first concentrate on I 1 and I 2 . We need the two following Lemma. Their proofs are postponed to Appendix B.1.
Furthermore, from the bound on the energy we directly have the following.

Lemma 4.6
There exists C > 0 such that We can in a same manner obtain the analogous Lemma with h replaced by ′ h . Using Lemma 4.5,4.6, and 4.1 (and their analogous counter part for ′ h ), one proves that I 1 and I 2 converges uniformly in time to 0 as h → 0.
The estimation of J rely on the following decomposition To estimate these terms, one needs the four following Lemma. To ease the reading, we postponed their proofs to Appendix B.1.

Lemma 4.10
Let us briefly describe how the seven terms of the decomposition of J are treated. The first term J 1 is estimated using Lemma 4.7. Lemma 4.8 gives an estimation of J 2 . The terms J 3 , J 4 , and J 6 are estimated using 4.1 and bounds on the discrete solution. Lemma 4.9 gives an estimation of J 5 . Finally, since The proof is postponed to Appendix B.1. We write Using convergences summarized at the beginning of Section 4.3 and Lemma 4.1, we can prove that Furthermore, Each of the terms converges to 0 as h → 0. Indeed, for The limit equation. To conclude, which means that (A, ) is a weak solution of the Maxwell-Klein-Gordon equation in the sense of Definition 1.

NUMERICAL RESULTS
In this section, we provide some numerical results to assess the theory. To this end, we first propose a fully discrete scheme and then study two types of test cases in two dimensions.

Fully discrete setting
We consider a time-discretization that consists of a uniform subdivision of N T + 1 (N T ∈ N * ) points of the interval [0, T]. The time step will be denoted Δt ∶= T N T . We propose the following simple time discretization of leap-frog type.
For a given sequence B k , we will also use the following In this work, we focus on first numerical results and postpone a more thorough fully-discrete numerical analysis for a future work.

Constraint preservation
One can straightforwardly check that the constraint is verified.
Proof. If one expresses the discrete differential, one obtains But from the definition k+ 1 2 h and the scalar product, So that Similar arguments as in Section 3.2.3 apply to prove that the constraint is preserved even in discrete time. This implies that if the constraint is verified at initial time, that is, for k = 1. ▪

Numerical tests
We have implemented the proposed fully-discrete scheme using a dedicated module that we developed using the Finite Element library Firedrake. 1 Let us point out that these results are preliminary and used as first assessments of the method. We consider a two dimensional domain given by a square A art (t, x, y) = (cos( with = √ 2 2 + 2 . Then we define In this way, (A art , art ) is a solution of the following set of equations: with initial conditions We then use the discretization proposed above, add the external currents (J A art , J art ) on the right hand side, and compute the corresponding solution (A n h , n h ) n∈{0, … ,N T } . Doing so, we are able to compute the L ∞ L 2 error as We computed the solution for 4 types of (m, ) couples: (0, 0), (1, 0), (0, 1) and (1, 1). In Figure 1, we plot for several values of the number of points, the values of the respective error defined in (83). The error analysis is performed using mesh parameters that fulfill Δt = CFLh. Doing so, we compute the minimum order in time and space of the scheme. In Figure 2, we plot the results in logarithmic scales. We infer an order of convergence of 1 (a linear regression would even give slopes of ≈ 1.5). Local orders are presented in Table 1.    Although this is out of the scope of this work to obtain a theoretical proof of the convergence orders, one can however try to find, through our analysis, some indication of first order convergence for this test case using regular fields. Indeed the approximation of the discrete covariant derivative D A h h by the LGT gauge compatible discretization come at the price of the estimate such as (4.5) where we do not hope to obtain more than first order of approximation of the covariant derivative.
Regarding the constraint , in this case, one has instead of (23), the following equation.
Hence, except if currents are chosen such that the right hand side vanishes, we do not expect preservation of the constraint in this artificial case. We will not compute the constraint for this precise test case and postpone the tests to the next section.

Second test case and preservation of the constraint
For this second test case, we propose to initialize the system with an electromagnetic plane wave for A and a 2D-Gaussian initial profile for centered at the center of the domain. In other words, we consider a focalized density profile for the modulus of the complex Klein-Gordon field. In Figure 3,  we observe that the constraint is preserved up to machine precision. This is in accordance with the theory developed in previous sections. The energy is not preserved by the time-integration scheme, however relative variations of the energy from its initial value are relatively small (see Figure 4) and this variation gets consistently smaller as the discretization parameters decrease (see Figure 5).
To conclude, we represent in Figure 6 some modulus profiles and their evolution with time (the initial profile is represented in Figure 7). At this stage of our investigations, we cannot really interpret these profiles since we work on a toy academic model. This will be part of future works to improve the physical relevance of the setting considered here and the corresponding test cases. The way we use these profiles here is to show the effect of the introduction of mass and self-coupling term. Even if, again, the physical relevance of the values chosen are not discussed here, we propose to test the use of nonunitary values of the mass and self-coupling term.   What we observe is that the effect of mass and/or self-coupling occurs after some elapsed time and that the effect of self-coupling seems to be weak. We therefore also tested a situation where self-coupling is ten times bigger than mass. In this case, we significantly see self-coupling effects over mass on the profiles patterns (see Figure 8).

CONCLUSION
In this work, we studied the Maxwell-Klein-Gordon equation in dimension two and focused on the semidiscrete analysis of a discretization scheme based on lattice gauge theory (LGT) combined with Nédélec finite elements. The special feature of this scheme is that it ensures gauge invariance at the discrete level therefore respecting the geometrical structure of the original equation through the preservation of the constraint. Our analysis is not an a priori analysis, but concentrates on sequential compactness arguments without the a priori knowledge of a solution of the continuous equations. We use the discrete energy principle combined to constraint preservation to extract bounds in the appropriate spaces and extract convergent subsequences. At last, we implement a fully discrete numerical scheme based on a leap frog type time integration. We provide first academical test cases to extract convergence orders and assess the method. Future works include more relevant physical test cases and the fully discrete analysis of the proposed scheme.  We will use the following decomposition of G A h h in terms of its gradient part and its nonlinear part, that upon requiring h to be sufficiently small, is close to iA h h . Let e = {m, n} denote an oriented edge of the mesh. The following decomposition holds Proof. Let e = {m, n} be an oriented edge of the mesh, one has,