Hugo Ricateau

Leticia F Cugliandolo

Numerical solutions of Hamiltonian pdes: a multi-symplectic integrator in light-cone coordinates

We introduce a novel numerical method to integrate partial differential equations representing the Hamiltonian dynamics of field theories. It is a multi-symplectic integrator that locally conserves the stress-energy tensor with an excellent precision over very long periods. Its major advantage is that it is extremely simple (it is basically a centered box scheme) while remaining locally well defined. We put it to the test in the case of the non-linear wave equation (with quartic potential) in one spatial dimension, and we explain how to implement it in higher dimensions. A formal geometric presentation of the multi-symplectic structure is also given as well as a technical trick allowing to solve the degeneracy problem that potentially accompanies the multi-symplectic structure.

Introduction

Partial differential equations (pdes) are involved in the description of many scientific problems of interest. Quite commonly, understanding the behavior of the latter starts by solving a pde. Unfortunately, the ability to find analytical solutions to such equations is the exception rather than the rule.

Usual strategies to tackle pdes are to consider an approximation of the initial equation (by removing some non-leading terms, for example) or a particular domain of the parameter space (perturbative approaches, . . . ). However, these partial pieces of information can be insufficient to understand the behavior of the system in a satisfactory way.

It then becomes relevant to focus on approximate solutions, but this time, of the original equation and for the full range of variation of the parameters. This is exactly what one tries to achieve by using numerical methods. The question then arises as to how to control the numerical approximation.

To be more precise, let us take a time dependent process ρ(t), governed by a differential equation f ρ, ρ , ρ , • • • = 0 , where the indicates time derivative. In a finite-difference representation of this equation, the approximation process is quite well controlled and, at each time step, we know the order of the error made. Yet, a priori, we cannot predict the accumulation of these errors over many time steps and we cannot control the approximation made on the time-dependent solution (especially in the long-time limit).

The question can then be rephrased as, why should we trust a solution obtained with a numerical solver?

To address this question the standard procedure is to test, as precisely as possible, all the known properties of the problem. Firstly:

i. If a particular solution is available, we can easily check whether the numerical solution is in agreement with it. ii. In the same spirit, we can compare a numerical solution to the exact one for some particular choices of the parameters (by turning off all the interaction terms, for example).

Nevertheless, these two kinds of tests are not robust enough and nothing ensures that the numerical approximation will behave in the same way in a different regime (where no exact solutions are available).

The second kind of test is based on symmetries and conservation laws:

i. If the theory admits a symmetry group (not spontaneously broken) we expect the numerical solutions to be (as closely as possible) in agreement with the discrete analogue of this symmetry group (and, especially, the discrete part of it). ii. Due to the symmetry group, the theory can exhibit some conserved quantities that the numerical solutions should preserve as closely as possible.

Using such tools lets us recover the ability to keep control 1 arXiv:1702.06863v1 [math.NA] 22 Feb 2017 over the numerical approximations and authorizes us to trust (or not) the numerical solutions obtained.

Obviously, the accessibility to such validations is closely related to the structure of the pde and the method we are going to introduce will only be applicable to the Hamiltonian ones (i.e. arising from De Donder -Weyl -Hamilton equations) that will be defined later. Let us simply state for the moment that they are Lagrangian pdes (i.e. arising from an Euler -Lagrange equation).

We also have to stress that we will focus on the particular class of finite-difference methods. Other classes (spectral methods, finite-element methods, . . . ) are not as natural as finite-difference ones to fulfil the constraints involved by the preservation of the multi-symplectic structure to be explained below.

A finite-difference method is like a cooking recipe composed of two ingredients:

i. Firstly, a lattice that samples a bounded region of the support (e.g. the space-time manifold). ii. Secondly, a set of discretization rules that translate the continuous quantities to their lattice analogues. The continuous unknowns, defined on the space-time manifold, are sampled through the lattice. The discretization rules specify how to combine these samples if we want to compute derivatives, force terms, . . . Applying these rules to the equation of motion toggles from a pde to a set of algebraic equations (governing the behavior of the quantities defined on the lattice). Solving these algebraic equations leads to a set of values on the lattice nodes. This is a sampling of the solution, and in adjunction with some interpolation rules, an approximate solution of the pde is thus constructed. However, it has to be noted that these samples are not necessarily exact and, both the samples and the interpolation process are responsible for dissimilarities with the exact solution.

Many standard finite-difference schemes already exist and are often adequate. Each method has its own preferred application field. In the kind of problems we will be interested in, we need to control the very long time behavior (with respect to a characteristic time-scale in the system) and we need a procedure that minimises the error accumulated over a huge number of steps. For this reason, we need to develop our own numerical scheme that performs well over long time scales (even though we may have to make some compromise on its short time quality).

Generally, most methods are able to behave rather correctly on short time-scales. Therefore, the simplest and the fastest the method, the better it is in this regime. However, the problem complicates at long times, since two phenomena conspire against the performance of most strategies:

i. On the one hand, a priori, the reduction of the time step improves the quality of the approximation and then the quality of the solution but this obviously inhibits reaching long times. ii. On the other hand, at each step, some numerical truncation errors are induced by the finite precision of the numbers' representation in a computer. Such errors are generally inflated when the size of the step decreases and accumulate as the number of steps increases.

It results that, for a given final time, the precision is bounded from above, hence the necessity to chose a numerical method designed to behave correctly whatever the number of steps to handle.

Regarding mechanical systems (support is of dimension one, e.g. just time), it exists a very particular class of finite-difference integrators: the symplectic ones. They are well known (especially by researchers in planetary evolution) because of their very good capability to preserve the energy of Hamiltonian systems with a high accuracy even over long times [START_REF] Mclachlan | Geometric integrators for ODEs[END_REF]. Such integrators are based on the conservation of a very important (even central) structure of mechanical systems: the symplecticity of the phase space.

Generalizations to field theories (pde) appends some difficulties since the conservation of the energy is no longer rigid enough. Actually, the correct fundamental quantity to be preserved is now the stress-energy tensor. Its conservation is local (by opposition to the conservation of the energy, which is through a space integral) and hence more fundamental. Therefore, the symplectic structure is no longer adapted and needs to be generalized.

Multi-symplectic numerical integrators, introduced by

Bridges and Reich at the beginning of the 21 st century [START_REF] Bridges | Multi-symplectic structures and wave propagation[END_REF][START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity[END_REF][START_REF]Numerical methods for Hamiltonian PDEs[END_REF], generalize to pdes the concept of symplectic integrators. Applied to conservative pdes, multisymplectic integrators exhibit excellent local conservation properties (especially of the stress-energy tensor) and a very stable behavior for long time integrations [START_REF] Frank | Linear PDEs and numerical methods that preserve a multisymplectic conservation law[END_REF].

In the past fifteen years the subject has been widely studied [START_REF] Hong | The multisymplecticity of partitioned Runge-Kutta methods for Hamiltonian PDEs[END_REF][START_REF] Islas | On the preservation of phase space structure under multisymplectic discretization[END_REF][START_REF]Backward error analysis for multisymplectic discretizations of Hamiltonian PDEs[END_REF][START_REF]Conservation properties of multisymplectic integrators[END_REF][START_REF] Liu | Multi-symplectic Runge-Kutta-type methods for Hamiltonian wave equations[END_REF][START_REF] Ma | Highorder compact splitting multisymplectic method for the coupled nonlinear Schrödinger equations[END_REF][START_REF] Mclachlan | High order multisymplectic Runge-Kutta methods[END_REF][START_REF] Moore | Multi-symplectic integration methods for Hamiltonian PDEs[END_REF] and successfully applied to a broad variety of problems including the non-linear Schrödinger equation [START_REF] Chen | Multi-symplectic Fourier pseudospectral method for the nonlinear Schrödinger equation[END_REF][START_REF] Chen | Symplectic and multi-symplectic methods for the nonlinear Schrödinger equation[END_REF][START_REF] Chen | Multisymplectic splitting method for two-dimensional nonlinear Schrödinger equation[END_REF][START_REF] Hong | Multisymplectic Runge-Kutta-Nyström methods for nonlinear Schrödinger equations with variable coefficients[END_REF][START_REF] Islas | Geometric integrators for the nonlinear Schrödinger equation[END_REF][START_REF] Zhu | Symplectic and multi-symplectic wavelet collocation methods for two-dimensional Schrödinger equations[END_REF], the nonlinear Dirac equation [START_REF] Hong | Multi-symplectic Runge-Kutta methods for nonlinear Dirac equations[END_REF], the Maxwell equations [START_REF] Sun | Symplectic and multisymplectic numerical methods for Maxwell's equations[END_REF], the Klein -Gordon equation [START_REF] Chen | Symplectic and multisymplectic Fourier pseudospectral discretizations for the Klein-Gordon equation[END_REF], the Kortewegde Vries (KdV) equation [START_REF] Ascher | Multisymplectic box schemes and the Korteweg-de Vries equation[END_REF][START_REF]On symplectic and multisymplectic schemes for the KdV equation[END_REF][START_REF] Hu | Multi-symplectic method for generalized fifth-order KdV equation[END_REF][START_REF] Wang | Multisymplectic Euler box scheme for the KdV equation[END_REF], the Boussinesq equation [START_REF] Chen | Multisymplectic geometry, local conservation laws and Fourier pseudospectral discretization for the[END_REF], as well as the Zakharov -Kuznetsov (ZK) equation [START_REF]Multisymplectic geometry, local conservation laws and a multisymplectic integrator for the Zakharov-Kuznetsov equation[END_REF].

We introduce here a new finite-difference multi-symplectic method based on the centered box scheme. The latter was one of the first multi-symplectic schemes introduced [START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity[END_REF] and it has been proved that it is stable and possesses a number of desirable properties [START_REF] Frank | Linear PDEs and numerical methods that preserve a multisymplectic conservation law[END_REF]. However, it is not well defined locally [START_REF] Ryland | On multisymplecticity of partitioned Runge-Kutta methods[END_REF][START_REF] Ryland | On the multisymplecticity of partitioned Runge-Kutta and splitting methods[END_REF], so it requires a global solver and hence is not scalable. The idea we introduce in this paper is to use a rotated lattice in the light-cone coordinates, that restores locality of the algorithm.

More precisely, the aim of this work is to introduce a multisymplectic integrator that is apt to efficiently obtain the long-time dynamics of a field theory with high precision.

We organize the presentation in a practical way, first showing how the method performs compared to other ones in the market, and next discussing its theoretical justification.

We also discuss, without any assumption on the dimension of space-time, the problem of the degeneracy of the multisymplectic structure and we show how to solve it in the particular case of the non-linear wave equation.

The outline of the paper is the following.

In the first part of this paper, section 1, we will compare our multi-symplectic integrator in the light-cone coordinates (msilcc) to two standard methods:

i. On the one hand, a very basic scheme based on the Newton approximation of derivatives (this method is widely used by a broad community and proves to be preserving the multi-symplectic structure too, although it is explicit and thus of lesser quality). ii. On the other hand, the method proposed by Boyanovsky, Destri and de Vega [START_REF] Boyanovsky | Approach to thermalization in the classical ϕ 4 theory in 1 + 1 dimensions: Energy cascades and universal scaling[END_REF] constructed such that it exactly conserves the energy of the system (nonlocal conservation).

The comparison will be performed using the so-called λ φ 4 field theory in dimension D = 1 + 1. We will, in particular, study the local conservation (or not) of the stressenergy tensor. This example will allow us to emphasize the strengths and weaknesses of our method.

The second part of this paper, section 2, will be devoted to the definition of the necessary concepts, to exhibit the multi-symplectic structure, to deduce from it the local conservation laws (as well as the global ones) and, finally, to present how to rewrite the equations to prepare the implementation of the msilcc method. Throughout this section, the concepts and results will be illustrated through the example of the non-linear wave equation (whose λ φ 4 theory, used in section 1 of the paper, is a particular case).

Finally, in section 3 we will introduce our method in detail. We will also demonstrate the conservation properties as well as a review of the concrete solving methods of the algebraic equations involved in the numerical approximation (at this point, the necessity to work in light-cone coordinates will clearly appear). Again, the non-linear wave equation will be our working example.

A short conclusions section will close the paper.

1 The multi-symplectic integrator in light-cone coordinates (msilcc) versus the standard methods

This first section will be devoted to the comparison of our multi-symplectic integrator in light-cone coordinates (msilcc) to two standard methods. The first one is a very basic scheme based on the Newton approximation of derivatives, directly implemented in the Lagrangian formulation of the pde. This method is the simplest and, generally, the fastest to implement, so it is widely used and it is an unavoidable starting point. The second method, developed by Boyanovsky, Destri and de Vega [START_REF] Boyanovsky | Approach to thermalization in the classical ϕ 4 theory in 1 + 1 dimensions: Energy cascades and universal scaling[END_REF], is designed such that the total energy (a non-local quantity) will be exactly conserved whatever the configuration of the field or the size of the integration step.

These two methods will be presented in detail throughout this section, while the msilcc method will be detailed in the next sections. Now, first of all, let us introduce the model which will support the comparison.

1.1 The λ φ 4 theory in 1 + 1 dimensions

1.1.

The equation of motion

The comparison will be preformed on the so-called

λ φ 4 model 1 in dimension D = d + 1 = 1 + 1.
The unknown is the real dynamic field, φ(x, t), governed by a second order, non-linear, pde (the equation of motion):

φ = ∂ 0 2 φ -∂ 1 2 φ = -V (φ) = -φ 1 + φ 2 , (1.1)
where x and t are space and time respectively,

∂ 0 = ∂ /∂ct, ∂ 1 = ∂ /∂x, c is a characteristic speed (e.
g. the speed of light) that we set to one, c = 1, and the derivative of the potential, V , is given by V (φ) = ∂V /∂φ. The λ appearing in the name of the model is the parameter of the non-linear term and it has been set to one. The other parameter in the potential, the one that accompanies the quadratic term, has also been set to one in such a way that the potential has only one absolute minimum at φ = 0.

There is no exact general solution to this equation. Nevertheless, some particular solutions can be obtained in terms of Jacobi elliptic functions [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF] and they can be useful as a first check of the accuracy of a numerical integrator (see section 1.5.3).

Boundary and initial conditions

As previously mentioned a finite-difference method can be decomposed in terms of two ingredients: the lattice and the discretization rules.

The notion of lattice is a bit ambiguous and needs to be clarified. First, let us suppose it to be a regular tiling (since there is, a priori, no reason to take a more complex structure). Moreover, the spatial part of the lattice should be finite. Otherwise the integrator would have to solve 1 which belong in the class of the non-linear wave equation with the potential

V (φ) = r 2 φ 2 + λ 4 φ 4 .
an infinite number of algebraic equations (with the same amount of unknowns), which is generally impossible.

Since the spatial part of the support is bounded, the solutions need to be constrained on the boundaries. In the following we will impose periodic boundary conditions (pbc) (even though this is not a requirement for our method) with a period of length L:

φ(x + L, t) = φ(x, t) , ∂ 0 φ(x + L, t) = ∂ 0 φ(x, t) . (1.2)
We use an initial condition that complies with the pbc:

φ(x, 0) = A sin 2π x L , ∂ 0 φ(x, 0) = 0 . (1.3)
Therefore, the total energy is

E exact = L 0 dx 1 2 (∂ 0 φ(x, 0)) 2 + 1 2 (∂ 1 φ(x, 0)) 2 + V (φ(x, 0)) , = A 2 π 2 L + L 8 + 3A 2 32 .
The initial amplitude, A, allows us to control the predominance of the non-linearity. For a sufficiently small amplitude, the non-linear part of the potential will be dominated by its quadratic part and the initial condition leads to a time-dependent solution that is close to the second eigenmode of the linear wave equation: A sin 2π x /L cos 2π t /L . Conversely, when A increases, the non-linear term becomes predominant and the behavior of the solution turns out to be much more complex.

Figure 1 represents the short time behavior of the solution obtained using the msilcc method for different values of A. As expected, for A = 0.1 the solution remains very close to the second eigenmode of the linear wave equation. For A = 3, the solution evolves in two ways: its characteristic time-scale decreases, and the amplitude of the oscillations becomes a little bit bigger than A (see in fig. 1, A = 3, the small circles at the center of the antinodes where the value of the field exceeds A). The impact of the non-linearity becomes significant. Then the larger is A, the shorter the characteristic time. The non-linearity is also destructing the structure of the eigenmode: when A increases the solution is more and more distorted.

The effect of A is twofold, it will allow us to explore the influence of the non-linearity and the effect of decreasing the quality of the sampling when the typical variation scale of the field becomes closer and closer to the lattice spacing.

The stress-energy tensor, its conservation and the charges

As previously mentioned, the most fundamental quantity that the theory shall preserve is the stress-energy tensor (see section 2.2.5 for its definition, the one of the charges as well as the proof of their conservation). For the λ φ 4 theory in 1 + 1 dimensions the stress-energy tensor is symmetric

T 00 = 1 2 (∂ 0 φ) 2 + 1 2 (∂ 1 φ) 2 + 1 2 φ 2 + 1 4 φ 4 , T 01 = T 10 = -∂ 0 φ ∂ 1 φ , T 11 = 1 2 (∂ 0 φ) 2 + 1 2 (∂ 1 φ) 2 - 1 2 φ 2 - 1 4 φ 4 .
Its local conservation is given by

∂ 0 T 00 + ∂ 1 T 10 = 0 , (1.4a) ∂ 0 T 01 + ∂ 1 T 11 = 0 , (1.4b)
or, in other words,

∂ 0 φ φ + φ 1 + φ 2 = 0 , ∂ 1 φ φ + φ 1 + φ 2 = 0 ,
which are satisfied as long as the equation of motion (1.1) holds.

Conversely, the numerical equivalents of these local conservation laws (eqs. (1.4a) and (1.4b)) will not be exactly satisfied. The violation comes from the fact that, in the discrete version of these equations, the term in brackets is not the discrete analogue of the equation of motion. This is precisely due to the fact that the discretization rules not always fulfil all the rules of differential calculus (Leibniz, . . . ).

Since these two quantities are non-vanishing they will allow us to control the quality of the numerical solution: a good numerical approximation should preserve, as closely as possible, the local conservation of the stress-energy tensor. These residues will be our first quantities of interest.

Let us now define the charges as

Q ν = L 0 T 0ν dx ,
where ν is either 0 or 1. These are global quantities. Integrating over space, the local conservation of the stressenergy tensor leads to the conservation of the charges,

∂ 0 Q ν = 0 .
Again, these quantities are not exactly conserved numerically and the resulting residues will be our second quantity of interest.

The testing conditions

In the previous sub-section we highlighted the quantities allowing us to examine the quality of a numerical approximation (of a Hamiltonian pde). Let us now introduce in which context they will be observed.

For each numerical method we will examine, through two situations, the error committed on the conservation of the stress-energy tensor (local) as well as the error committed on the evaluation of the charges (global). Firstly, we will have a look at how these errors behave as functions of A. We recall that A has a twofold effect: it affects the weight of the non-linearity, but also the quality of the sampling (since when A increases, the characteristic timescale decreases while the time-step remains fixed). So, we expect the errors to be lower at small A than at large A.

Secondly, we will fix A = 10 and observe how the errors behave as a function of time. We will explore both short and long time behavior.

Let us now define the symbolic operator ∆. When applied on a continuous equation, its effect is to extract the residue from the discrete analogue of the equation, and divide this residue by a characteristic quantity such that the result is not dimensional and the error can thus be compared to 1.

Before entering into the evaluation of the quality of the numerical methods let us present their construction in detail.

The Newton method

The Newton method is probably the simplest finitedifference method one can develop. Its ease of use and its efficiency make it a classic. However, we will see that it can be inaccurate and even unstable. For now, let us describe its construction.

Sampling the space-time manifold

The support of the theory, M, is a flat 1 + 1 dimensional Minkowski space-time manifold. Taking into account the boundary conditions and the fact that the method will be used as an integrator, M becomes a flat half cylinder:

M = T 1 × R + , x = j δ t = n δ δ δ
Figure 3: Lattice description of the space-time manifold in the Newton method. The lattice spacing in space and time are chosen to be equal to preserve covariance.

where T 1 = S 1 = R /L Z is the flat one-dimensional torus of length L, and where, without loss of generality, the initial condition is supposed to be given at t = 0.

The lattice, M , will then be taken as a regular tiling of M with, as generator, a square of width δ aligned with the space and time coordinates. Therefore, the lattice is defined by

M = δ Z /N Z × δ N = δ Z N × δ N ,
where N δ = L. The geometry of M is represented in fig. 3 and is nothing else than a square lattice. Now, the field, φ : M → R, can be sampled through the lattice as ϕ : M → R such that

ϕ j n = φ(x = j δ, t = n δ) ,
where n ∈ N and j ∈ 0, N .

By this sampling process, at a given time, we have switched from the infinite number of degrees of freedom of the dynamic field to a representation with only a finite number (N ) of degrees of freedom that can be used in a computer. This achieves the first step of the construction of the finite-difference approximation.

The Newton scheme

In order to complete the construction of the finitedifference scheme, the second step is to provide the rules that will indicate how to combine the samples of the field (i.e. the elements of ϕ) in order to obtain the physical quantities (and especially the equation of motion).

The derivatives of the field will be approximated using the Newton's rule which reads

∂ 0 φ(x = j δ, t = n δ) ≈ D ± 0 ϕ j n ≡ ± ϕ j n±1 -ϕ j n δ , ∂ 1 φ(x = j δ, t = n δ) ≈ D ± 1 ϕ j n ≡ ± ϕ j±1 n -ϕ j n δ ,
where the + (respectively -) stands for the forward (respectively backward) approximation. These two definitions (+ or -) are inequivalent and they express the fact that a finite difference can either represent the derivative at the end or at the beginning of the interval. This is exactly why Newton's method is so easy to implement but is also the reason for its inaccuracy (a finite difference should only represent the derivative at the center of the interval). However, the centered Newton's rule for the first order derivative (that combines D + and D -to involve the points # + 1 and # -1) will not be used since it leads to an inconsistent approximation of the second order derivative 2 .

The equation of motion can be approximated in two ways: using forward then backward rules or vice versa, i.e.

∂ 2 ≈ D -D + or ∂ 2 ≈ D + D -(either forward -forward
or backward -backward leads to an inconsistent approximation of the second derivatives 3 ). In both cases the discrete version of the equation of motion at time t = n δ and position x = j δ reads

ϕ j n+1 -2ϕ j n + ϕ j n-1 δ 2 - ϕ j+1 n -2ϕ j n + ϕ j-1 n δ 2 = -ϕ j n 1 + ϕ j n 2 .
This algebraic equation is explicit in ϕ j n+1 . The evolution of a given state can then be efficiently obtained using

ϕ j n+1 = ϕ j n-1 + ϕ j-1 n + ϕ j+1 n -δ 2 ϕ j n 1 + ϕ j n 2 .
The nodes of the lattice that appear in this equation are highlighted in fig. 3: the vertex in the left hand side of the equation is represented as while the vertices involved in the right hand side of the equation are represented as .

This concludes the definition of the Newton method. Let us now use these rules to obtain the discrete formulation of the stress-energy tensor.

2 To show this fact, let us define the centered Newton's rule as

D C ϕn = D + ϕn + D -ϕn 2 = ϕn+1 -ϕn-1 2 δ .
The second order derivative, that reads

D C D C ϕn = ϕn+2 -2ϕn + ϕn-2 4 δ 2 ,
thus leads to two independent sub-lattices (the odd one and the even one).

3 since

D + D + ϕn = ϕn+2 -2ϕn+1 + ϕn δ 2 ,
involves two unknowns (ϕn+2 and ϕn+1) and since

D -D -ϕn = ϕn -2ϕn-1 + ϕn-2 δ 2 ,
requires to have solved the neighbouring equation in space to get ϕn (which is incompatible with the periodic boundary conditions).

The energy and the stress-energy tensor

As the derivatives can be approximated in two ways (D + or D -), the stress-energy tensor can be defined in two ways too:

T 00 ± = 1 2 (D ± 0 ϕ) 2 + 1 2 (D ± 1 ϕ) 2 + 1 2 ϕ 2 + 1 4 ϕ 4 , T 01 ± = T 10 ± = -D ± 0 ϕ D ± 1 ϕ , T 11 ± = 1 2 (D ± 0 ϕ) 2 + 1 2 (D ± 1 ϕ) 2 - 1 2 ϕ 2 - 1 4 ϕ 4 ,
where the space and time labels were omitted. These two definitions are inequivalent but both of them are valid and lead to a residue (again omitting the n and j indexes):

D ∓ 0 T 00 ± + D ∓ 1 T 10 ± = 0 ± , D ∓ 0 T 01 ± + D ∓ 1 T 11 ± = 1 ± .
In practice the two definitions of these residues behave in the same way and the results will only present 0 = 0 + and 1 = 1 + . Obviously, the same reasoning can be applied to the conservation of the charges but we do not detail it here.

Energy conservation

The general treatment of the energy conservation will be exposed later, but let us briefly show how the total energy behaves. At time t = n δ (again, index n will be omitted) it can be defined in two ways

E+ = Q 0 + , E-= Q 0 -,
where Q 0 ± is the charge defined as

Q 0 ± = δ N -1 j=0 T 00 ± .
One can also envisage to combine these two definitions as

E ave = E+ + E- 2 .
These three definitions are represented as a function of time for A = 10 in fig. 4.

We observe that both E+ and Evary in time with an amplitude of the order of 5% of their time-averaged value. The amplitude of the variations for E ave is reduced to ∼ 1% due to a compensation of the errors in E+ and E-.

Nevertheless, E ave does not correspond to any discretization rule: it is the average of the energies obtained using different rules, which differs from the energy that would be obtained from the combination of the forward and backward rules (that, as already mentioned, leads to an incorrect approximation of the second order derivatives). We also notice that there is no apparent change in the amplitude of the deviations as time elapses up to a timescale at which the energies rapidly diverge. The divergence of the energies is directly due to the divergence of the solution which has been destabilised by the integration method.

We conclude that, although Newton's method is straightforward to implement, it is inaccurate (the total energy conservation up to 1% is not acceptable in most applications) and can even become unstable. Therefore, Newton's method will not be a good choice to integrate a field theory over a long time. These observations will be confirmed by the study of the local conservation laws presented in the following sections. Before presenting the local analysis, let us first introduce another finite-difference method.

1.3 The Boyanovsky -Destride Vega (BDdV) method

The Boyanovsky -Destride Vega (BDdV) method [START_REF] Boyanovsky | Approach to thermalization in the classical ϕ 4 theory in 1 + 1 dimensions: Energy cascades and universal scaling[END_REF] that we present here has been developed such that it exactly preserves the total energy of the system, making it a good candidate for long time integrations. However, as previously said, the conservation of the total energy is not the most fundamental principle for a field theory that should foremost locally preserve the stressenergy tensor.

The lattice

In the BDdV method the space-time manifold is rotated by π /4. More precisely, the space-time manifold, M, is unchanged, and the lattice, M , is still taken as a regular tilling of M. The generator is still a square of width δ, but aligned with the light-cone coordinates. Therefore, the lattice is defined by

x = j + 1+σn 4 √ 2 δ t = n δ / √ 2 δ δ j n
M = √ 2 δ j + 1 + σ n 4 , n δ √ 2 n ∈ N , j ∈ Z N ,
where

σ n = 2(n mod 2) -1 ≡ ±1 .
M is represented in fig. 5 and is nothing else than a square lattice in the light-cone coordinate system which correctly respects the boundary conditions.

Finally, in the same way as for the Newton method, the field, φ : M → R, can be sampled through the lattice as ϕ : M → R such that

ϕ j n = φ x = √ 2 δ j + 1 + σ n 4 , t = n δ √ 2 ,
where n ∈ N and j ∈ 0, N .

Exact energy preserving approximation

We now provide the rules that allow one to express the discrete analogues of the physical quantities, respecting the directions imposed by the lattice. Under this constraint, the derivatives are written along the light-cone coordinates as

(∂ 0 -∂ 1 )φ √ 2 x = √ 2 δ j + 1 + 3 σ n 4 , t = n δ √ 2 ≈ ϱ 0 ϕ j+ σn 2 n = ± ϕ j n±1 -ϕ j+ σn±1 2 n δ , (∂ 0 + ∂ 1 )φ √ 2 (• • •) ≈ ϱ 1 ϕ j+ σn 2 n = ± ϕ j n±1 -ϕ j+ σn∓1 2 n δ ,
where the • • • indicate that the field is evaluated at the same point on the space-time manifold as in the first equation.

The method now differs from Newton's since the discretization rules are not applied to the equation of motion but to the energy and the constraints imposed by its conservation are used to derive a modified discrete evolution equation. In the continuum limit this equation would be identical to the equation of motion, but in the discrete formulation it differs from the one we would have obtained had we directly applied the rules to the equation of motion. As a matter of fact, the idea introduced here is very deep but we will have the opportunity to come back to this later.

The local energy density (the 00 component of the stressenergy tensor) will then be approximated in two ways (following the same principle as for the two possible Newton approximations of the derivatives) and is given by

T 00 ± j+ σn 2 n = 1 2 ϱ 0 ϕ j+ σn 2 n 2 + 1 2 ϱ 1 ϕ j+ σn 2 n 2 - 1 4 + 1 8 1 + ϕ j n±1 2 2 + ϕ j n 2 + ϕ j+σn n 2 , ( 1.5)
where the n and j indices can no longer be omitted since they are not obvious. The difference between these two possible definitions reads

T 00 + j+ σn 2 n -T 00 - j+ σn 2 n = ϕ j n+1 -ϕ j n-1 δ 2 R j+ σn 2 n
, where

R j+ σn 2 n = ϕ j n+1 + ϕ j n-1 1 + δ 2 8 2 + ϕ j n 2 + ϕ j+σn n 2 -ϕ j n -ϕ j+σn n .
On the other hand, the total energy is given by

Q 0 ± n = δ N -1 j=0 T 00 ± j+ σn 2 n .
It can be shown (using periodic boundary conditions) that these two definitions are equivalent,

Q 0 + n = Q 0 -n+1 = E n ,
defining the total energy at time t = n δ / √ 2 with no ambiguity. Now, this energy is exactly conserved if

Q 0 + n = E n = E n-1 = Q 0 -n
, that will be satisfied as soon as

T 00 + j+ σn 2 n = T 00 - j+ σn 2 n , that is to say, if R j+ σn 2 n = 0 .
Since R involves the samples of the field at different times this equation is a sort of equation of motion. Moreover, since it is explicit in ϕ j n+1 , the evolution of a given state can be efficiently followed using

ϕ j n+1 = -ϕ j n-1 + ϕ j n + ϕ j+σn n 1 + δ 2 8 2 + ϕ j n 2 + ϕ j+σn n 2 .
The nodes of the lattice that are involved in this equation are highlighted in fig. 5: the vertex in the left hand side of the equation is represented as and the ones in the right hand side are represented as .

The stress-energy tensor

The 00 component of the stress-energy tensor was defined in eq. (1.5). The two remaining independent components are

T 01 ± j+ σn 2 n = T 10 ± • • = 1 2 ϱ 0 ϕ j+ σn 2 n 2 - 1 2 ϱ 1 ϕ j+ σn 2 n 2 , T 11 ± j+ σn 2 n = 1 2 ϱ 0 ϕ j+ σn 2 n 2 + 1 2 ϱ 1 ϕ j+ σn 2 n 2 + 1 4 - 1 8 1 + ϕ j n±1 2 2 + ϕ j n 2 + ϕ j+σn n 2 .
These two definitions (+ and -) are inequivalent (they only match once integrated over space) but both of them are valid and each one leads to two residues.

Ď∓ 0 T 00 ± -T 10 ± j n + Ď∓ 1 T 00 ± + T 10 ± j n = √ 2 0 ± j n , Ď∓ 0 T 01 ± -T 11 ± j n + Ď∓ 1 T 01 ± + T 11 ± j n = √ 2 1 ± j n ,
where

ϱ 0 T µν ∓ j n = ± T µν ∓ j-σn 2 n±1 -T µν ∓ j± 1 2 n δ , ϱ 1 T µν ∓ j n = ± T µν ∓ j-σn 2 n±1 -T µν ∓ j∓ 1 2 n δ ,
with both µ and ν being either 0 or 1. In practice, the two definitions of the residues behave in the same way and we will only present 0 = 0 + and 1 = 1 + . The same reasoning can be applied to the conservation of the charges and will not be detailed.

Energy conservation

Although the energy is defined without ambiguity (since the two definitions of the stress-energy tensor are equivalent once integrated over space), we still define the total energy at time t = n δ / √ 2 in two ways

E± n = Q 0 ± n ,
and we follow their time evolution independently. The numerical outcome is shown in fig. 6.

We observe, first of all, that the two definitions of the energy behave exactly in the same way, confirming that there is no ambiguity. Then we stress that the value of the energy differs from the exact one (the difference is of order 2% ). This is not surprising and is due to the discretization process (what can be astonishing is that the difference is so small). We also observe that the total energy is exactly conserved, as expected.

At this point, one could reasonably conclude that the BDdV method is a very good choice for the short and long time integration of conservative field theories. However, as we shall see in section 1.5, this conclusion would be premature. Unfortunately, the stress-energy tensor is not conserved locally as it is the total energy.

The msilcc method: a short review of expected properties

The detailed description of the multi-symplectic integrator in light-cone coordinates (msilcc) will be given in the two last sections of the paper.

We just want to stress here that the lattice is the same as the one used in the BDdV method (see fig. 5). In brief, the difference lies in the rules, which are closer to the ones employed in Newton's method.

The msilcc method is designed such that the discretization process exactly preserves the multi-symplectic structure of the phase space. It is also implemented in such a way to respect, as much as possible, the rules of differential calculus. The direct consequence is that the local conservation of the stress-energy tensor is remarkably good even on long time-scales. However, the method is not engineered to conserve the global charges and we do not expect to have the same kind of "magic" compensation of local errors that ensures the BDdV method.

Let us finish this very brief presentation of the expected properties of the msilcc method by showing how the total energy (which is here uniquely defined) behaves in time in fig. 7. We first observe that there seem to be two interlaced curves. Actually, this is not the case, there is only one energy that jumps from one carrier curve to the other. This "double" structure is due to the lattice geometry in combination with the discretization rules. More precisely, when the time index is odd there is a shift of the space index and hence the field is not sampled at the same places, leading to a different energy. Therefore, there are "two curves", one for odd times and the other one for even ones.

Having clarified the effect of the time-discretization we now describe the actual time variation of the total energy. Firstly, over short time-scales, the deviations are around 1% of its value. Secondly, there is no long term trend to increase this deviation. Accordingly, these two remarks allow us to promote the msilcc method as a good candidate for the long time integration of conservative field theories. In the following we analyse the local conservation laws.

Comparison

In this Section we compare the performance of the three numerical integrators discussed so far in a complementary way following what we have already discussed in section 1.1.4.

Influence of the non-linearity

The first situation will explore the influence of the nonlinearities (in coordination with the influence of the quality of the lattice spacing). Figures 8 and9 represent the error committed on the conservation laws by the different methods as a function of A (the amplitude of the initial condition). For this test, the system is integrated up to a time t /L = 1 (i.e. the solution is obtained over a square), and the error, denoted ∆(y = 0), is taken as the largest deviation from the identity y = 0 ever encountered (in absolute value and divided by a characteristic quantity, as introduced in section 1.1.4).

Let us start by discussing fig. 8. We first observe that all methods improve their performance for smaller A.

We emphasize that there are some missing data-points for the Newton method (behind A ∼ 20). This is due to the fact that the approximation becomes unstable before the final integration time for too strong non-linearity. Beyond that point, the solution diverges and the errors as well. Behind this feature there is a first important remark: the largest the effect of the non-linearity, and the worse the quality of the sampling, the quickest the Newton approximation becomes unstable. This fact is worrying since the parameter region we want to explore is precisely the one in which the non-linearities are relevant. Concomitantly, we want to reach long times and it is not desirable to have to oversample the field in time with a too small time spacing.

To pursue the remarks on the Newton method, we stress that it behaves quite well while it remains stable (local errors are between 10 -2 and 10 -1 ). However, it is at minimum 3 orders of magnitude worse than the msilcc method. The difference can be attributed to the fact the Newton's method is explicit, while the msilcc one is implicit (explicit schemes are known to be worse than implicit ones).

Concerning the BDdV method, the violation of the local conservation laws is very important with an error that ranges between 10 -1 and 10 +1 . Quite surprisingly, the exact conservation of the energy is only due to the compensation of these large errors once integrated over space.

Finally, the msilcc method produces errors that range from 10 -9 to 10 -2 . They disappear very abruptly when the non-linearity becomes negligible. This is actually due to the fact that the method is exact for a linear problem (as we will show in section 3.4.1). So, the msilcc method appears, for now, as a very good choice to integrate conservative field theories over long times. Let us now look at how the errors on the charges behave and, in particular, the energy one (see fig. 9). First of all, we remark that the errors committed on the conservation of the energy are in agreement with what we observed earlier when we showed their evolution in time. The figure shows that the BDdV prescription is better (actually almost exact since 10 -14 is of the order of the machine precision) than the msilcc method which is itself better than the Newton method. The conservation of the second charge is almost exact for both the BDdV and the msilcc schemes and is of the order of the conservation of the energy for the Newton method.

In conclusion, Newton's method presents a not so bad local conservation of the stress-energy tensor. However, by accumulation of these errors, the charges are not conserved and the violation of their conservation, though not too large, is not sufficiently good for high precision measurements. Concerning the BDdV method, it presents very poor local conservation properties that, quite surprisingly, lead to a very good conservation of the charges (due to a deceptive cancelation of the errors). Finally, the msilcc method behaves more like Newton's but with much better local conservation features that yield an acceptable conservation of the charges.

10 0 A ∆(∂ 0 Q ν = 0) msilcc BDdV Newton max t,x ∆ ∂ 0 Q 1 = 0 max t,x ∆ ∂ 0 Q 0 = 0

Long time behavior

We now explore the long time properties. Figures 10 and 11 represent the error committed on the conservation laws as a function of time. Note that for each method the integration is performed using A = 10 and two errors are displayed: the first one is the largest deviation ever encountered (in absolute value), and the second one is the largest deviation at time t. The comments made on figs. 8 and 9 still hold and we will only describe the time behavior here.

Firstly, we observe that the Newton method rapidly becomes unstable (after t /L = 4) and is no longer able to describe the evolution of the field.

Secondly, we remark that the instantaneous error evolves in time (over several orders of magnitude). So, it is preferable to consider, instead of the instantaneous error, the worst one ever encountered from the beginning.

Finally, the most interesting comment that applies to the BDdV and the msilcc method as well is that the worse errors occur during the short time behavior: after a rapid evolution (of the order of the characteristic time-scale of the system, as we observe by comparing figs. 8 and 9 with fig. 1 for A = 10), the error stabilizes to a value which can be hopefully considered as definitive (of the order of 10 0 for the BDdV method and 10 -5 for the msilcc one). This is particularly true for the local conservation of the stress-energy tensor but less clear for the conservation of the charges even though they seem to reach a constant too.

Preliminary conclusion

A first element of conclusion is that we need to be extremely wary of methods that possess remarkable properties on some observables but not necessarily the most fundamental objects of the theory.

Within the three methods here presented, the msilcc is the only one that one could trust to integrate a conservative field theory over a long time interval. However, its implementation has a cost: the discrete equations of motion are implicit and more expensive to solve (in terms of computational time) than the other two methods. Fortunately, the scheme remains well-defined locally (i.e. there is no need to solve the set of algebraic equations globally) and it can be easily scaled to larger volumes and/or extended to theories defined on higher dimensions.

Up to now, we have eluded the concrete results in term of the (numerical) solution of the pde and one can imagine that all these elements of conservation only have a negligible influence. This is not true. As an example, over an integration time as short as t /L = 1, we observe differences of the order of 1% between the solutions obtained with the different methods (under the same conditions as described in section 1.1.2 and for A = 10). In some situations the differences can become dramatically larger (up to 20%) as we show on fig.

12 that represents the field after an integration time t /L = 1, still with periodic boundary conditions, but with an initial state that corresponds to a particular solution of eq. (1.1) in term of a Jacobi elliptic function [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF]. in terms of Jacobi elliptic functions after a time t /L = 1 (a particular solution for different initial condition from the ones used so far). The datapoints represent the field evolved by the different numerical methods (triangles for Newton, diamonds for BDdV and circles for msilcc). On this scale we see no difference between the exact and the numerical solutions obtained with the Newton and msilcc methods.

Symmetry breaking potential

So far we have not considered the influence of r (the parameter that accompanies the quadratic term in the potential) and we want to show that the msilcc method behaves just as well for a potential in double well. Figure 13 shows the error committed on the conservations of the charges and the stress-energy tensor as a function of r. First of all, we observe that the msilcc method has the same conservation properties whatever the shape of the potential. Secondly, we remark that the errors do not depend on r (except in the very large r limit). This is a direct consequence of a feature that will be proved in section 3.4.1: the deviations from the conservation of the stress-energy tensor only arise with the non-linear part of the Hamiltonian. Actually we no longer observe this feature in the large r limit since for A = 10 and for such values of r the field can only "explore" the quadratic part of the potential. In other words, when |r| increases the non-linear part of the potential disappears and the errors are modified. Now, the rest of this paper will be devoted to the construction of the msilcc method in a more general setup

∆(∂ µ T µν = 0), ∆(∂ 0 Q ν = 0) -0.1 0 0.1 max t,x ∆ ∂ 0 Q 1 = 0 max t,x ∆ ∂ 0 Q 0 = 0 max t,x ∆ ∂µT µ1 = 0 max t,x ∆ ∂µT µ0 = 0 0.1 1 10 100 r
Figure 13: Error committed by the msilcc method on the conservations of the charges and stress-energy tensor as a function of the parameter that accompanies the quadratic term in the potential: r. The initial condition follows eq. ( 1.3) with A = 10 and the errors are accumulated over an integration time of t /L = 1 while we still have L / √ 2 δ = 128. The vertical axis is cut between 10 -5 and 10 -15 while the horizontal axis is in logarithmic scale from -100 to -0.1 and from 0.1 to 100 (scale is linear between -0.1 and 0.1). Positive values of r mean that the potential has only one minimum at φ = 0 while negative values of r mean that the potential is a double well with two minimum at φ = ± √ -r (so far r was setted to 1).

which will be illustrated through the example of the nonlinear wave equation. In the next section we will introduce a number of necessary concepts and we will discuss how the method is concretely applied, while in the last section we will finally present its construction.

Applicability and preliminaries

This section introduces the concept of symplecticity (for mechanical systems) and its generalization to multisymplecticity (for field theories). We then introduce the idea of Hamiltonian pdes through the De Donder -Weyl formalism. Afterwards, we show how to rewrite the equation of motion (when possible) in a way that allows one to use the msilcc method. Finally, we will exhibit (both in a general setup and with an example) the multisymplectic structure, prove its conservation, and discuss the definition and properties of the stress-energy tensor.

Most of the points presented in this section are just reminders except for two of them. Firstly, the link between the De Donder -Weyl formulation and multisymplectic geometry is not so common (usually, the De Donder -Weyl formulation is treated through the formalism of the poly-symplectic geometry). Secondly, as far as we know, the discussion of the degeneracy of the multi-symplectic structure, and especially its resolution, is completely new.

2.1 From symplectic to multi-symplectic structure

Emergence of the symplectic structure

Let us consider a mechanical problem, constituted by N particles, in its Hamiltonian formulation (for an exhaustive presentation of the Hamiltonian formalism see, for example, [START_REF] José | Classical Dynamics: A Contemporary Approach[END_REF]). The i th body is described by its position and momenta, or Darboux coordinates, q i , p i , and the dynamics are fully characterized by the Hamiltonian, H q i , p i . The phase space, Ω, is a differential manifold locally parametrized by the union of the Darboux coordinates.

The symplectic character of the phase space of a mechanical system is rather well known, so we will only remind its central aspects [START_REF] José | Classical Dynamics: A Contemporary Approach[END_REF]. The symplectic manifolds are differentiable manifolds equipped with a closed non-degenerate 2-form. We now recall how this 2-form is constructed.

Let us describe Ω using the generalized coordinates {ζ a } ≡ q i , p i with a ∈ 1, 2N (the union of the positions and momenta). We define ∂ a = ∂ /∂ζ a and {d a = dζ a }, respectively, as a basis for T Ω, the tangent space of Ω, and T * Ω, its dual (see [START_REF] Nakahara | Geometry, topology, and physics[END_REF] for the necessary concepts of geometry). Now, we provide to Ω a 2-form

ω = ω ab d a ∧ d b ,
where we have used the Einstein summation convention, as we will do (unless explicitly stated) in the rest of this paper. In particular, for a mechanical system and in Darboux coordinates this 2-form reads

ω = dq i ∧ dp i .
It is obviously closed (dω = 0), and it can be shown to be non-degenerate (det ω = 0, the proof will be given later). Therefore, (Ω, ω) defines a symplectic manifold where ω is the symplectic structure.

The conservation of the symplectic structure under a Hamiltonian flow leads to Liouville's theorem (since the volume form on Ω is obtained from ω).

The Poisson bracket is directly related to ω since for any pair of differentiable functions of the phase space variables, * ({ζ a }) and ({ζ a }), we have

{ * , } = ω(χ * , χ ) = ω [ab] ∂ a * ∂ b , (2.1)
where χ * is the vector field associated to * by ω defined as

ω(χ * , •) = d * [•] , (2.2)
and where [ab] denotes the antisymmetric combination of the indices

ω [ab] = ω ab -ω ba 2 . ( 2
{ * , } = ∂ * ∂p i ∂ ∂q i - ∂ * ∂q i ∂ ∂p i .
Symplecticity is an essential structure of the phase space of mechanical systems and numerical methods that preserve it are very popular [START_REF] Mclachlan | Geometric integrators for ODEs[END_REF]. For field theories, the symplectic structure needs to be generalized if we want to develop the same class of methods for pdes.

Generalization to the multi-symplectic structure

A multi-symplectic manifold, (Ω, {ω µ }), is a differential manifold equipped with several independent symplectic structures 4 (in our case, one per space-time direction). As we shall see later, this structure is natural for a field theory subject to a conservation law that fits the invariance under the Hamiltonian flow.

The next section will illustrate how the multi-symplectic aspect of the phase space of a field theory emerges. We will first recall how to construct a Hamiltonian formulation of a field theory. Then, we will identify the phase space, highlight the underlying multi-symplectic structure, and establish its conservation. Finally, we will discuss the properties of the stress-energy tensor.

2.2

The De Donder -Weyl (DW) Hamiltonian formulation of field theories

From Lagrangian to DW Hamiltonian formulation

Let us start with a space-time Lorentzian manifold, M, of dimension D = 1 + d. We assume M to be non dynamic (the metric is not subject to an equation of motion). Without lose of generality, M will be supposed to be flat with metric η ≡ diag (1, -1, • • • , -1). We define a local coordinate system {x µ }, with ∂ µ = ∂ /∂x µ a basis of T M, and µ ∈ 0, d .

Next, we consider a field theory on M described by the action S φ i , where φ i with i ∈ 1, N is a collection 4 Let us complete this definition with some vocabulary remarks.

One can also encounter in the literature the concept of polysymplectic (or n-plectic) manifold which is actually different: a n-plectic manifold is a differential manifold with a closed, nondegenerate, (n + 1)-form (poly-symplectic is for any n > 1 while symplectic stands for 1-plectic). So, a multi-symplectic manifold is poly-symplectic too since we can define = µ ω µ , but the reverse is not necessarily true. Finally, we want to stress that the meaning of poly-symplectic and multi-symplectic can be exchanged depending on authors.

of dynamic fields. This action originates in a Lagrangian density L (which is assumed to depend only on the field and its first derivatives) and reads

S = d D x L φ i , ∂ µ φ i ,
where d D x is the measure over M and µ ∈ 0, d . The stationarity of S leads to the Euler -Lagrange equations (i.e. the equations of motion) for the fields

δS δφ i = 0 = ∂L ∂φ i -∂ µ ∂L ∂(∂ µ φ i ) .
This is the Lagrangian formulation of a field theory.

The idea of the Hamiltonian formulation of classical mechanics is to substitute the generalized speed ( qi ) by a conjugate momentum (p i ). Obviously the same reasoning can be applied to a field theory but, unfortunately, it breaks the Lorentz covariance of the theory. The idea of De Donder and Weyl is to reestablish the covariance by introducing one conjugate momentum per direction of space-time (such that they are treated on an equal footing). Thus, we define

ψ i µ = ∂L ∂(∂ µ φ i ) ,
as the conjugate momentum of the field φ i along the µ th direction of M. Then, provided that the following Legendre transform is not singular,

H = ψ i µ ∂ µ φ i -L ,
defines the De Donder -Weyl Hamiltonian density. Henceforth, the unknowns are the fields ( φ i ) with their conjugate momenta in each direction of space-time ({ψ i µ }). Together, they are the local coordinates of a differential manifold, Ω, which is the De Donder -Weyl definition of phase space (a multi-symplectic manifold as we will prove later). The Hamilton equations generalize to

∂ µ ψ i µ = - ∂H ∂φ i , (2.4a)
∂ µ φ i = ∂H ∂ψ i µ . (2.4b)
This is the DW Hamiltonian formulation of a field theory (see [START_REF] Kanatchikov | On the canonical structure of the De Donder-Weyl covariant Hamiltonian formulation of field theory I. Graded Poisson brackets and equations of motion[END_REF][START_REF] Paufler | De Donder-Weyl equations and multisymplectic geometry[END_REF][START_REF] Rieth | The Hamilton-Jacobi theory of De Donder and Weyl applied to some relativistic field theories[END_REF] for a detailed review of this formulation of classical field theories).

Equations (2.4a) and (2.4b) can be rewritten in a more symmetrical way (as in the case of classical mechanics) by introducing a Poisson bracket for each space-time direction. This collection of Poisson brackets is intimately related to the multi-symplectic structure and each bracket is associated to a symplectic 2-form.

Let us now describe Ω in the generalized coordinates {ζ a } ≡ φ i ∪ {ψ i µ } where a is an index conveniently chosen to sweep the collection. Following [START_REF] Bridges | Multi-symplectic structures and wave propagation[END_REF][START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity[END_REF][START_REF]Numerical methods for Hamiltonian PDEs[END_REF], we define on T Ω the vector ζ = ζ a ∂ a that allows one to rewrite eqs. (2.4a) and (2.4b) in the abstract form

M µ • ∂ µ ζ = ∇H , (2.5)
where {M µ } is a set of constant skew-symmetric matrix of T Ω. The gradient is defined as

∇ : C ∞ (Ω) → T Ω ∇H → (dH) * = δ ab ∂ b H ∂ a ,
and we recall that {∂ a } is a basis of T Ω.

The multi-symplectic structure will directly arise from eq. (2.5). We will come back to this point later but we already emphasize that, as long as a pde can be written in such a form, the msilcc method will apply.

The next step is to exhibit the multi-symplectic structure and to prove its conservation. This will be done in the next section, but let us first clarify the procedure presented above on an example.

The non-linear wave equation

We consider the dynamics of a real scalar field, φ, whose Lagrangian density is given by

L = 1 2 ∂ µ φ ∂ µ φ -V (φ) .
The Euler -Lagrange equation reads

∂ µ ∂ µ φ + V (φ) = φ + V (φ) = 0 .
We introduce now

ψ µ = ∂L ∂(∂ µ φ) = ∂ µ φ ,
which is the conjugate momentum of φ in the µ th direction. Then, the DW Hamiltonian density reads

H = ψ µ ∂ µ φ -L = 1 2 ψ µ ψ µ + V (φ) . Defining ζ T = φ ψ 0 ψ 1 • • • ψ d ,
the equation of motion is given by eq. (2.5) provided that

M µ ab = δ µ a-1 δ 0 b -δ µ b-1 δ 0 a ,
where a, b ∈ 0, D .

In dimension D = 1 + 0 (d = 0), aliasing q ≡ φ, p ≡ ψ 0 and H ≡ H, we recover the expected Hamilton equation of motion for a mechanical problem

0 -1 1 0 • q ṗ =     ∂H ∂q ∂H ∂p     = V (q) p .
In dimension D = 1 + 1 (d = 1) the M matrices read

M 0 =   0 -1 0 1 0 0 0 0 0   , M 1 =   0 0 -1 0 0 0 1 0 0   ,
and we stress that, for both of them, the eigenvalues are 0 and ±i.

In the general case all the M matrices have the same eigenvalues: ±i and 0 (d times degenerate). As we will see in the next section this fact causes some difficulties.

For the moment, we are going to highlight the multisymplectic structure of Ω relying on eq. (2.5). Then, we will return to the vanishing eigenvalues and we will explain how to treat them in the particular case of the non-linear wave equation.

The multi-symplectic structure

This section will be devoted to the construction of the multi-symplectic structure on Ω and, again following [START_REF] Bridges | Multi-symplectic structures and wave propagation[END_REF][START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity[END_REF][START_REF]Numerical methods for Hamiltonian PDEs[END_REF], we are going to show how it directly emerges from eq. (2.5).

In the previous section we have shown that for a single scalar field theory the Hamilton equation is fully characterized by a Hamiltonian density and a collection of D (constant and skew-symmetric) matrices ({M µ }) that can be used to define D 2-forms {ω µ }:

ω µ = - 1 2 M µ ab d a ∧ d b .
The action of these 2-forms on a pair of vectors is

ω µ ( * , ) = M µ • * , = - * , M µ • ,
for any * , ∈ T Ω, where ␣, ␣ is the scalar product on T Ω.

All the {ω µ } are closed (dω µ = 0) since all the {M µ } are independent of the fields. On the other hand, let us assume for the moment (we will come back to this at the end of this section) that they are all non-degenerate. Finally, since the {M µ } are linearly independent (as we will shortly show) the {ω µ } are linearly independent as well. So, (Ω, {ω µ }) is a multi-symplectic manifold.

Let us first prove the independence of the {M µ } matrices. In order to lighten the following computation, let us introduce space-time indices which behave differently under the Einstein summation rule. From now, the index takes just one value and does not imply summation (even if repeated), while the σ index (σ = ) behaves in the standard way. All other indices are unaffected. Let us now start by supposing that one of the {M µ } is linearly dependent on the others: M = α σ M σ . Then, in eq. (2.5), the operator on the left hand side becomes

M µ • ∂ µ = M • ∂ + M σ • ∂ σ = α σ M σ • ∂ + M σ • ∂ σ = M σ • (α σ ∂ + ∂ σ ) = M σ • ∂σ ,
where the ˜refers to a different coordinate system. As the direction ∂ has just disappeared from the differential operator, it means that the dynamics along this direction are trivial. Thus the {M µ } are linearly independent.

To complete the proof of the multi-symplecticity of Ω it remains to discuss the question of the degeneracy of the {ω µ }. The 2-forms will be non-degenerate as long as they all satisfy det ω µ = det M µ = 0 .

In other words the {M µ } matrices should not have any zero eigenvalue. Nevertheless we have seen in the previous section that for the example of the non-linear wave equation all the {M µ } matrices have the same eigenvalues and especially a d times degenerate zero. So, as it is, the phase space of the non-linear wave equation is multi-symplectic only for d = 0 (that is to say, for the mechanical problem).

We are now going to present the resolution of this problem on the particular example of the non-linear wave equation.

The same reasoning can be applied to other theories.

Let us first recall that the aim of the DW Hamiltonian formalism was to preserve covariance. Thus, it treats all the space-time directions on an equal footing by introducing a conjugate momentum for each. Still, time is not space, and the conjugate momentum in time will have a different status. We will call it the canonical one.

Degeneracy comes from the existence of non-canonical conjugate fields. So, the idea is to ensure that all the conjugate momenta be canonical (i.e. the conjugate momentum along time of a dynamic field). The solution is to add extra fields, interacting with each other, such that the new field theory allows a conjugate momentum to be shared by several fields and thus to be canonical for someone. Finally, we still want to preserve covariance and, for each direction of space-time, each momentum needs to be the conjugate of a field.

Putting these elements together, we can modify the Lagrangian density such that it describes an equivalent, but non-degenerate, problem. This construction can be graphically represented by placing both the fields and their conjugate momenta on the vertices of a D-dimensional hypercube: each direction stands for a space-time direction and a line means that one of the fields on the edges is the conjugate momentum (along this direction) of the other. So, on any path, a dynamic field alternates with a conjugate field. We have drawn it for the three first dimensions in fig. 14. Now, to construct the new field theory explicitly, we first consider the collection of fields

∂ 0 d = 0 φ ψ 0 ∂ 0 ∂ 1 φ ψ 0 ψ 1 γ d = 1 ∂ 0 ∂ 1 ∂ 2 φ ψ 0 ψ 1 ψ 2 γ 01 γ 02 γ 12 Γ d = 2
Φ (i) µ 1 •••µ i = Φ (i) [µ 1 •••µ i ] µ 1 ,••• ,µ i ∈ 0,d i∈ 0,D , (2.6)
where the pair of square brackets denote the antisymmetrization defined as

[µ 1 •••µ i ] = 1 i! p∈P i ({µ 1 •••µ i }) σ(p) p ,
where P i is the symmetric group of i symbols and σ(p) is the signature of p. In particular, eq. ( 2.3) gives this definition for two indices.

Before going further, we stress that the Φ (i) will not be treated as a tensor field but as a collection of scalar fields, conveniently assembled in the same object. The collection (2.6) is separated in two (equal) parts: {Φ (2i) } contains the dynamical fields, while {Φ (2i+1) } contains the conjugate fields. As a final remark, this collection is composed of

D i=0 D i = 2 D
elements and will indeed be suitable to populate the vertices of the D-dimensional hypercube introduced earlier as a graphic representation of this construction.

Next, we consider the Lagrangian density

L 0 = 1 2 ≤ D /2 j=0 1 (2j)! ∂ ν 0 Φ (2j) ν 1 •••ν 2j ∂ ν 0 Φ ν 1 •••ν 2j (2j) - ≤ D /2 j=0 1 (2j)! V [ν 1 •••ν 2j ] (2j) Φ (2j) ν 1 •••ν 2j ,
where the ν indices are in 0, d , which is nothing else than the concatenation of 2 d independent theories 5 (note that all these theories need to belong in the class of the nonlinear wave equation but it is not required that they have the same potential). Then we add to this Lagrangian two vanishing coupling terms. First, we obviously have

0 = + 1 2 ≤ D /2 j=1 1 (2j -1)! Φ (2j) µν 1 •••ν 2j-1 ∂ µ ∂ ν Φ νν 1 •••ν 2j-1 (2j) - 1 2
≤ D /2 j=1 1 (2j -1)! Φ (2j) µν 1 •••ν 2j-1 ∂ µ ∂ ν Φ νν 1 •••ν 2j-1 (2j)
, and then, using an integration by parts (and ignoring the boundary terms), we claim that

L + = - 1 2 ≤ D /2 j=1 1 (2j -1)! ∂ µ Φ νν 1 •••ν 2j-1 (2j) ∂ ν Φ (2j) µν 1 •••ν 2j-1 + 1 2 ≤ D /2 j=1 1 (2j -1)! ∂ ν Φ νν 1 •••ν 2j-1 (2j) ∂ µ Φ (2j) µν 1 •••ν 2j-1
will not modify the equations of motion since it is just 0 rewritten in a convenient way. Secondly, we have

0 = - ≤ D /2-1 j=0 1 (2j)! Φ (2j) ν 1 •••ν 2j ∂ µ ∂ ν Φ µνν 1 •••ν 2j (2j+2)
, since it contains the contraction of a symmetric tensor (∂ µ ∂ ν ) with an anti-symmetric one (Φ µν••• ). Again, using an integration by parts (and forgetting the boundary terms), we define

L -= ≤ D /2-1 j=0 1 (2j)! ∂ µ Φ µνν 1 •••ν 2j (2j+2) ∂ ν Φ (2j) ν 1 •••ν 2j
and this term does not affect the dynamics either.

We now consider the theory described by the Lagrangian density

L = L 0 + L + + L -,
which should be equivalent to the simultaneous treatment of 2 d independent problems (that belong in the class of the non-linear wave equation). Let us now construct the DW Hamiltonian formulation of this theory. We first introduce, for each dynamic field in the collection (2.6) and for each direction of space-time, a conjugate momentum

Ψ µ 0 •••µ 2i (2i+1) = Ψ µ 0 [µ 1 •••µ 2i ] (2i+1) = ∂L ∂ ∂ µ 0 Φ (2i) µ 1 •••µ 2i = ∂ µ 0 Φ µ 1 •••µ 2i (2i) + ∂ µ Φ µµ 0 •••µ 2i (2i+2) + 2i η µ 0 [µ 1 ∂ µ Φ |µ|µ 2 •••µ 2i ] (2i) -2i ∂ [µ 1 Φ |µ 0 |µ 2 •••µ 2i ] (2i) + 2i(2i -1) η µ 0 [µ 1 ∂ µ 2 Φ µ 3 •••µ 2i ] (2i-2)
, 5 The anti-symmetrization of the indices in the potential are required since the fields are anti-symmetric objects and in particular V 01

(2) Φ

(2)

01

should coincide with V 10

(2) Φ

(2)

10

.

where we have used that

∂ ∂ ν 0 Φ (j) ν 1 •••ν j ∂ ∂ µ 0 Φ (i) µ 1 •••µ i = j! δ j i δ µ 0 ν 0 δ µ 1 [ν 1 • • • δ µ i ν j ] .
Then, defining

Φ µ 0 •••µ 2i (2i+1) = Ψ [µ 0 •••µ 2i ] (2i+1) = (2i + 1) ∂ [µ 0 Φ µ 1 •••µ 2i ] (2i) + ∂ µ Φ µµ 0 •••µ 2i (2i+2)
, we can prove 6 that

Ψ µ 0 [µ 1 •••µ 2i ] (2i+1) = Φ µ 0 •••µ 2i (2i+1) + 2i η µ 0 [µ 1 Φ µ 2 •••µ 2i ] (2i-1)
, and in consequence, that all the conjugate momenta are indeed already defined in the collection (2.6).

The Euler -Lagrange equation reads

∂L ∂Φ (2i) µ 1 •••µ 2i = ∂ µ 0 ∂L ∂ ∂ µ 0 Φ (2i) µ 1 •••µ 2i = ∂ µ 0 Ψ µ 0 [µ 1 •••µ 2i ] (2i+1) = ∂ µ ∂ µ Φ µ 1 •••µ 2i (2i) + ∂ µ ∂ ν Φ νµµ 1 •••µ 2i (2i+2) + 2i ∂ [µ 1 ∂ µ Φ |µ|µ 2 •••µ 2i ] (2i) -2i ∂ [µ 1 ∂ µ Φ |µ|µ 2 •••µ 2i ] (2i) + 2i(2i -1) ∂ [µ 1 ∂ µ 2 Φ µ 3 •••µ 2i ] (2i-2) = ∂ µ ∂ µ Φ µ 1 •••µ 2i (2i) = Φ µ 1 •••µ 2i (2i)
. Summarizing, we first consider the concatenation of 2 d independent theories (through the Lagrangian density L 0 ).

Based on what we introduced in the previous sections, 6 By substituting the expression of

Φ µ 0 •••µ 2i (2i+1) and Φ µ 2 •••µ 2i (2i-1)
and after a straightforward identification it remains to prove that

(2i + 1) ∂ [µ 0 Φ µ 1 •••µ 2i ] (2i) = ∂ µ 0 Φ µ 1 •••µ 2i (2i) -2i ∂ [µ 1 Φ |µ 0 |µ 2 •••µ 2i ] (2i)
.

On the other hand, we successively have

(2i + 1) ∂ [µ 0 Φ µ 1 •••µ 2i ] (2i) = = (2i + 1) 1 2i + 1 2i+1 j=1 ∂ [µj Φ µ j+1 •••µ 2i |µ 0 |µ 1 •••µ j-1] (2i) = ∂ µ 0 Φ µ 1 •••µ 2i (2i) + 2i j=1 ∂ [µj Φ µ j+1 •••µ 2i |µ 0 |µ 1 •••µ j-1 ] (2i) = ∂ µ 0 Φ µ 1 •••µ 2i (2i) + 2i j=1 ∂ [µj Φ |µ 0 |µ 1 •••µ j-1 µ j+1 •••µ 2i] (2i) = ∂ µ 0 Φ µ 1 •••µ 2i (2i) - 2i j=1 ∂ [µ1 Φ |µ 0 |µ j µ 2 •••µ j-1 µ j+1 •••µ 2i] (2i) = ∂ µ 0 Φ µ 1 •••µ 2i (2i) - 2i j=1 (-1) j-2 ∂ [µ1 (-1) j-2 Φ |µ 0 |µ 2 •••µ j-1 µ j µ j+1 •••µ 2i] (2i) = ∂ µ 0 Φ µ 1 •••µ 2i (2i) - 2i j=1 ∂ [µ 1 Φ |µ 0 |µ 2 •••µ 2i ] (2i) = ∂ µ 0 Φ µ 1 •••µ 2i (2i) -2i ∂ [µ 1 Φ |µ 0 |µ 2 •••µ 2i ] (2i)
.

we know that the DW Hamiltonian formulation of this theory requires the collection {Φ

(2i) , Ψ (2i+1) } of 2 d (D + 1)
fields to construct the phase space (2 d dynamical fields and D 2 d conjugate momenta). We also know that this phase space is not a multi-symplectic manifold since the multi-symplectic structure of each "sub"-theory is, independently, degenerate. As we claimed earlier, this comes from the existence of a non-canonical conjugate field. To fix this issue, we introduce the additional coupling terms L + and L -. Above all, these couplings do not affect the dynamics and each "sub"-theory remains independent. However, these couplings have an interesting side effect: they allow a conjugate momentum to be shared by several dynamical fields and thus enable all the conjugate fields to be canonical. Thus, correctly chosen, L + and L -lead to the closed collection

{Φ (2i) , Φ (2i+1) } = {Φ (i) } of 2 D fields,
equality composed of dynamic fields and canonical conjugate momenta, and where "closed" has two meanings. On the one hand, all the dynamic fields have all their conjugate momenta in the collection. On the other hand, all the conjugate fields are, for all the directions of space-time, the conjugate momentum of a dynamic field that belong in the collection. Now, from the Euler -Lagrange equations and the definitions of the conjugate momenta, we get the De Donder -Weyl -Hamilton equations

-∂ µ Φ µµ 1 •••µ 2i (2i+1) -2i ∂ [µ 1 Φ µ 2 •••µ 2i ] (2i-1) = ∂H ∂Φ (2i) µ 1 •••µ 2i (2.7a) (2i + 1) ∂ [µ 0 Φ (2i) µ 1 •••µ 2i ] + ∂ µ Φ (2i+2) µµ 0 •••µ 2i = ∂H ∂Φ µ 0 •••µ 2i (2i+1) (2.7b) where H = 1 2 ≤ (D-1) /2 j=0 1 (2j + 1)! Φ ν 0 •••ν 2j (2j+1) Φ (2j+1) ν 0 •••ν 2j + ≤ D /2 j=0 1 (2j)! V [ν 1 •••ν 2j ] (2j) Φ (2j) ν 1 •••ν 2j ,
is the Hamiltonian density associated to L.

In order to identify the form of eq. (2.5), we need to flatten all these indices (i.e. for every configuration of values of all these indices we associate one, and only one, index). So, if the list (µ 1 • • • µ i ) is sorted and free of duplicates, we define

fl(i, µ 1 • • • µ i ) = 1 + i-1 j=0 D j + i j=1 µ j -1 ν j = µ j-1 +1 d-i+j+1 ν j+1 = ν j +1 • • • d ν i = ν i-1 +1 1 ∈ 1, 2 D ,
with µ 0 = -1 and where the stacking of sums can be re-expressed in term of generalized harmonic numbers as well.

In the following, we assume that both (α 1 • • • α i ) and (β 1 • • • β j ) are sorted and duplicate free. Thus, defining the vector state ζ as

ζ a = ζ a=fl(2i,α 1 •••α 2i ) = Φ (2i) α 1 •••α 2i = ζ a=fl(2i+1,α 1 •••α 2i+1 ) = Φ α 1 •••α 2i+1 (2i+1)
, eqs. (2.7a) and (2.7b) can be rewritten in the form of eq. (2.5), i.e.

M µ a b ∂ µ ζ b = ∂ a H , provided that M µ a=fl(i,α 1 •••α i ) b=fl(j,β 1 •••β j ) = (i mod 2 -1)× δ i+1 j j k=1 (-1) k-1 δ α 1 β 1 • • • δ α k-1 β k-1 δ µ β k δ α k β k+1 • • • δ α i β j + δ i-1 j i k=1 (-1) k-1 δ α 1 β 1 • • • δ α k-1 β k-1 η µα k δ α k+1 β k • • • δ α i β j + (i mod 2)× δ j+1 i i k=1 (-1) k-1 δ β 1 α 1 • • • δ β k-1 α k-1 δ µ α k δ β k α k+1 • • • δ β j α i + δ j-1 i j k=1 (-1) k-1 δ β 1 α 1 • • • δ β k-1 α k-1 η µβ k δ β k+1 α k • • • δ β j α i .
These {M µ } matrices are skew-symmetric, linearly independent and non-degenerate (they all have two eigenvalues, ±i, 2 d times degenerate). The phase space of the theory, Ω, is now a multi-symplectic manifold and we have finally obtained a correct covariant Hamiltonian formulation of the non-linear wave equation.

Let us illustrate how this construction works in the particular dimension D = 1 + 1. We start by considering the collection (2.6) and for notational convenience we create aliases for these 2 2 fields as

Φ (0) = φ (2.8a) Ψ 0 (1) = Φ 0 (1) = ψ 0 (2.8b) Ψ 1 (1) = Φ 1 (1) = ψ 1 (2.8c) Φ (2) 01 = -Φ (2)
10 = γ (2.8d)

Ψ 001 (3) = -Ψ 010 (3) = Φ 1 (1) = ψ 1 Ψ 101 (3) = -Ψ 110 (3) = Φ 0 (1) = ψ 0 .
Here, we have introduced one new dynamic field, γ, independent of φ. See fig. 14 for a graphical representation of these fields. Then we consider the concatenation of these two theories described by the Lagrangian density

L = 1 2 (∂ 0 φ) 2 - 1 2 (∂ 1 φ) 2 - 1 2 (∂ 0 γ) 2 + 1 2 (∂ 1 γ) 2 + ∂ 0 φ ∂ 1 γ -∂ 0 γ ∂ 1 φ -V (φ) -Ṽ (γ) ,
where Ṽ is the potential of the extra theory which can be freely chosen. The (canonical) conjugate momenta are defined by

ψ 0 = ∂L ∂(∂ 0 φ) = ∂L ∂(∂ 1 γ) = ∂ 0 φ + ∂ 1 γ , ψ 1 = ∂L ∂(∂ 1 φ) = ∂L ∂(∂ 0 γ) = -∂ 0 γ -∂ 1 φ .
The Hamiltonian density is

H = 1 2 ψ 0 2 - 1 2 ψ 1 2 + V (φ) + Ṽ (γ) ,
and defining the vector state

ζ T = φ ψ 0 ψ 1 γ ,
the dynamics are fully described by eq. (2.5) provided that

M 0 =     0 -1 0 0 1 0 0 0 0 0 0 1 0 0 -1 0     , M 1 =     0 0 -1 0 0 0 0 1 1 0 0 0 0 -1 0 0     .
These two matrices are no longer degenerate (their eigenvalues are +i, +i, -i and -i) and they define a multisymplectic structure on the phase space.

Recaping, we started with two theories (one for φ and one for γ). While isolated they break the multi-symplectic structure of phase space. Joined together they complete each other such that the theory of the coalition restores multi-symplecticity.

Using this construction, we reduced the number of unknowns for each independent "sub"-theory from 3: φ, ψ 0 , ψ 1 and γ, ψ 1 , ψ 0 to 2: φ, ψ 0 and γ, ψ 1 . In a general settup we reduced the number of unknowns per independent "sub"-theory from D + 1 (the field and all its conjugate momentums) to 2 (the field and its canonical conjugate momentum) as in the standard non-covariant Hamiltonian formulation of field theory.

Finally, we want to stress that from the numerical point of view, to integrate a D-dimensional theory, we actually need to integrate 2 d D-dimensional theories. These "extra" theories can be used in two ways (or a mix of the two):

i. By considering a theory the solution of which is known we get an error estimate of the integration process. This feature comes from the fact that the integration is performed through the conjugate momentum which is itself shared between different dynamic fields. If an error occurs during the integration of one of the fields, it will reverberate on the others and will be caught by the control field(s).

ii. They can be used to integrate, at the same time, several replicas of the theory (in a statistical approach for example) or even different theories.

In the present section we have introduced a construction that leads to a phase space that is a multi-symplectic manifold. In the next sections, we will first prove the conservation of the multi-symplectic structure under the Hamiltonian flow. Then we will define the stress-energy tensor, the charges and discuss their properties.

Conservation of the multi-symplectic structure

To prove the conservation of the multi-symplectic structure, we consider the dual of eq. (2.5) which reads

M µ ab ∂ µ ζ b d a = ∂ a H d a i.e. ω µ (∂ µ ζ, •) = dH .
(2.9)

Now, taking the exterior derivative of it, we successively have

d(ω µ (∂ µ ζ, •)) = ddH = 0 = d M µ ab ∂ µ ζ b d a = M µ ab d ∂ µ ζ b ∧ d a = M µ ab ∂ µ d b ∧ d a ,
where we have used that d and ∂ µ commute since they act in different spaces (Ω does not depend on the position on M).

On the other hand, we have

∂ µ ω µ = 1 2 M µ ab -∂ µ d a ∧ d b -d a ∧ ∂ µ d b = 1 2 M µ ab -∂ µ d a ∧ d b + ∂ µ d b ∧ d a = M µ ab ∂ µ d b ∧ d a = 0 .
Thus, we have proved (on-shell) the local conservation of the multi-symplectic structure

∂ µ ω µ = 0 . (2.10)
By definition a multi-symplectic integrator is a numerical method that exactly preserves the discrete version of eq. (2.10).

The stress-energy tensor, its conservation and the charges

We define now the stress-energy tensor as the symmetric 2-tensor

T µν = 1 2 (ω µ (∂ ν ζ, ζ) + ω ν (∂ µ ζ, ζ) + η µν ω κ (ζ, ∂ κ ζ)) + η µν H (2.11) = T νµ .
Before proving that it is subject to a local conservation law, we need to make some remarks. First of all, performing an integration by parts and provided that boundary terms vanish (space-time is unbounded, subject to periodic boundary conditions, . . . ) we have

d D x ω µ (∂ ν ∂ ν ζ, ζ) = -d D x ω µ (∂ ν ζ, ∂ ν ζ) = 0 . (2.12) Since ω µ is skew-symmetric in its two arguments ω µ (∂ ν ζ, ∂ ν ζ) identically vanishes.
We stress that eq. (2.12) holds for any ζ (i.e. off-shell), so the integrands are equal, and hence

ω µ (∂ ν ∂ ν ζ, ζ) = -ω µ (∂ ν ζ, ∂ ν ζ) = 0 .
Secondly, following the same reasoning we have

ω ν (ζ, ∂ ν ∂ µ ζ) = -ω ν (∂ µ ζ, ∂ ν ζ) = -ω ν (∂ ν ζ, ∂ µ ζ) = ω ν (∂ µ ζ, ∂ ν ζ) = 0 .
We now have all the necessary ingredients to prove the local conservation of the stress-energy tensor. Using the relations established above, and eqs. (2.9) and (2.10), it is straightforward to deduce

∂ ν T µν = ∂ ν T νµ = 0 . (2.13)
Hence, the stress-energy tensor is locally conserved (onshell).

Let us now sketch why we have defined the stress-energy tensor as eq. (2.11). We suppose that the action of our theory,

S = d D x L φ i , ∂ µ φ i ,
is invariant under the local transformation

x µ → x µ = x µ + µ , (2.14)
where 1 is constant while may depend on x. The local variation of the field is

δφ i (x) = φ i x -φ i (x) ,
while the global variation is

∆φ i (x) = φ i (x) -φ i (x) .
The local variation only depends on the nature of the field (scalar, vector, tensor, . . . ) while the global variation includes the effect of the transformation. They are related by

∆ = δ -µ ∂ µ + O( ) .
On the other hand, the measure is affected by a Jacobian and reads

d D x = ∂x ∂x d D x = 1 + ∂ µ µ + O 2 d D x .
Then, the variation of the action is

δS = d D x δL + L δd D x = d D x (δL + L ∂ µ µ ) = d D x (∆L + ∂ µ (L µ )) .
The global variation of the Lagrangian density reads

∆L = ∂L ∂φ i ∆φ i + ∂L ∂(∂ µ φ i ) ∆∂ µ φ i = ∂ µ ∂L ∂(∂ µ φ i ) ∆φ i + ∂L ∂(∂ µ φ i ) ∂ µ ∆φ i = ∂ µ ∂L ∂(∂ µ φ i ) ∆φ i .
(using the Euler -Lagrange equation and the fact that the global variation commutes with the space-time derivatives). Hence, the variation of the action

δS = 0 = d D x ∂ µ ∂L ∂(∂ µ φ i ) ∆φ i + L µ ,
defines a conserved Noether current associated to the symmetry of the action under transformation (2.14):

j µ = ∂L ∂(∂ µ φ i ) ∆φ i + L µ = ∂L ∂(∂ µ φ i ) δφ i - ∂L ∂(∂ µ φ i ) ∂ ν φ i -η µν L ν .
Now, if we restrict ourselves to global transformations only ( becomes constant), the Noether current associated with the global translational invariance along is

j µ = η µν L - ∂L ∂(∂ µ φ i ) ∂ ν φ i ν ,
since all the fields (whatever their nature) have no local variation (the Jacobian is the identity). Now, the stressenergy tensor is naturally defined as the collection of the Noether currents, each associated with the global translational invariance along a direction of space-time. It is then defined as

Θ µν = η µν L - ∂L ∂(∂ µ φ i ) ∂ ν φ i ,
and is, by construction, subject to a local conservation law

∂ µ Θ µν = 0 .
Note that this tensor is not necessarily symmetric but it can be symmetrized following the popular Belinfante procedure. See, e.g. [START_REF] Peter | Primordial cosmology, Oxford graduate texts[END_REF].

The stress-energy tensor is expressed in terms of the Hamiltonian density (up to a sign) as

Θ µν = η µν H + ψ i µ ∂ ν φ i -η µν ψ i κ ∂ κ φ i . ω µ (∂ ν ζ, ζ)/2 produces terms of the form ψ i µ ∂ ν φ i -∂ ν φ i ψ i µ /2
, which, after integration by parts give ψ i µ ∂ ν φ i . Then, adding a term to restore its symmetry, the stress-energy tensor is indeed given by eq. (2.11).

To conclude this section, let us now define the charges and prove their conservation. The charges are the conserved quantities associated to the Noether currents. Since the stress-energy tensor is a collection of D Noether currents, it will define D charges as

Q µ = d d x T 0µ .
The stress-energy tensor is locally conserved (eq. (2.13)) and hence the charges are subject to a global conservation law (using Stokes' theorem and again assuming there are no boundary terms):

dQ µ dx 0 = d d x ∂ 0 T 0µ = -d d x ∂ j T jµ = 0 ,
where j runs on 1, d .

Summary

Throughout this section, we have first introduced the notion of symplecticity and highlighted its central rôle in Hamiltonian mechanics. We have then presented the De Donder -Weyl covariant Hamiltonian formulation of field theory and we have shown that it leads to a multisymplectic phase space. We have also introduced (relying on the particular example of the non-linear wave equation) a construction that allows one to obtain a theory, equivalent to the original one, but which does not lead to a degeneracy of the multi-symplectic structure of the phase space. Finally, we have proved the local (on-shell) conservation of multi-symplecticity. We have also constructed the stress-energy tensor and the charges, and we proved their local and global conservation (on-shell), respectively.

To do so, following [START_REF] Bridges | Multi-symplectic structures and wave propagation[END_REF][START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity[END_REF][START_REF]Numerical methods for Hamiltonian PDEs[END_REF] we introduced the general form (2.5) of the De Donder -Weyl -Hamilton equations. Now we claim that, whenever a pde can be written in the form (2.5), and has a non-degenerate multi-symplectic structure, the msilcc method can be applied. We are going to present how to implement it in the next section.

Construction and properties

The msilcc method is a centered box (finite-difference) scheme [START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity[END_REF] except that we do not implement it on the traditional hypercubic lattice. Instead, we use a lattice based on light-cone coordinates, that has the advantage to restore the locality of the method (i.e. there is the same number of unknowns as equations in each cell).

In this section we will present the implementation of the msilcc method as well as a review of some of its interesting properties. We will illustrate it on the example of the λ φ 4 theory mainly in 1 + 1 dimensions.

The lattice: sampling the space-time manifold

Let us consider the new coordinate system

x0 = x 0 √ 2 - 1 √ 2 d j=1 x j , xj = √ 2 x j + δ j d x0 , (3.1)
the inverse of which is

x 0 = 1 √ 2 d µ=0
xµ ,

x j = 1 √ 2 xj -δ j d x0 .
Defining ∂µ = ∂ /∂ xµ , the associated vector basis is

∂0 = ∂ 0 -∂ d √ 2 , ∂j = ∂ 0 + ∂ j √ 2 ,
and its inverse

∂ 0 = ∂0 + ∂d √ 2 , ∂ j = √ 2 ∂j -∂ 0 .
One can remark that (3.1) is not the usual light-cone coordinate system (except in dimension D = 1 + 1). The difference is mainly that the set ∂µ is not an orthogonal basis (while it is for the usual definition).

Let us now sample the space-time manifold M using the lattice

M = n = o + δ ňµ ∂µ ňµ ∈ Z, n ∈ M , (3.2)
where o is arbitrary (chosen such that M respects as much as possible the boundaries of M) and δ is the lattice spacing.

At each point on the lattice (n ∈ M ) we define the elementary cell (the definition can be extended to each point

n C ∈ M + δ C ) cell(n) = n + δ ∂ 0 + σ ∂ ρ √ 2 σ = ±1, ρ ∈ 0, d ,
if all the vertices of the cell belong in M (possibly using periodic boundary conditions). Figure 15 represents how the lattice looks like in low dimensions. On this lattice, the approximation of the derivatives along the original directions of space-time (highlighted in color) will simply be obtained using the midpoint rule. Hence, they only involve two points of M and this makes the method local.

∂ 0 • • • • • D = 0 + 1 • • • • • • • • • ∂ 1 ∂ 0 • • • • D = 1 + 1 • • • • • • • • • • • • • • • ∂ 1 ∂ 2 ∂ 0 • • • • • • D = 2 + 1
Let us now highlight that the direction ∂µ selects in cell(n) one (and only one) square, of width δ, with vertices

n R -| ∂µ = n + δ R -| ∂µ = n n R +| ∂µ = n + δ R +| ∂µ = n R -| ∂µ + δ ∂µ = n + δ ∂ 0 + σ ∂ ρ √ 2 n L -| ∂µ = n + δ L -| ∂µ = n R +| ∂µ -σ δ √ 2 ∂ ρ = n + δ ∂ 0 -σ ∂ ρ √ 2 n L +| ∂µ = n + δ L +| ∂µ = n L -| ∂µ + δ ∂µ = n + δ √ 2 ∂ 0 ,
where the second equation is used to determine σ = ±1 and ρ ∈ 0, d . The center of this square (which is the center of the cell as well),

n C = n + δ C = n + δ ∂ 0 √ 2 , ( 3.3)
is the point where all the approximations are made.

The numerical approximation scheme

Definition

The approximation rules are constructed by a concatenation of the centered box scheme, itself a concatenation of midpoint rules, applied on a square of the elementary cell. Indeed, this scheme is the simplest way to construct a symplectic integrator. Therefore, by combining midpoint rules in the most symmetric possible way, we expect to conserve this property and thus obtain a simple multisymplectic integrator (locally well defined thanks to the lattice).

For ζ or one of its components and for f smooth enough, the approximation rules are given by

f {ζ a } n C ≈ f ζ a n C , ∂µ f {ζ a } n C ≈ ϵ f {ζ a } n C ,
with n C defined in eq. (3.3),

ζ a n C = 1 2D m∈cell(n)
ζ a (m) , and (3.4)

ϵ f {ζ a } n C = 1 δ f 1 2 ζ a n L +| ∂µ + ζ a n R +| ∂µ -f 1 2 ζ a n L -| ∂µ + ζ a n R -| ∂µ .
For the moment we do not know whether these approximation rules respect the rules of differential calculus (we will explore this issue in the following). A priori, the algebraic manipulations done in the continuous formulation will not be equivalent to the ones done on the discrete representation. Hence, the msilcc scheme should be applied only in the light-cone coordinate system (i.e. all the derivatives ∂ µ have to be re-expressed in term of the derivatives ∂µ before applying the scheme). When and only when directly applied on a field, we have

∂ 0 ζ a n C = ∂0 + ∂d √ 2 ζ a n C ≈ [• • •] , ∂ j ζ a n C = 2 ∂j -∂0 -∂d √ 2 ζ a n C ≈ [• • •] .
After some straightforward algebraic manipulations, we obtain

∂ µ ζ a n C ≈ D µ ζ a n C = (3.5) 1 √ 2 δ ζ a n C + δ ∂ µ √ 2 -ζ a n C -δ ∂ µ √ 2 .
Equation (3.5) defines the derivatives of the field along the original directions of space-time as nothing else than the midpoint rule. Nevertheless, remember that this is true only for a linear function of the field, otherwise it is necessary to return to ϵ .

The discrete analogue of the equation of motion (2.5)

in cell(n) at n C is M µ • ζ n + δ ∂ 0 + ∂ µ √ 2 -ζ n + δ ∂ 0 -∂ µ √ 2 = √ 2 δ ∇H 1 2D σ=±1 ρ∈ 0,d ζ n + δ ∂ 0 + σ ∂ ρ √ 2 .
As expected, the approximation of the equation of motion is indeed a concatenation of midpoint rules. Let us illustrate how the procedure works with an example.

Application to the λ φ 4 theory in 0 + 1 dimension

The mechanical problem is described by the two unknowns q and p that only depend on time. We first sample them through M q n = q(t = n δ) ,

p n = p(t = n δ) .
Then, applying the msilcc scheme, we get the discrete version of the equations of motion

p n -p n+1 = δ q n+1 + q n 2 1 + q n+1 + q n 2 2 , q n+1 -q n = δ p n+1 + p n 2 .
In this particular case, one could write down the explicit expressions for (q n+1 , p n+1 ) as functions of (q n , p n ). But these expressions are quite cumbersome and it is not worth presenting them here. This is a symplectic approximation.

The λ φ 4 theory in 1 + 1 dimensions

The lattice is now the same as for the BDdV method (see fig. 5). Again, defining

σ n = 2(n mod 2) -1 ≡ ±1 ,
we sample the fields through M as

ζ a j n = ζ a x = √ 2 δ j + 1 + σ n 4 , t = n δ √ 2 ,
where n ∈ N and j ∈ 0, N . Therefore the discrete version of the equation of motion is given by the set of algebraic equations

ψ 0 j n+1 -ψ 0 j n-1 + σ n ψ 1 j+σn n -σ n ψ 1 j n = - √ 2 δ φ j n-1 + φ j n + φ j+σn n + φ j n+1 4   1 + φ j n-1 + φ j n + φ j+σn n + φ j n+1 4 2   , φ j n+1 -φ j n-1 -σ n γ j+σn n + σ n γ j n = √ 2 δ ψ 0 j n-1 + ψ 0 j n + ψ 0 j+σn n + ψ 0 j n+1 4 , γ j n+1 -γ j n-1 + σ n φ j+σn n -σ n φ j n = - √ 2 δ ψ 1 j n-1 + ψ 1 j n + ψ 1 j+σn n + ψ 1 j n+1 4 , ψ 1 j n+1 -ψ 1 j n-1 -σ n ψ 0 j+σn n + σ n ψ 0 j n = 0 .
These are the equations used to integrate the λ φ 4 theory with the msilcc method in the first section. Again, they are not implicit but too much complicated to write down in an explicit form. Hence they were treated as implicit equations and solved using the Levenberg -Marquardt algorithm for non-linear least squares [START_REF] Nocedal | Numerical Optimization[END_REF].

In this way we generated the data exposed in section 1.

Conservation properties

Leibniz's product rule for quadratic forms

Let us now explore how the approximation rules behave with respect to the rules of differential calculus. We first apply the discrete derivative to a quadratic form. After a straightforward but tedious calculation we obtain

∂µ ζ a ζ b n C ≈ ϵ ζ a ζ b n C = = ζ a ∂µ n C ϵ ζ b n C + ζ b ∂µ n C ϵ ζ a n C ≈ ζ a ∂µ ζ b n C + ζ b ∂µ ζ a n C , (3.6)
where the average value on the square selected by ∂µ is

ζ a ∂µ n C = 1 4 ζ a n L +| ∂µ + ζ a n R +| ∂µ + ζ a n L -| ∂µ + ζ a n R -| ∂µ . ( 3.7)
First of all, eq. (3.6) defines the approximation rule for ζ a ∂µ ζ b such that the Leibniz's product rule for quadratic forms holds (actually the msilcc scheme was designed for that purpose since it is a simple way to construct an approximation that preserves the multi-symplectic structure). As a second remark, the Leibniz's product rule remains valid on the discrete scheme for quadratic forms only. Finally, this is not true for ∂ µ (except in D = 1 + 1 dimensions since ∂µ coincides with by definition). Hence, the necessity to work in the light-cone coordinate system (all the derivative have to be re-expressed in terms of ∂ before making approximations).

To conclude this section, we stress 7 that the same relations hold on the lattice M + δ C .

Preservation of cross derivatives equality

Now, we define the collection z k . Each z k lives on M + δ C and is linear in the field. So the collection z k is 7 When applied on objects of the collection z k (defined in section 3.3.2) the average value on the square selected by ∂µ as well as the full average value are defined on M + δ C at n as

z k ∂µ (n) = 1 4 z k n L +| ∂µ -δ C + z k n R +| ∂µ -δ C + z k n L -| ∂µ -δ C + z k n R -| ∂µ -δ C z k (n) = 1 2D m∈cell(n-δ C ) z k (m) . limited to ζ a , ζ a ∂µ , ϵ ζ a .
The average values were defined in the previous section.

We add here the definition of the derivatives

∂µ f z k (n) ≈ ϵ f z k (n) = 1 δ f 1 2 z k n L +| ∂µ -δ C + z k n R +| ∂µ -δ C -f 1 2 z k n L -| ∂µ -δ C + z k n R -| ∂µ -δ C .
Following the same reasoning as when the field was ζ a defined on M , we have (again, only when directly applied on a field)

∂ µ z k (n) ≈ D µ z k (n) = 1 √ 2 δ z k n + δ ∂ µ √ 2 -z k n -δ ∂ µ √ 2 .
Using these definitions, one can give meaning to the second derivative of the field. We have

D µ D ν ζ a (n) = = 1 √ 2 δ D ν ζ a n + δ ∂ µ √ 2 -D ν ζ a n -δ ∂ µ √ 2 = 1 2 δ 2 ζ a n + δ ∂ µ + ∂ ν √ 2 -ζ a n + δ ∂ µ -∂ ν √ 2 -ζ a n -δ ∂ µ -∂ ν √ 2 + ζ a n -δ ∂ µ + ∂ ν √ 2 = 1 2 δ 2 ζ a n + δ ∂ ν + ∂ µ √ 2 -ζ a n + δ ∂ ν -∂ µ √ 2 -ζ a n -δ ∂ ν -∂ µ √ 2 + ζ a n -δ ∂ ν + ∂ µ √ 2 = 1 √ 2 δ D µ ζ a n + δ ∂ ν √ 2 -D µ ζ a n -δ ∂ ν √ 2 = D ν D µ ζ a (n) ,
proving the identity of the cross derivatives in discrete space-time. Using the relation between D µ and ϵ , we find that the same applies on the light-cone coordinates:

ϵ Ďν ζ a (n) = Ďν ϵ ζ a (n) .

Exact conservation of the multi-symplectic structure

Let us now prove the conservation of the multi-symplectic structure. We first perform the change of coordinates in the left hand side operator of the equation of motion (2.5)

M µ • ∂µ = M µ • ∂ µ ,
and we obtain the set of skew-symmetric matrices in the new coordinate system

M 0 = 1 √ 2 M 0 - 1 √ 2 d j=1 M j , M j = √ 2 M j + δ j d M 0 .
From M µ we define ωµ that actually behaves as a component of a D-vector in

M ωµ = ∂ xµ ∂x ρ ω ρ .
The set { ωµ } defines the multi-symplectic structure in the light-cone coordinate system and is subject to the same conservation law

∂µ ωµ = ∂x ρ ∂ xµ ∂ ρ ∂ xµ ∂x σ ω σ = ∂ ρ ω ρ = 0 .
By taking the exterior derivative of the equation of motion (2.5) we have

M µ ab ∂µ d b = ∂ a ∂ b H(ζ) d b .
Then, the local conservation of multi-symplecticity is, numerically,

∂µ ωµ ≈ ϵ ωµ = = - 1 2 M µ ab ϵ d a ∧ d b ∂µ + d a ∂µ ∧ ϵ d b = M µ ab ϵ d b ∧ d a ∂µ = ∂ a ∂ b H( ζ ) d b ∧ d a ∂µ = 0 ,
since the contraction of the symmetric object ∂ a ∂ b with the skew-symmetry of the wedge product vanishes. So, the multi-symplectic structure is indeed exactly preserved by the msilcc scheme which is hence a multisymplectic integrator.

Conservation of the stress-energy tensor

In this section we investigate the effect of the msilcc scheme on the stress-energy tensor.

Local approximate conservation of the stress-energy tensor

Let us first obtain two preliminary results. On the one hand, we have the commutativity of the oriented average value (eq. (3.7)) with itself

ζ a ∂µ ∂ν (n) = = 1 4 σ=± X={L,R} ζ a ∂µ n -δ C + δ X σ| ∂ν = 1 16 σ=± σ =± X={L,R} X ={L,R} ζ a n -2 δ C + δ X σ| ∂ν + δ X σ | ∂µ = 1 4 σ =± X ={L,R} ζ a ∂ν n -δ C + δ X σ | ∂µ = ζ a ∂ν ∂µ (n) .
On the other hand, we have the commutativity of the oriented average value with the derivative

ϵ ζ a ∂ν (n) = = 1 4 σ=± X={L,R} ϵ ζ a n -δ C + δ X σ| ∂ν = 1 8 δ σ=± σ =± X={L,R} X ={L,R} σ ζ a n -2 δ C + δ X σ| ∂ν + δ X σ | ∂µ = 1 2 δ σ =± X ={L,R} σ ζ a ∂ν n -δ C + δ X σ | ∂µ = ϵ ζ a ∂ν (n) .
Now we consider the non-symmetrized part of the stressenergy tensor

T µν = 1 2 [ω ν (∂ µ ζ, ζ) -η µν ω κ (∂ κ ζ, ζ)] + η µν H(ζ) .
Since T is a tensor (η is a tensor and ω and ∂ are vectors), we have

Ť µν = ∂ xµ ∂x ρ ∂ xν ∂x σ T ρσ = 1 2 ων ∂µ ζ, ζ -ηµν ωκ ∂κ ζ, ζ + ηµν H(ζ) ≈ Ť µν , where ηµν = ∂ ρ xµ ∂ ρ xν .
The numerical version of Ť is defined (using the approximation rules introduced earlier) as Then, using all the preliminary results of this section, the approximation of the local conservation of the stress-energy tensor reads

Ť µν = 1 2 ων ϵ ζ, ζ ∂µ -ηµν ωκ Ďκ ζ, ζ ∂κ + ηµν H( ζ ) .
Ďν Ť µν = 1 2 0 -ϵ ωκ Ďκ ζ, ζ ∂κ + ϵ H( ζ ) = ϵ H( ζ ) - 1 2 ωκ ϵ Ďκ ζ, ζ ∂κ ∂µ + ωκ Ďκ ζ ∂µ , ϵ ζ ∂κ = ϵ H( ζ ) - 1 2 ωκ Ďκ ϵ ζ, ζ ∂µ ∂κ + ωκ Ďκ ζ ∂µ , ϵ ζ ∂κ = ϵ H( ζ ) -dH ζ ϵ ζ = ϵ H I ( ζ ) -dH I ζ ϵ ζ , (3.8)
where

H I = H -H Q ,
is the non-quadratic part of H (H Q is the quadratic part of the Hamiltonian density).

Accordingly, the msilcc scheme exactly preserves the local conservation of the stress-energy tensor for any linear Hamiltonian pde. When applied on a non-linear problem, the msilcc method breaks the conservation of the stress-energy tensor only because the chain rule does not hold on the discrete space-time. Nevertheless, if the sampling is good enough we expect the msilcc integrator not generate large violations of this conservation law.

One can remark that there is no longer any second derivative in eq. Obviously, this operation is strictly forbidden. Nevertheless, eq. (3.9) is a very good estimator of eq. (3.8). This can be understood if we remember that eq. (3.8) mainly evaluates how much the chain rule is violated for nonquadratic functions on the lattice. Therefore, increasing the averaging is not an essential element.

In practice, on the example of the λ φ 4 theory in 1 + 1 dimensions the difference between eq. (3.8) (or explicitly eq. (3.10)) and eq. (3.9) (explicitly eq. (3.11)) is negligible and it is almost impossible to distinguish the two on the numerical results.

The substantial advantage of the estimator (3.9) is that it is simpler to compute, but first and foremost, that it is more local (it involves only the current cell). Thus, the accuracy of the integration can be checked regardless of the neighbouring cells. This ensures a better scalability of the method by reducing the number of communications.

Motivation to use the light-cone coordinates

Solving the numerical equation of motion requires the simultaneous solution of a set of algebraic equations. This is equivalent to finding the root of a vector function. It is generally preferable to see this problem as the minimization of the square norm of a vector function since algorithms for this kind of optimization are more robust and more diversified than the ones for finding roots. Indeed, this can be achieved by using standard optimization methods such as the Levenberg -Marquardt, Powell's Dog Leg, etc. [START_REF] Nocedal | Numerical Optimization[END_REF].

We alluded to this feature earlier, but we now want to stress the importance of the lattice. It has been chosen such that in each cell there is only one point at the latest time. Thus, in each cell we have as many algebraic equations as unknowns. The method is well defined locally. Usually, the centered box scheme is implemented on a hypercubic lattice which is indeed simpler but leads to more unknowns than equations in each cell (except in dimension D = 0 + 1). The method is still globally well defined since each unknown is involved in the equations of the neighbouring cells. However, at each time step, it requires to solve the whole system in one block. Therefore, if we want to dispatch the problem on several process units a huge number of communications are needed (known to be a bottleneck for high performance computations).

The main advantage of the msilcc method, is that it restores the locality of the algorithm while most of the expressions (the equation of motion, the conservation of the stress-energy tensor, . . . ) remain quite simple as we have shown through the example of the λ φ 4 theory.

We finally want to make a remark concerning the initial conditions: the lattice of the msilcc method is such that a cell involves three levels of time. Therefore, at the initial time, in each cell, we have two unknowns for only one equation. The idea to solve this tricky problem is to assume (only at the initial point) that the average in space is equal to the average in time (i.e. the average over all the points of the cell at t = 0 is equal to the average of the two points at t = 0 ± δ C ). In this way we have removed the superfluous unknowns. Nevertheless, it requires that the equation of motion contains a derivative along time of all the fields of the vector state. Hence the necessity to work with a formulation of the problem that will not lead to a degeneracy of the multi-symplectic structure.

Higher dimensions

In the present section we discuss the limits of the msilcc method and possible ways of improvement.

The lattice M , defined in eq. (3.2), is an attempt to generalize to higher dimensions the one introduced in section 1.3.1 for D = 1 + 1 (see fig. 5). However, we have experienced some instabilities of the method in dimension D > 1 + 1. This change in the behavior of our method when the dimension of space-time increases is a bit astonishing. We suspect two reasons for that. When the dimension of space-time becomes higher than 1 + 1:

i. On the one hand, the ensemble of the cells of the lattice is no longer a tessellation of the space-time manifold (i.e. there are points in space-time which are not contained in any cell). ii. On the other hand, the oriented average, • • • ∂µ (see eq. (3.7)), no longer coincides with the full average • • • (eq. (3.4)).

One can imagine another generalization of the lattice introduced in section 1.3.1 which avoids the two problems mentioned above. This is a hypercubic lattice, oriented in such a way that there is only one unknown in each cell (cells are now hypercubes). It consists in starting with another light-cone coordinate system:

∂µ = ∂x ν ∂ xµ ∂ ν = R ν µ ∂ ν ,
where R is a rotation matrix (i.e. R ∈ SO(D)) such that the direction (1, 1, • • • , 1) is mapped to (1, 0, • • • , 0) (in that way, each cell will only contain one unknown). So,

R ν 0 = 1 √ D .
Then, it remains to orthogonalize the remaining rows of R, which can be achieved by defining

R ν µ = 1 µ(µ + 1) ν < µ , R ν µ = - µ µ + 1 ν = µ , R ν µ = 0 ν > µ ,
where µ ∈ 1, d and ν ∈ 0, d .

The space-time is now discretized using an hypercubic lattice rotated by R, and the discretization rules are simple concatenations of midpoint rules as introduced in [START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity[END_REF].

Such a method would have the same property of multisymplecticity and locality while, hopefully, remaining stable in any dimension. But we leave the complete development of it for future work.

Conclusions

The purpose of this paper was to introduce a new numerical method to integrate partial differential equations stemming from the Hamiltonian dynamics of field theories. The method is a centered box scheme, implemented on the light-cone coordinates, in such a way to restore the locality of the algorithm without losing its multi-symplectic properties.

Our method has local conservation properties (and therefore global conservation properties as well) in agreement with what is generally achieved by multi-symplectic integrators. The errors committed do not strongly accumulate, remaining very small over very long periods of time. This is important in applications in which the long-time limit of evolution should be reached with good confidence.

In the process of comparing the performance of our algorithm to other ones in the literature we showed that exact global conservation properties, as the ones imposed in the BDdV technique, do not necessarily guarantee small errors in the local conservation laws.

We highlighted the link between the De Donder -Weyl formalism of field theories and the multi-symplectic structure of phase space, and we treated the latter on a rigorous geometric way. We developed the construction of the stress-energy tensor in the Hamiltonian formalism. We showed that it is exactly conserved in the continuum and we derived the error committed by the algorithm in its discrete implementation. In particular, we showed that it is exactly preserved for a linear equation.

Interestingly, depending on the model that we considered, the multi-symplectic structure was found to be degenerate in spatial dimension larger than zero. We showed how to solve this problem in any dimension using the particular case of the wave equation as an example. The generalization to other field equations should follow similar steps.

Figure 1 :

 1 Figure 1: Space-time plots of the solutions of eq. (1.1) with initial and boundary conditions given by eq. (1.3) and eq. (1.2) respectively. Different panels show data for different values of A, obtained with the msilcc method and L / √ 2 δ = 1024. Lines are iso-levels of the field while color is constant in between.Figure 2 represents a crosssection of these space-time plots for the smallest values of A.

Figure 2 :

 2 Figure 2: Cross-sections of the space-time plots in fig. 1 along the axis x = L /4. Red line is the cross-section of the second eigenmode of the linear theory: A cos 2π t /L .

Figure 4 :

 4 Figure 4: Time evolution of the total energy for A = 10 with the Newton integrator and L /δ = 128.

Figure 5 :

 5 Figure 5: Lattice description of the space-time manifold used in the BDdV method.

Figure 6 :√ 2 δ

 62 Figure 6: Time evolution of the total energy for A = 10 using the BDdV method with L / √ 2 δ = 128. Beyond t /L = 1 the horizontal axis is shown in a different linear scale and data at 255 consecutive instants are skipped between two successive points.

Figure 7 :

 7 Figure 7: Time evolution of the total energy for A = 10 using the msilcc method with L / √ 2 δ = 128. The horizontal axis is the same as in fig. 6. The second part of the graph also shows a histogram of the occurrences of the energy during the integration process (from t /L = 0 to t /L = 100) with 128 bins uniformly distributed on the interval [0.999, 1.001].

Figure 8 :

 8 Figure 8: Error committed by the different methods on the local conservation of the stress-energy tensor as a function of the initial amplitude (discretization is still on 128 points).

Figure 9 :

 9 Figure 9: Error committed by the different methods on the conservation of the charges as a function of the initial amplitude (discretization is still on 128 points).

∆Figure 10 :

 10 Figure10: Long time behavior of the error committed by the different methods on the local conservation of the stress-energy tensor for A = 10 (discretization is still on 128 points). The upper pair of curves are for the BDdV method (triangles), the intermediate ones for the Newton method (diamonds), and the lower ones for the msilcc method (circles). Open and closed symbols show different ways of measuring the error as defined in the text. Beyond t /L = 1 the horizontal axis is shown in a different linear scale and the curves with open symbols are not plotted since they vary too rapidly with respect to this new time-scale (these represent instantaneous errors that in any case are not relevant on this time-scale).

Figure 11 :

 11 Figure 11: Long time behavior of the error committed by the different methods on the conservation of the charges for A = 10 (discretization is still on 128 points). Same symbol convention and time axis as in fig. 10.

Figure 12 :

 12 Figure12:The red line is the exact solution of eq. (1.1) in terms of Jacobi elliptic functions after a time t /L = 1 (a particular solution for different initial condition from the ones used so far). The datapoints represent the field evolved by the different numerical methods (triangles for Newton, diamonds for BDdV and circles for msilcc). On this scale we see no difference between the exact and the numerical solutions obtained with the Newton and msilcc methods.

Figure 14 :

 14 Figure14: Graphic representation of the construction of a non-degenerate theory in dimension D = 1 + 0, = 1 + 1 and = 1 + 2 (the canonical direction is shown with a thick line while the others are not). The fields in the d = 1 case will be introduced in eqs. (2.8a) to (2.8d).

Figure 15 :

 15 Figure 15: A portion of the lattice M in D = 0 + 1, 1 + 1 and 2 + 1 dimensions where we have highlighted the elementary cell.

Now, we use

  the exact conservation of the multisymplectic structure Ďν ων ϵ ζ, ζ ∂µ = 0 = = ων Ďν ϵ ζ, ζ ∂µ ∂ν + ων ϵ ζ ∂ν , Ďν ζ ∂µ , and the dual of the equation of motion (2.5) to prove that ων Ďν ϵ ζ, ζ ∂µ ∂ν = ων Ďν ζ ∂µ , ϵ ζ ∂ν = dH ζ ϵ ζ .

  (3.8). Hence, let us approximate Ďν Ť µν by removing the lowest level average value: Ďν Ť µν ϵ H I (ζ) -dH I ( ζ ) ϵ ζ . (3.9)

3. 4 . 2

 42 The λ φ 4 theory in 1 + 1 dimensionsIn the case of the λ φ 4 theory in 1 + 1 dimensions, and assuming that the extra field γ is free (i.e. used as a One can imagine here to proced in the same way by applying the discretization rules on the stress-energy tensor and then use its conservation as an equation of motion (hence an exact conservation of the stress-energy tensor). However, it would become necessary to evaluate the error committed on the original equation of motion. This would leads to evaluate the quantity:ϵ H I ( ζ ) -dH I ζ ϵ ζ (i.e. eq. (3.8)). The dH I ζ ϵ ζ term is the right hand side of the equation of motion while ϵ H I ( ζ ) arises with the derivatives of the stress-energy tensor. So, whether the discretization is performed on the equation of motion or on the stress-energy tensor, to estimate the quality of the approximation we have to evaluate how much dH I ζ ϵ ζ differs from ϵ H I ( ζ ) in both cases.

  .3) Actually, eqs. (2.1) and (2.2) mean that the Poisson bracket (parametrized by ω ab ) is the inverse of the symplectic 2-form (parametrized by ω ab ). In particular, for a mechanical system and in Darboux coordinates the Poisson bracket reads

Acknowledgements

We thank G. Biroli and A. Ciocchetta for useful discussions. L. F. C. is a member of Institut Universitaire de France.

control parameter), eq. (3.8) explicitly becomes

where ± selects µ = 0 or 1 and

while eq. (3.9) explicitly becomes

(3.11)

Note on the possibility of an exact conservation of the stress-energy tensor

It is actually possible to obtain an exact conservation of the stress-energy tensor. Let us remember the idea of the BDdV method: the discretization rules were applied on the energy instead of the equation of motion, then the constrains on the conservation of the energy were used as an equation of motion. These two procedures are equivalent in the continuum limit, but they are not on the lattice since the rules of differential calculus are no longer fulfilled in the latter setting.