Chrysoula, Constantin Alessandro Charles

Frédéric, Johannes, Kang, Luc, Matthieu, Maxence, Oscar Ruben Pierre

Thibault, Thomas Yifan Sophie

Je Remercie En

Keywords: 4.5 Numerical results

Je remercie en premier lieu ma directrice de thèse, Leticia Cugliandolo, pour ses conseils avisés et la pertinence de ses remarques. Elle a su m'accorder la liberté de réflexion nécessaire à tout travail de recherche et elle m'a laissé l'autonomie propice au développement de mes propres idées, qu'elles fussent bonnes ou moins bonnes. Son dynamisme a été une aide précieuse pour mener ces travaux à leur terme. Ce fut une rencontre enrichissante.

Je remercie mes rapporteurs, Claudio Destri et Malte Henkel, d'avoir accepté cette responsabilité. Je les remercie plus particulièrement pour leur retour constructif qui m'a permis d'améliorer, je l'espère, ce manuscrit. Je remercie également les membres du jury : David Dean, Claude Godrèche et Pascal

Introduction

I Numerical integration of classical conservative field theories

1 Introduction and preliminaries • • • • • • • • • • • • • • • • • • • • • • • • 1.1 Introduction • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 1.2 Geometry preliminaries • • • • • • • • • • • • • • • • • • • • • • • • • •

Introduction

Complex systems gather a broad variety of problems from many different fields (in Physics but also in Biology, in Computer Sciences, in Finance . . . ). They are systems constituted of a large number of interacting degrees of freedom, and therefore, any attempt of an exhaustive description would be out of reach (and irrelevant). However, they can be described in a satisfactory way using a statistical approach. This is precisely the purpose of statistical physics: understand the overall behaviours of a system from its microscopic description, ie how and when collective behaviours emerge.

A statistical system is said at equilibrium when the probability distribution of states (ie the probability that the system being in a given state) coincide with the equilibrium one (eg the Boltzmann distribution). In fact, the concept of equilibrium concerns the statistical properties of the system: it is a characterisation in probability, not of the exact state of the system at a given time. Although this particular situation is quite well understood, in nature, it is the exception rather than the rule. A system which is not in equilibrium is said to be out-of-equilibrium, and, unfortunately, this cannot be tackled with the same framework; such a situation is generally much more difficult to describe. Among all the possibilities for a system to be out-of-equilibrium, there is a particular situation that has especially drawn attention over the last few years, and that is the focus of intense research. The post-quench evolution is one the simplest out-ofequilibrium situation: it consists of preparing the system in a particular state that satisfies a chosen statistics; then, at the initial time, it is suddenly left free to evolve. This procedure can also be seen as abruptly changing a control parameter in a system being at equilibrium. Afterwards, the evolution in time of the statistical observables is obtained by averaging over different copies of the same procedure. The aim of studying the out-of-equilibrium dynamics induced by a quench is to bring out some generic behaviours, and to allow for a better understanding of such systems.

This thesis is devoted to study the post-quench dynamics in chosen circumstances. It is organised in two independent chapters; in the first one, we discuss the numerical tools to deal with conservative field theories in such a (non-perturbative) far-fromequilibrium regime. In the second chapter, we discuss the effects of quenching the control parameter with a finite changing rate in the particular case of the ferromagnetic kinetic Ising models (ie Ising models endowed with microscopic stochastic dynamics for the individual spins).

Chapter I

Numerical integration of classical conservative field theories 1 Introduction and preliminaries

Introduction

While the equilibrium properties of statistical systems are quite well understood (although the computations are not always feasible in practice, especially analytically), there is no such a general framework to deal with their out-of-equilibrium behaviour. Yet, whether they are blocked in a metastable state (for instance, in metallurgy, the state of a metal after a quench hardening process), or they never reach an equilibrium state (because of some external driving phenomenon like in meteorology), inequilibrium systems are rather the exception than the rule.

The out-of-equilibrium dynamics of classical conservative field theories was the initial subject of this thesis. In particular, on the last few years, it has been observed that a field theory, suddenly quenched in a far-from-equilibrium state, may present a remarkable dynamics in some circumstances [START_REF] Aarts | On thermalization in classical scalar field theory[END_REF]: after a fast and complex evolution, the system reach a quasi-equilibrium that does not correspond to the thermal equilibrium state of the problem (observables have converged to an unexpected value). Later, on a much longer time-scale, the system slowly evolves to the correct equilibrium state (it can possibly stay in the pre-thermal state in case the system is an integrable model).

The most interesting feature of this pre-thermalisation dynamics [START_REF] Berges | Nonequilibrium quantum fields and the classical field theory limit[END_REF][START_REF] Berges | Prethermalization[END_REF][START_REF]Isotropization far from equilibrium[END_REF][START_REF] Borsanyi | Thermal features far from equilibrium: Prethermalization[END_REF][START_REF] Boyanovsky | Approach to thermalization in the classical ϕ 4 theory in 1 + 1 dimensions: Energy cascades and universal scaling[END_REF][START_REF] Destri | Ultraviolet cascade in the thermalization of the classical φ 4 theory in 3+ 1 dimensions[END_REF] is that the pre-thermal state can be described by a slightly adjusted version of the equilibrium tools. Afterwards, the evolution from the pre-thermal state to the equilibrium one is sufficiently slow to be considered as infinitely slow; once again, we are able to treat it with the tools of the equilibrium statistical physics [START_REF] Gring | Relaxation dynamics and pre-thermalization in an isolated quantum system[END_REF][START_REF] Kollar | Generalized Gibbs ensemble prediction of prethermalization plateaus and their relation to nonthermal steady states in integrable systems[END_REF]. In such a situation, the part of the dynamics that cannot be described using the equilibrium framework is reduced to only a few times the characteristic time-scale of the problem. See [4, 7-9, 11, 16, 17] for further details.

Understanding in which circumstances such a pre-thermalisation dynamics may happen was the first question of interest. Unfortunately, while trying to address the latter, it appeared that the numerical integration of conservative field theories is not a straightforward task; especially when trying to reach long-times. In particular, the basic approaches are simply not working. This first chapter of the thesis is devoted to the development of a reliable numerical integrator for classical conservative field theories.

Our goal is to solve the equation of motion: a non-linear partial differential equation1 (pde) possessing a particular structure. 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 behaviour 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 magnitude of the error that has been done. Yet, a priori, we cannot predict the accumulation of these errors over many time steps and therefore, 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 observe how reducing the discretisation steps affects the solution; in particular, the idea is to check whether the solution converges to something fix when the steps cancel. This simple approach is an excellent estimator when we try to evaluate how good a numerical method performs on short-or intermediate-time solutions. However, this brute-force approach is inapplicable for long-time solutions.

Another standard option to address the question of the quality of a numerical solver 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, 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 latter, the theory can exhibit some conserved quantities that the numerical solutions should preserve as closely as possible.

Introduction and preliminaries

Symmetries and conserved quantities are very good error estimators. They provide an overall control on the numerical approximations and authorise us to trust (or not) in a numerical solution.

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 (ie arising from De Donder -Weyl -Hamilton equations), that will be defined later. Let us simply state for the moment that they are Lagrangian pdes (ie arising from an Euler -Lagrange equation).

The symmetry group of such a theory generally decomposes as follows: G = T × S × I, where i. T is the space-time symmetry group (eg the Poincaré group). Unfortunately, the discretisation process will generally break this symmetry; however, we shall carefully pay attention to preserve the surviving part: the Lorentz covariance. ii. S is the part associated with the multi-symplectic structure to be explained below.

Its preservation will be in depth explored all through this chapter. iii. I is the internal symmetry group of the theory (eg Z 2 ). This part will be neglected since it does not cause any difficulties in practice. Now, we have to stress that we shall focus on the particular class of finite-difference methods. 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 (eg the space-time manifold). ii. Secondly, a set of discretisation 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 discretisation rules specify how to combine these samples with the aim 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 behaviour 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 [START_REF] Morton | Numerical Solution of Partial Differential Equations: An Introduction[END_REF], and are often adequate. Each method has its own preferred application field. In the kind of problems we shall be interested in, we need to control the very long-time behaviour (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 simpler and faster 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, decreasing the time-step reduces the part of the errors due to the approximation scheme but, in the meantime, it increases the part coming from the accumulation of the truncation errors (inevitably occurring at each step); consequently, the precision cannot be indefinitely improved. This has to be taken into account when trying to reach long-times; hence, the necessity to chose a numerical method designed to behave correctly whatever the number of steps to handle (and which does not require a too small time-step).

Regarding mechanical systems (support is of dimension one, eg just time), there exists a very particular class of finite-difference integrators: the symplectic ones. They are well known (especially by researchers in planetary orbits 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.

Generalisations to field theories (pde) raises further 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 generalised.

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], generalise to pdes the concept of symplectic integrators. Applied to conservative pdes, multi-symplectic integrators exhibit excellent local conservation properties (especially of the stress-energy tensor) and a very stable behaviour 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 multi-symplecticity 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 | High-order 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 [47-49, 54, 56, 71], the non-linear 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]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]Multisymplectic geometry, local conservation laws and Fourier pseudospectral discretization for the[END_REF], as well as the Zakharov -Kuznetsov (ZK) equation [START_REF] Chen | Multisymplectic geometry, local conservation laws and a multisymplectic integrator for the Zakharov-Kuznetsov equation[END_REF].

However, none of these methods fulfil all of our requirements; they all only consider the S part of the symmetry group, completely forgetting about the space-time part of it (T ). The only one which fortuitously does not break the covariance of the theory (the centred box scheme [START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian {PDEs} that conserve symplecticity[END_REF]) suffers from a severe scalability issue.

The aim of this work is to introduce a new finite-difference multi-symplectic method, based on the centred 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]; including compliance with the Lorentz covariance. 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, causing a severe scalability issue. The idea we introduce in this chapter, inspired by [START_REF] Boyanovsky | Approach to thermalization in the classical ϕ 4 theory in 1 + 1 dimensions: Energy cascades and universal scaling[END_REF], is to use a rotated lattice in the light-cone coordinates; this restores the locality of the algorithm, drastically enhancing its efficiency without affecting its high precision.

We organise the presentation in a pedagogical way; first discussing the theoretical 1 Introduction and preliminaries justification of the method, and next showing how it performs compared to other ones in the market. We also discuss, without any assumption on the dimension of spacetime, the problem of the degeneracy of the multi-symplectic structure and we show how to solve it in the particular case of the non-linear wave equation.

The outline of the chapter is the following.

The remaining of this preliminary part, sections 1.2 to 1.4, will be devoted to some reminders of differential geometry and Hamiltonian mechanics (both in its theoretical and numerical aspects).

The second part of this chapter, 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 our multi-symplectic 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 4, is a particular case).

Next, in section 3, the method we have developed: the multi-symplectic integrator in the light-cone coordinates (msilcc), will be introduced in detail. We shall also demonstrate its conservation properties. Again, the non-linear wave equation will be our working example.

Finally, in section 4, we shall compare the msilcc integrator to two standard methods:

i. On the one hand, a very basic scheme based on the Euler approximation of derivatives (this method is widely used by a broad community and proves to be preserving the multi-symplectic structure too). 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 (non-local conservation); we will show that this is unfortunately not a guarantee of quality.

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

A short conclusions section will close the first chapter of this thesis.

Geometry preliminaries

The constructions presented throughout this chapter will require some elements of symplectic geometry. This section introduces the necessary concepts of differential geometry, without being intended to be an exhaustive presentation; it is inspired from [START_REF] Nakahara | Geometry, topology, and physics[END_REF] to which we refer for a complete review.

First, the concept of tangent space to a differentiable manifold is introduced. Next, the notion of dual space, of cotangent space and objects called co-vectors are defined. Secondly, the tensor product and, in particular, a special kind of tensors called forms will be defined. Later, another noteworthy tensor: the metric, that connects the objects of a space to the objects of its dual will be introduced. Afterward, the exterior derivative, the gradient, the interior product as well as the Lie derivative will be defined.

Note that, unless explicitly stated, Einstein's summation convention will be used.

Numerical integration of classical conservative field theories

Ω Z T z Ω z ∂ 2 ∂ 1 Figure I.
1 -Illustration of the tangent space to a 2-dimensional manifold. T z Ω is the tangent space to Ω at z. The basis vectors, ∂ 1 and ∂ 2 , can be seen as velocity vectors for a point moving respectively along the curvilinear coordinate z 1 or z 2 . z 1 , z 2 is a local coordinate system that parametrises Z, an open neighborhood of z on Ω. In this example, the coordinate system z 1 , z 2 is in fact globally defined on Ω.

Tangent space: vectors

Let us consider a differentiable manifold 1 , Ω, of finite dimension n.

The tangent space to Ω, at any point z ∈ Ω, denoted 2 T z Ω, is the set of all the vectors tangent to the manifold at z. T Ω is isomorphic to R n ie is a n-dimensional vector space.

Since Ω is a manifold, it can be described, at least locally, by a coordinate system, {z a }, where a runs from 1 to n. Hence, the set

∂ a := ∂ ∂z a a∈ 1,n (1.2.1)
is a vector basis for T Ω. This basis can be assumed as orthonormal provided that the local coordinate system {z a } is correctly chosen on the manifold.

To see this, let us consider Z, an open neighborhood of z on Ω. Since this open set can be as small as necessary, and since Ω is a manifold, Z is isomorphic to R n ; Z is therefore covered by a well-defined coordinate system {z a }. On the other hand, the derivative with respect to the curvilinear coordinate z a results in a velocity vector which is tangent to the trajectory swept by z a . Hence, ∂ a at z is a tangent vector to Z at z. By repeating this construction for each coordinate, we construct a set of n tangent vectors. This set is linearly independent provided that the coordinate system {z a } is well-defined at z. Hence, these n tangent vectors, the set (1.2.1), form a basis for T z Ω. Having defined a basis for the tangent space, any vector v ∈ T Ω can be decomposed in coordinates as v = v a ∂ a .

(1.2.2)

1 also called a smooth manifold, it is a topological manifold coming with a notion of differentiability [START_REF] Milnor | Topology from the Differentiable Viewpoint[END_REF]. 2 if the point where the tangent space is taken is not a relevant information in the context, it will be omitted, and the tangent space to Ω, at any point, will be denoted T Ω. This is an abuse of notation, and T Ω should not be confused with the tangent bundle, commonly denoted by mathematicians in this way (the tangent bundle being the disjoint union of the tangent spaces ie z∈Ω TzΩ).

Let ∂ a := ∂ /∂z a and ∂a := ∂ /∂ za be two basis of the tangent space, arising from two curvilinear coordinate system on the manifold: {z a } and {z a }. One can pass from one to the other using the chain rule:

∂ a = ∂ ∂z a = ∂ zb ∂z a ∂ ∂ zb = ∂ a zb ∂b , and (1.2.3a) ∂a = ∂a z b ∂ b . (1.2.3b)
Therefore, the coordinates of the vector v = v a ∂ a = ṽa ∂a ∈ T Ω transform as v a = ∂b z a ṽb , and (1.2.4a)

ṽa = ∂ b za v b , (1.2.4b)
which is known as how vectors transform.

Dual space

Let us now consider a vector space, V , of finite dimension n, and let K be the field of the scalars. V always admits a dual space, denoted V * , defined as the vector space of all the linear maps from V to K.

Among other properties, as long as V is finite-dimensional, V * has the same dimension. The dual of the dual is isomorphic to the vector space itself, again, as long as it is finite-dimensional. This isomorphism is in practice strongly restrictive: the possible bijections between (V * ) * and V are almost reduced to the identity.

Let B := {e a } be a basis in V , where a runs from 1 to n, and let us denote B * := {e a } a basis in V * . One can show that, for any basis B in V , there exists a particular basis B * in V * , namely the dual basis, such that

e b (e a ) := δ b a , (1.2.5a)
where δ is the Kronecker symbol. Since the dual of the dual is almost the vector space itself, the reciprocal holds: for any basis B * in V * , there exists a particular basis B in V such that e b (e a ) := δ a b .

(1.2.5b)

From now, let B * := {e a } be the dual basis to B := {e a }. Thus, for any pair v = v a e a ∈ V and w = w a e a ∈ V * , the action of the one on the other expands in coordinates as

w(v) = v(w) = v a w b δ b a = v a w a . (1.2.6)
However, at this stage, there is no relation between v and w; they live in different spaces. The connection between the elements of V and the elements of its dual relies on the existence of a metric tensor (which will be defined later).

Cotangent space: co-vectors

Let us return to the notion of tangent space.

Since T Ω is a finite-dimensional vector space, it possesses a dual, called the cotangent space, which has thus the same dimension. This cotangent space is denoted (T Ω) * =: T * Ω and its elements are called co-vectors.

Let now B := {∂ a } be the basis in T Ω arising from the curvilinear coordinate system {z a } (on Ω, or, at least, on an open neighborhood of z). From now, B * , the dual basis in T * Ω (ie the basis for co-vectors), will be denoted

{d a := dz a } a∈ 1,n . (1.2.7)
The action of a co-vector on a vector (or vice versa) is respectively obtained from eqs. (1.2.5a) and (1.

∂b (d a ) = T a c ∂b ( dc ) = T a c δ c b = T a b = ∂b z c ∂ c (d a ) = ∂b z c δ a c = ∂b z a , and ∂ b ( da ) = T a c ∂ b (d c ) = T a b = ∂ b zc ∂c ( da ) = ∂ b za .
Hence, eqs. (1.2.9a) and (1.2.9b) read d a = ∂b z a db , and (1.2.10a)

da = ∂ b za d b , (1.2.10b)
meaning that co-vectors transform in an inverse way, compared to vectors, under a reparameterisation of the manifold. Therefore, a co-vector w = w a d a = wa da ∈ T * Ω has its coordinates that transform as w a = ∂ a zb wb , and (1.2.11a)

wa = ∂a z b w b . (1.2.11b)
Tensor product: tensors and p-forms

The tensor product 1 of vector spaces is a commutative bilinear product that allows one to obtain a bigger vector space in such a way that the dimension of the product is the product of the dimensions.

Let us consider N vector spaces, {V i }, where i runs from 1 to N , and M co-vector spaces, {W j }, where j runs from 1 to M , and let

e (i)
a i a i ∈ 1,dim V i and e b j (j) b j ∈ 1,dim W j (1.2.12)
1 or direct product.

1 Introduction and preliminaries be, respectively, any basis in V i and W j .

The tensor product of these spaces, denoted

U M N := N i=1 V i ⊗ M j=1 W j = V 1 ⊗ V 2 ⊗ • • • ⊗ V N ⊗ W 1 ⊗ W 2 ⊗ • • • ⊗ W M , (1.2.13a) is a vector space of dimension dim U M N = N i=1 dim V i M j=1 dim W j . (1.2.13b)
The objects living in U M N are called tensors of order (N, M ). A tensor of order (1, 0) is simply a vector; a tensor of order (0, 1) is a co-vector; while a tensor of order (1, 1) is a matrix. In fact, in the coordinates of the tensor, N is the number of upper indices while M is the number of lower indices.

Since U M

N is a finite-dimensional vector space, it possesses a basis:

   N i=1 e (i) a i ⊗ M j=1 e b j (j) = e (1) a 1 ⊗ e (2) a 2 ⊗ • • • ⊗ e (N ) a N ⊗ e b 1 (1) ⊗ e b 2 (2) ⊗ • • • ⊗ e b M (M )    (1.2.14)
where a i ∈ 1, dim V i for each i ∈ 1, N , and where b j ∈ 1, dim W j for each j ∈ 1, M . Hence, any tensor t ∈ U M N can be decomposed in coordinates as

t = t a 1 a 2 ••• a N b 1 b 2 ••• b M N i=1 e (i) a i ⊗ M j=1 e b j (j) . (1.2.15)
In the following, we shall only consider tensor products of the tangent and the cotangent spaces

T M N Ω := (T Ω) ⊗N ⊗ (T * Ω) ⊗M = N i=1 T Ω ⊗ M j=1 T * Ω , (1.2.16)
and a tensor t ∈ T M N Ω will be decomposed in coordinates as

t = t a 1 a 2 ••• a N b 1 b 2 ••• b M N i=1 ∂ a i ⊗ M j=1 d b j . ( 1
∂ a i ⊗ M j=1 d b j = N i=1 ∂ a i zã i M j=1 ∂b j z b j N i=1 ∂ã i ⊗ M j=1
db j , and (1.2.18a)

N i=1 ∂ã i ⊗ M j=1 db j = N i=1 ∂ã i z a i M j=1 ∂ b j zb j N i=1 ∂ a i ⊗ M j=1 d b j . (1.2.18b)
Then, the transformation rules for the coordinates of t are straightforward: in

t a 1 a 2 ••• a N b 1 b 2 ••• b M
, the upper indices are contravariant (ie they transform as the coordinates of a vector), while the lower indices are covariant (ie they transform as the coordinates of a covector).

Let us now define the symmetric product: .2.19b) as well as the anti-symmetric product, or wedge product:

N i=1 ∂ (i) a i = ∂ (1) a 1 ∨ ∂ (2) a 2 ∨ • • • ∨ ∂ (N ) a N := 1 N ! p∈P N (a 1 ••• a N ) N i=1 ∂ (i) p i , and (1.2.19a) M j=1 d b j (j) := 1 M ! p∈P M (b 1 ••• b M ) M j=1 d p j (j) , ( 1
N i=1 ∂ (i) a i = ∂ (1) a 1 ∧ ∂ (2) a 2 ∧ • • • ∧ ∂ (N ) a N := 1 N ! p∈P N (a 1 ••• a N ) σ(p) N i=1 ∂ (i) p i , and (1.2.20a) M j=1 d b j (j) := 1 M ! p∈P M (b 1 ••• b M ) σ(p) M j=1 d p j (j) , (1.2.20b)
where P M is the symmetric group of M symbols, and σ(p) is the signature of the permutation p. Note that ∨ and ∧ are not equivalent to ⊗ since the product on the right hand side is not a commutation of the product on the left hand side; the permutation apply on the indices, not on the vectors (see the indices in gray parenthesis1 ).

Since the permutation apply on indices, symmetrisation and anti-symmetrisation can be equivalently2 defined through the coordinates of the tensor. Symmetrisation of contravariant indices, and respectively, of covariant indices, is defined through the symmetric product as

t (a 1 a 2 ••• a N ) N i=1 ∂ a i := t (a 1 a 2 ••• a N ) N i=1 ∂ a i = t a 1 a 2 ••• a N N i=1
∂ a i , and (1.2.21a)

t (b 1 b 2 ••• b M ) M j=1 d b j := t (b 1 b 2 ••• b M ) M j=1 d b j = t b 1 b 2 ••• b M M j=1 d b j , (1.2.21b)
while anti-symmetrisations are defined from

t [a 1 a 2 ••• a N ] N i=1 ∂ a i := t [a 1 a 2 ••• a N ] N i=1 ∂ a i = t a 1 a 2 ••• a N N i=1 ∂ a i , and (1.2.22a) t [b 1 b 2 ••• b M ] M j=1 d b j := t [b 1 b 2 ••• b M ] M j=1 d b j = t b 1 b 2 ••• b M M j=1 d b j . (1.2.22b)
Therefore, the symmetrisation brackets,

t (a 1 a 2 ••• a N ) := 1 N ! p∈P N (a 1 ••• a N )
t p , and (1.2.23a)

t (b 1 b 2 ••• b M ) := 1 M ! p∈P M (b 1 ••• b M ) t p , (1.2.23b)
while the anti-symmetrisation brackets,

t [a 1 a 2 ••• a N ] := 1 N ! p∈P N (a 1 ••• a N )
σ(p) t p , and (1.2.24a)

t [b 1 b 2 ••• b M ] := 1 M ! p∈P M (b 1 ••• b M ) σ(p) t p . (1.2.24b)
These brackets symmetrise/anti-symmetrise all enclosed indices. We also define the protection brackets, Among all possible tensors, there is a noteworthy type, called p-forms, that are antisymmetric tensors of order (0, p). The formalism later presented mainly relies on 2-forms:

ω = ω [ab] d a ⊗ d b = ω [ab] d a ∧ d b = ω ab d a ∧ d b , (1.2.25)
where ω ab = ω [ab] is anti-symmetric in its two indices.

Metric and dual vectors

Let us consider now another particular kind of tensors. A metric tensor1 is a symmetric non-degenerated tensor of order (0, 2).

Let η = η (ab) d a ⊗ d b = η (ab) d a ∨ d b = η ab d a ∨ d b be a metric on 2 Ω.
First, since the metric is a symmetric non-degenerated application that maps two vector to a scalar, it defines a scalar product. Let u = u a ∂ a and v = v a ∂ a be two vectors of T Ω. Their scalar product is defined as

u, v := η(u, v) = η ab u a v b . (1.2.26)
Since the metric is never degenerated, it can be inverted; its inverse,

η -1 = η ab ∂ a ∨ ∂ b (such that η ab η bc = δ c a )
, defines a scalar product for co-vectors. On the other hand, while the metric acts on a vector to produce a co-vector, its inverse acts on a co-vector to produce a vector. They connect the objects of the tangent space with the objects of its dual in a bijective way defining the notion of dual for vectors/co-vectors. They act on the basis vectors/co-vectors as

η(∂ a ) = η bc d b ∨ d c (∂ a ) = η bc δ c a d b = η ab d b =: (∂ a ) * , and (1.2.27a) η -1 (d a ) = η bc ∂ b ∨ ∂ c (d a ) = η bc δ a c ∂ b = η ab ∂ b =: (d a ) * . (1.2.27b) Hence, v = v a ∂ a ∈ T Ω and w = w a d a ∈ T * Ω are dual provided that v * = (v a e a ) * = v a (e a ) * = v a η ab e b = w b e b = w , ie (1.2.28a) v b = η ab v a = w b , (1.2.28b)
which means that the dual of a vector, is the vector of the dual space having the same coordinates in the dual basis. In coordinates, we observe that the metric, η ab , is used to lower indices while its inverse, η ab , raise indices. In particular, η ab η bc = η c a = δ c a .

Numerical integration of classical conservative field theories

Gradient and exterior derivative

The exterior derivative, d, generalises to forms the notion of differential. It maps a p-form,

f = f a 1 a 2 ••• ap p j=1 d a j , (1.2.29)
to a p + 1-form:

df := df a 1 a 2 ••• ap ∧ p j=1 d a j , (1.2.30a)
where df a 1 a 2 ••• ap is the differential: applied on a function f : Ω → K (a scalar field on Ω), it reads

df := ∂f ∂z a d a = ∂ a f d a . (1.2.30b)
Using its definition and Schwarz's theorem1 , one can show that the exterior derivative of the exterior derivative vanishes: dd := 0. A p-form that vanishes under exterior derivative (df = 0) is said to be closed.

We shall also use the gradient of the scalar fields on Ω; it is the dual of the differential

∇f := (df ) * = ∂ a f (d a ) * = ∂ a f η ab ∂ b . (1.2.31)

Interior product and Lie derivative

Let us finally introduce the notion of interior product and Lie derivative.

Let f be a p-form on Ω while v is a vector. The interior product of v and f is the

p -1-from (ι v f )• := f (v, •) , (1.2.32)
that is the contraction of the first index in the form with the index of the vector.

Let now v = {v(z) ∈ T z Ω|∀z ∈ Ω} be a vector field 2 on Ω. The Lie derivative with respect to the vector field v is defined following the Cartan's identity as the anti-commutator of the interior product and the exterior derivative

L v • := (ι v , d)• = ι v (d•) + d(ι v •) . (1.2.33)
The Lie derivative measures the variations along the integral curves of the vector field.

This concludes the definition of the necessary geometric tools. Let us now recall the essential steps of the construction of the Hamiltonian formulation of mechanics.

1 Introduction and preliminaries

Reminders of Hamiltonian mechanics

This section will be devoted to remind some essential aspects of the Hamiltonian formulation of mechanics. This presentation is inspired from [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF][START_REF] Goldstein | Classical Mechanics, Addison-Wesley series in physics[END_REF][START_REF] José | Classical Dynamics: A Contemporary Approach[END_REF][START_REF] Landau | Mechanics[END_REF] to which we refer for further details. Generalisations to field theories will be presented later.

We first recall the Hamilton's equations of motion and how they emerge from the least action principle. Next, we define the phase space and recall the proof of the conservation of its symplectic structure. Secondly, we recall the Poisson formulation of the Hamilton's equations. Then, the notion of symplectic invariant is introduced via the canonical transformations of the phase space. Later, we shall recall the proof of Liouville's theorem for statistical mechanics of Hamiltonian systems. Afterwards, we shall briefly comment on how to generalise this construction in case of a timedependent Hamiltonian.

Lagrangian, action and Euler -Lagrange equation of motion

Consider a conservative mechanical problem parametrised by a time t and that consists of N degrees of freedom (eg a d-dimensional problem with N/d particles, . . . ). Let q i (t) and qi := dq i /dt respectively be the generalised position and velocity of the i th degree of freedom (i running from 1 to N ).

The dynamics are characterised by a Lagrangian, L q i , qi , that obviously depends on the generalised positions and velocities but, as an assumption, not explicitly on time (the time dependence is only thought q i and qi ):

∂L ∂t = 0 . (1.3.1)
The action reads

S q i := L q i , qi dt . (1.3.2)
Given the dependencies of the Lagrangian, the variation of the action is

δS = ∂L ∂q i - d dt ∂L ∂ qi δq i , (1.3.3)
where δS is the variation of the action arising from a variation in the trajectory, δq i (t) (see [START_REF] Gelfand | Calculus of Variations[END_REF] for a rigorous definition and a complete presentation of the calculus of variations). The least action principle, that is mainly δS = 0, leads to the following Euler -Lagrange equations1 of motion:

∂L ∂q i - d dt ∂L ∂ qi = 0 , (1.3.4)
where i ∈ 1, N .

Hamiltonian formulation

For each generalised position one may define a conjugate momentum

p i := ∂L ∂ qi . (1.3.5)
The Hamiltonian is defined by the following Legendre transformation that substitute the conjugate momentum p i to the generalised velocity qi provided that it is not a singular transformation:

H q i , p i := p i qi -L q i , qi . ( 1
Ḣ = 0 . (1.3.9c)
The Legendre transformation (1.3.6) substitutes the N 2 nd order differential Euler -Lagrange equations of motion by a set of 2N differential equations of first order (the third Hamilton's equation, the one that concerns H, is in facts trivial).

Phase space and symplectic structure

The phase space, denoted Ω, is a differentiable manifold of dim Ω = n = 2N that is parametrised by the Darboux coordinates: {z a } :≡ q i ∪ {p i } (the union of the positions and momenta).

Since Ω is a differentiable manifold we can define on it the 2-from

ω = ω ab d a ∧ d b := dq i ∧ dp i -dp i ∧ dq i = 2 dq i ∧ dp i . (1.3.10) It is a closed (dω = 0) non-degenerated (det ω = 0) 2-from. Hence, (Ω, ω) is what is called a symplectic manifold
where ω is the symplectic form on it. This symplectic structure is an essential part of the Hamiltonian formulation of mechanics;

explicit physical examples will be given below, page 25.

Let us now define

z ∈ T Ω z = z a ∂ a := q 1 • • • q N p 1 • • • p N t . (1.3.11)
As long as {z a } is a global coordinate system on Ω, assumption which we shall admit from now on, for each z ∈ T Ω there is an unique associated point z ∈ Ω. z is called the state vector since it selects a unique state of the system. In fact, with this particular vector field we provide some kind of vector space structure to Ω.

Introduction and preliminaries

Using the symplectic form (1.3.10) and the state vector (1.3.11), the Hamilton's equations of motion are equivalent to

ι ż ω = ω( ż, •) = dH . (1.3.12)
Indeed, since For any smooth function f on Ω there exists a vector field, χ f , defined as

ι ż ω = ω( ż, •) = qi dp i -ṗi dq i ,
ω(χ f , •) := df . (1.3.13)
Hence, the Hamiltonian admits such a vector field, the Hamiltonian vector field, χ H , and from eq. (1.3.13) it follows that the equation of motion also reads

ż = χ H . (1.3.14)
We also define the Hamiltonian flow, Φ, as the integral curves of the Hamiltonian vector field

∇Φ := (dΦ) * = χ H . (1.3.15)
A fundamental property of the symplectic form is that it is conserved under the Hamiltonian flow: using the closeness of the symplectic form as well as eqs. (1.3.13) and (1.3.14), one successively obtains

L χ H ω = ι χ H dω + d(ι χ H ω) = d(ω(χ H , •)) = d(dH) = 0 . (1.3.16)

Poisson bracket

Since ω is non-degenerated it can be inverted and its "inverse" reads

ω := η -1 ⊗ η -1 (ω) = ω ab ∂ a ∧ ∂ b , (1.3.17)
where ω ab ω bc = -δ c a . Here we have used the symbol instead of the usual -1 since ω (ω) = ω ω = -n (due to the sign, they are not truly inverses of each other). Let f q i , p i , g q i , p i and h q i , p i be three smooth functions on Ω.

The Poisson bracket of f and g,

{f, g} = -{g, f } := ω(χ f , χ g ) = L χ f g = ω (df, dg) = ω ab ∂ a ∧ ∂ b (df, dg) = ω [ab] ∂ a f ∂ b g = ∂f ∂q i ∂g ∂p i - ∂f ∂p i ∂g ∂q i , (1.3.18)
is related to (it is sometimes referred to as being) the inverse of the symplectic form.

It is an anti-symmetric bilinear application that respects the product rule {f g, h} = f {g, h} + {f, h} g , (1.3.19) and that is subject to a Jacobi identity:

{{f, g}, h} + {{g, h}, f } + {{h, f }, g} = 0 . (1.3.20)
The Poisson bracket notoriously plays a fundamental role in Hamiltonian mechanics. Most of the properties that give the Poisson bracket such a fame are in fact properties of the symplectic form, however these two objects are so closely related that they merge into the same concept. Accordingly, in the following, these two objects will be interchangeably used depending on which one is the more suitable to highlight a given property.

Earlier we have introduced the action of the Hamiltonian flow in terms of symplectic form, let us now express it in terms of the Poisson bracket. We consider a smooth function f := f q i , p

Introduction and preliminaries

As a final remark, we stress that the Poisson brackets of the Darboux coordinates have a noteworthy form q i , q j = 0 , (1.3.24a) p i , p j = 0 , and (1.3.24b)

q i , p j = δ i j , (1.3.24c)
usually called canonical Poisson brackets.

Canonical transformations: symplectic geometry

Canonical transformations are special deformations of the phase space that preserve the structure of the Hamilton's equations. Let us consider the generic transformation

q i → Q i q j , p j , t , (1.3.25a)
p i → P i q j , p j , t , and (1.3.25b)

H q i , p i , t → K Q i , P i , t . (1.3.25c)
Recalling that the Hamilton's equations in the Darboux coordinates read

qi = ∂H ∂p i , (1.3.26a) ṗi = - ∂H ∂q i , and (1.3.26b)
Ḣ = ∂H ∂t , (1.3.26c)
δS = 0 = δ p i qi -H dt = δ P i Qi -K dt . (1.3.28)
Here we have written the least action principle in Hamiltonian form1 for both the original and the transformed coordinates (and Hamiltonian). This constraint holds provided that the integrands only differ by a total derivative or a scaling factor:

λ p i qi -H = P i Qi -K + dF dt . (1.3.29)
The scaling transformation, represented by λ, merely stands for a change of units and it can be assumed that λ = 1 without loss of generality. Now, observing eq. (1.3.29) we remark that F := F 1 q i , Q i , t is necessarily a function of q i , Q i and t. Its differential reads dF = p i dq i -

P i dQ i + (K -H) dt (1.3.30a) = ∂F 1 ∂q i dq i + ∂F 1 ∂Q i dQ i + ∂F 1 ∂t dt . (1.3.30b)
Thus, we are able to generate a canonical transformation provided that there exists a function F 1 such that

p i = ∂F 1 ∂q i , (1.3.31a
) There exists four kinds of generating functions: here, F 1 is a function of q i , Q i and t. However, by Legendre transformation, we can substitute to F 1 a function that depends on t and on any couple selected in q i , Q i , p i , P i and use it to generate a canonical transformation. The four kinds of generating functions are listed in table I.1.

P i = - ∂F 1 ∂Q i ,

Generating function Transformation equations

F := F 1 q i , Q i , t p i = ∂F 1 ∂q i and P i = - ∂F 1 ∂Q i F := F 2 q i , P i , t -Q i P i p i = ∂F 2 ∂q i and Q i = ∂F 2 ∂P i F := F 3 p i , Q i , t + q i p i q i = - ∂F 3 ∂p i and P i = - ∂F 3 ∂Q i F := F 4 ({p i , P i }, t) + q i p i -Q i P i q i = - ∂F 4 ∂p i and Q i = ∂F 4 ∂P i
Table I The eq. (1.3.31c) has a particular status and is called the Hamilton -Jacobi equation

K = λ H + ∂F ∂t , (1.3.32)
where we have restored the λ parameter (previously set to 1).

Let us now consider a time-independent deformation of the phase space, ie the generating function does not depend explicitly on time:

∂F ∂t = 0 . (1.3.33)
1 Introduction and preliminaries Hence, the transformation (1.3.25a) to (1.3.25c) becomes q i → Q i q j , p j , (1.3.34a) p i → P i q j , p j , and (1.3.34b)

H q i , p i → K Q i , P i . (1.3.34c)
Using the Hamilton -Jacobi eq. (1.3.32), and setting again λ = 1, one obtains K = H. Meaning that the transformed Hamiltonian is the original one expressed in the new coordinates and hence, that K and H are functions of q i , p i as well as of Q i , P i . Using this, the equations of motion for the new coordinates successively read

Qi = ∂Q i ∂q j qj + ∂Q i ∂p j ṗj = ∂Q i ∂q j ∂H ∂p j - ∂Q i ∂p j ∂H ∂q j = ∂Q i ∂q j ∂H ∂Q k ∂Q k ∂p j + ∂H ∂P k ∂P k ∂p j - ∂Q i ∂p j ∂H ∂Q k ∂Q k ∂q j + ∂H ∂P k ∂P k ∂q j = ∂H ∂Q k ∂Q i ∂q j ∂Q k ∂p j - ∂Q i ∂p j ∂Q k ∂q j + ∂H ∂P k ∂Q i ∂q j ∂P k ∂p j - ∂Q i ∂p j ∂P k ∂q j = ∂H ∂Q j Q i , Q j qp + ∂H ∂P j Q i , P j qp ,
Q i , Q j qp = 0 , (1.3.37a)
{P i , P j } qp = 0 , and (1.3.37b)

Q i , P j qp = -P j , Q i qp = δ i j . (1.3.37c)
This is a necessary and sufficient condition.

Let us now consider two smooth functions, f := f q i , p i = f Q i , P i and g :=

g q i , p i = g Q i , P i . One successively obtains {f, g} qp = ∂f ∂q i ∂g ∂p i - ∂f ∂p i ∂g ∂q i = ∂f ∂Q j ∂Q j ∂q i + ∂f ∂P j ∂P j ∂q i ∂g ∂Q k ∂Q k ∂p i + ∂g ∂P k ∂P k ∂p i - ∂f ∂Q j ∂Q j ∂p i + ∂f ∂P j ∂P j ∂p i ∂g ∂Q k ∂Q k ∂q i + ∂g ∂P k ∂P k ∂q i = ∂f ∂Q j ∂g ∂Q k ∂Q j ∂q i ∂Q k ∂p i - ∂Q j ∂p i ∂Q k ∂q i + ∂f ∂P j ∂g ∂P k ∂P j ∂q i ∂P k ∂p i - ∂P j ∂p i ∂P k ∂q i + ∂f ∂Q j ∂g ∂P k ∂Q j ∂q i ∂P k ∂p i - ∂Q j ∂p i ∂P k ∂q i + ∂f ∂P j ∂g ∂Q k ∂P j ∂q i ∂Q k ∂p i - ∂P j ∂p i ∂Q k ∂q i = ∂f ∂Q i ∂g ∂Q j Q i , Q j qp + ∂f ∂P i ∂g ∂P j {P i , P j } qp + ∂f ∂Q i ∂g ∂P j Q i , P j qp + ∂f ∂P i ∂g ∂Q j P i , Q j qp (1.3.38) = ∂f ∂Q i ∂g ∂P i - ∂f ∂P i ∂g ∂Q i = {f, g} QP , (1.3.39)
where the last line assumes the canonicity of the transformation that maps q i , p i to Q i , P i . Hence, the Poisson bracket is invariant under a time-independent canonical transformation. This is also true for the symplectic structure since the Poisson bracket and the symplectic form are two equivalent objects. A change of coordinates that preserves the symplectic form is called a symplectic transformation. The set of all symplectic transformations forms the symplectic group.

However, we must stress an important point on the terminology: canonical and symplectic transformations are not the same thing. A canonical transformation is a deformation of phase space that preserves the structure of Hamilton's equations. On the other hand, a symplectic transformation preserves the symplectic form (ie the Poisson bracket). The two concepts collapse only for the particular case of a timeindependent transformation. For a change of coordinates that depends explicitly on time, they are not equivalent anymore, at least not in the framework introduced until now.

Canonical and symplectic transformations are conciliated by the time-dependent Hamiltonian formalism that modifies the definition of the phase space as well as the definition of the symplectic form to handle the explicit dependence on time; this will be briefly discussed at the end of the section. In this framework, canonical and symplectic transformations stand for the same concept.

Liouville's theorem

The Liouville's theorem plays a fundamental role in the statistical description of Hamiltonian mechanical systems. It describes the evolution of a density over phase space (relying on the conservation of the volume form under Hamiltonian flow1 ). This density can be seen either as a probability density or as a density of independent 1 Introduction and preliminaries particles. Let us first prove the conservation of the volume form; then we shall explicitly state the theorem.

The natural volume form in the phase space is the top wedge power of the symplectic 2-form: .3.40) where • can be either the floor ( • ), or the ceil ( • ) function. The variation under Hamiltonian flow 1 of this 2N -form is

dV ∝ ω ∧N = 2 N N N/2 dq 1 ∧ dp 1 ∧ dq 2 ∧ dp 2 ∧ • • • ∧ dq N ∧ dp N . ( 1
L χ H dV ∝ L χ H ω ∧N = N (L χ H ω) ∧ ω ∧N -1 = 0 , (1.3.41)
where we have used the invariance of the symplectic 2-from (see eq. (1.3.16)), as well as the fact that the Lie derivative respects the Liebnitz rule with respect to the wedge product. Hence, the natural volume form (1.3.40) is conserved under Hamiltonian flow. We are now going to use this fact to prove the Liouville's theorem.

Let ρ q i , p i , t be a density on the phase space: whether it is a probability density or a density of independent particles. The number of particles/the probability to be, at time t, in the infinitesimal volume dV, centred in q i , p i , is N := ρ q i , p i , t dV.

The time variation of this quantity is only due to the Hamiltonian flow (no disappearing nor creation of particles/probability), hence one successively obtains

∂N ∂t = -L χ H (N ) (1.3.42) = -L χ H (ρ dV) = -(L χ H ρ) dV -ρ L χ H dV = -(ι χ H dρ) dV = -dρ(χ H ) dV = -ω(χ ρ , χ H ) dV = ω(χ H , χ ρ ) dV = {H, ρ} dV , (1.3.43)
where we have used the conservation of the volume form, eq. ( This is Liouville's theorem.

Numerical integration of classical conservative field theories

Time-dependent Hamiltonian

The framework introduced until now assumes that the Lagrangian (as well as the Hamiltonian) does not explicitly depend on time. To handle its possible dependence, the objects introduced so far require some extensions. Let us consider a mechanical system described by a Lagrangian L q i , qi , t that explicitly contains the time. The Hamiltonian is still defined by the Legendre transformation

H q i , p i , t := p i qi -L q i , qi , t , (1.3
} :≡ {t} ∪ q i ∪ {H} ∪ {p i }.
Here, the ordering in the Darboux coordinates is not arbitrary: p i is the conjugate momentum to q i and the same holds for H which is the conjugate momentum of t. Indeed, H is the conserved quantity that correspond to the translational invariance along t as it is for the usual conjugate momentums.

The symplectic 2-form requires some extensions too. It is now defined as ω = ω ab d a ∧ d b := 2 dq i ∧ dp idt ∧ dH , and is naturally closed (dω = 0) and non-degenerated (det ω = 0). Hence, (Ω, ω) is still a symplectic manifold.

The state vector, z ∈ T Ω, becomes

z = z a ∂ a := t q 1 • • • q N H p 1 • • • p N t . (1.3.48)
The equation of motion henceforth reads

ω( ż, •) = 0 , (1.3.49) since ω( ż, •) = Ḣ dt -ṗi dq i + qi dp i -dH = Ḣ dt -ṗi dq i + qi dp i - ∂H ∂t dt - ∂H ∂q i dq i - ∂H ∂p i dp i = 0 .
Hence, the evolution of the state vector is explicitly given by

ż = χ H + ω (dH, •) , (1.3.50)
where χ H is defined as earlier but using the actual extended symplectic form.

In this way, the framework previously introduced can be extended to take into account an explicit time dependence in the Hamiltonian, or in the deformations of phase space.

Throughout this section, we have recalled some elementary constructions of the Hamiltonian mechanics. Later in the chapter (section 2) we shall introduce the De Donder -Weyl framework that generalise these concepts to the field theory. Before that, let us review the numerical integration methods for conservative mechanics.

1 Introduction and preliminaries

Numerical integration of Hamiltonian mechanical systems

The methods for integrating field theories are widely inspired from the one of mechanics, and most of the observations valid in the latter case will remain valid.

The aim of this section is to give an overview of the existing methods for numerically solving conservative ordinary differential equations (odes); keeping in mind our final goal, we shall carefully analyse if the methods can correctly handle the long-time dynamics of the problem.

There is a broad variety of methods to tackle odes numerically, but, as already stated, here we focus on finite-difference methods. Besides, this presentation is not intended to be exhaustive; our objective is to bring out some general concepts, and we only introduce a restricted selection of the most known technics. This presentation is inspired from [START_REF] Collatz | The Numerical Treatment of Differential Equations, Grundlehren der mathematischen Wissenschaften[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential -Algebraic Problems[END_REF][START_REF] Stoer | Introduction to Numerical Analysis[END_REF] where one can find a complete review on numerical methods in general and, in particular, on the subject of numerical integration of odes.

First, two mechanical systems are introduced to be used as application examples. Secondly, we present the construction, as well as a proof of the main properties of some of the most renowned methods. Next, their strengths and weaknesses will be analysed, as well as their accuracy, on the two cases mentioned just above.

The harmonic oscillator and the simple pendulum

For the purpose of evaluating the performances of the numerical integrators, we consider two Hamiltonian mechanical problems whose exact solution is known.

We consider a particle with only one degree of freedom subject to a time-independent Hamiltonian flow. The dynamics of the system is therefore fully described by the state vector

z := q p , (1.4.1a)
subject to the equation of motion

ω ab żb = ∂ a H , (1.4.1b)
where the symplectic form ω = ω ab d a ∧ d b is defined as

ω ab :≡ 0 1 -1 0 . (1.4.1c) i.
The first considered problem is the linear case of the harmonic oscillator whose Hamiltonian is

H := p 2 2 + q 2 2 . (1.4.2)
The equation of motion is of course q = p , and (1.4.3a) ṗ = -q .

(1.4.3b)

Assuming the initial configuration q(t = 0) = q 0 , and (1.4.4a)

p(t = 0) = p 0 , (1.4.4b)
the solution reads q(t) = q 0 cos t + p 0 sin t , and (1.4.5a) p(t) = p 0 cos tq 0 sin t , (1.4.5b)

and is obviously T := 2π-periodic.

ii. The second problem of interest is the so-called simple pendulum. It is a nonlinear exactly solvable model [START_REF] Baker | The Pendulum: A Case Study in Physics[END_REF][START_REF] Beléndez | Exact solution for the nonlinear pendulum[END_REF][START_REF] Suzuki | Physics of Simple Pendulum: A Case Study of Non-Linear Dynamics[END_REF] whose dynamics follows the flow of the Hamiltonian

H := p 2 2 + 1 -cos q . (1.4.6)
The equation of motion is thus q = p , and (1.4.7a) ṗ =sin q .

(1.4.7b)

Let us define the modulus 1 , k, as

k 2 := H 2 = p 0 2 4 + 1 -cos q 0 2 , ie (1.4.8a) k = ± p 0 2 4 + sin 2 q 0 2 . (1.4.8b)
The system can present two sorts of behaviour. First, for k 2 ≤ 1 the system is oscillating, and the exact evolution for such an initial configuration follows q(t) = q 0 + 2 sgn sn t + τ k On the other hand, for k 2 ≥ 1 the system is whirling 2 :

q(t) = q 0 + 2 sgn p 0 am k(t + τ ) 1 k 2 -am kτ 1 k 2
, and (1.4.10a)

p(t) = 2k dn k(t + τ ) 1 k 2 , where (1.4.10b) τ = 1 k dn -1 p 0 2k 1 k 2 . (1.4.10c)
In both cases, the motion has a periodicity 3 .4.11) that depends on the orbit. Here, cn (z|m), dn (z|m) and sn (z|m) are the three principal Jacobi elliptic functions, am (z|m) is the amplitude function and K(m) is the complete elliptic integral of the first kind [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF].

T (q 0 , p 0 ) := 4 K k 2 = p 0 2 4 + sin 2 q 0 2 , ( 1
1 an integral of motion.

2 note that in the following expression sgn p0 cannot vanishes since p0 = 0 necessarly corresponds to the oscillating behaviour (k 2 ≤ 1). 3 since the phase space is compact in the direction of the first Darboux coordinate (q).

Introduction and preliminaries

The Euler's methods

This section describes the construction and properties of the Euler's methods. They are a collection of possible ways to approximate the derivative of a function and consequently an ode.

The first step of the method is to discretise the time as a uniform uni-dimensional lattice that samples the state vector:

t → t n := n δ , and (1.4.12a)

z(t) → z n := z(t n ) . (1.4.12b)
The continuous vector field z(t) is substituted by the samples z n on each of the lattice's nodes t n . This representation reduces the amount of unknowns from a non-enumerable set to a countably infinite. The vector field can then be approximately reconstructed from the samples using an interpolation rule.

The Euler's methods approximate the derivative of a function at point t n+α := (n + α)δ by the finite difference 1 żn+α ≈: .4.13) This approximation differs from the exact definition of a derivative only by the lack of the lim δ→0 ; however, this is precisely this limit that cannot be taken on a computer and that explains the errors produced by numerical integrators.

z n+1 -z n δ . ( 1
In this definition the parameter α can be freely chosen but is generally bounded into [0, 1]. Nevertheless, there are a few noteworthy choices for this parameter:

i. α = 0 is the so called Euler explicit rule, ii. α = 1 /2 is the Midpoint rule, iii. while α = 1 is the Euler implicit rule. Using approximation (1.4.13) the discrete equivalent of the equation of motion (1.4.1b) characterises the dynamics of the discretised problem:

ω ab z b n+1 -z b n δ = ∂ a H(z n+α ) , (1.4.14)
where

z n+α := α z n+1 + (1 -α) z n . (1.4.15)
Thus, the discrete evolution of the harmonic oscillator is given by

q n+1 -q n δ = α p n+1 + (1 -α) p n , and (1.4.16a) p n+1 -p n δ = -α q n+1 -(1 -α) q n , (1.4.16b)
1 note that the centred Euler's rule, żn ≈ zn+1 -zn-1 2δ , produces numerical errors of higher order in δ but leads to an incorrect approximation of the second order derivative:

zn ≈ zn+2 -2zn + zn-2 4δ 2 ,
which consists of two independent sub-lattices (the odd one and the even one).

while the evolution of the simple pendulum by q n+1q n δ = α p n+1 + (1α) p n , and (1.4.17a)

p n+1 -p n δ = -sin (α q n+1 + (1 -α) q n ) . (1.4.17b)
These equations define a map from z n to z n+1 that allows one to obtain an approximate solution of the problem. For α = 0 this application is explicitly defined while in general 1 , for α = 0, the map is implicit. Solving such an implicit map requires finding the root of a multi-dimensional vector function 2 . There exists a broad variety of algorithms for roots finding [START_REF] Madsen | Methods for non-linear least squares problems[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF] but from first principles they are all based on the Newton -Raphson method or the Gauss -Seidel one 3 ; however, these algorithms suffer from their lack of robustness. A more robust approach 4 is to substitute the problem of finding a root by the problem of finding a minimum of the square. This can be achieved using the standard gradient descendent, gradient conjugate or Newton methods; nevertheless, to obtain the numerical results that will be presented later in the section, we have used the Levenberg -Marquardt algorithm that is slightly more sophisticated.

The symplectic form ω n is defined at time t n as follows and should have the same structure at later times if the evolution map preserves it:

ω n = ω ab d a n ∧ d b n , ie (1.4.18a
)

ω n+1 = ω ab d a n+1 ∧ d b n+1 . (1.4.18b)
By rewriting eq. (1.4.14), and applying the exterior derivative on it, one obtains

ω ab z b n+1 = ω ab z b n + δ ∂ a H(z n+α )
, and then (1.4.19a)

ω ab d b n+1 = ω ab d b n + δ α ∂ a ∂ b H d b n+1 + (1 -α) ∂ a ∂ b H d b n , (1.4.19b)
where ∂ a ∂ b H is taken at z n+α . Hence, one successively obtains

[ω ab -δ α ∂ a ∂ b H] d b n+1 = [ω ab + δ (1 -α) ∂ a ∂ b H] d b n , ie (1.4.20) [ω ab -δ α ∂ a ∂ b H] d a n+1 ∧ d b n+1 = [ω ab + δ (1 -α) ∂ a ∂ b H] d a n+1 ∧ d b n ω ab d a n+1 ∧ d b n+1 = ω n+1 = = [ω ab + δ α ∂ a ∂ b H] d a n+1 ∧ d b n + δ (1 -2α) ∂ a ∂ b H d a n+1 ∧ d b n = -[ω ba -δ α ∂ b ∂ a H] d a n+1 ∧ d b n + δ (1 -2α) ∂ a ∂ b H d a n+1 ∧ d b n = -[ω ba + δ (1 -α) ∂ b ∂ a H] d a n ∧ d b n + δ (1 -2α) ∂ a ∂ b H d a n+1 ∧ d b n
1 for non-linear problems.

2 ie from R n to R m .

3 the line search method (ie dichotomy) is not appliable in dimension higher than one. 4 it is generally not recommended to search for a root by a minimization approach since optimisation algorithms may converge to a secondary minimum that actually does not correspond to a root of the function. However, in the present case and if δ is small enough, the minimization algorithm starts in a close vicinity of the root; preventing it from falling into the valley of a secondary minimum.

1 Introduction and preliminaries

= ω ab d a n ∧ d b n + δ (1 -2α) ∂ a ∂ b H d a n+1 ∧ d b n ω n+1 = ω n + δ (1 -2α) ∂ a ∂ b H d a n+1 ∧ d b n . (1.4.21)
The symplectic form is preserved by the evolution map if and only if i. δ = 0, which is the continuous limit and has no practical interest, ii. or α = 1 /2, which is the Midpoint rule, iii. or

∂ a ∂ b H d a n+1 ∧ d b n = 0 . (1.4.22)
Using eq. (1.4.20), d a n+1 can be explicitly expressed by inverting the term in brackets on the left-hand side:

d a n+1 = D a b d b n , (1.4.23)
where D a b is a complicated expression but we do not need it explicitly. Thus, condition (1.4.22) becomes

D c a ∂ c ∂ b H d a n ∧ d b n = 0 , (1.4.24)
and since d a n ∧ d b n is a basis for the 2-forms, this is equivalent to

D c a ∂ c ∂ b H = 0 . (1.4.25)
This equation is satisfied provided that i. H is a solution of this differential equation but this is not true in general, ii. or ∂ c ∂ b H = 0 which is a particular case too, iii. or finally, if D c a = 0 which implies d a n+1 = 0 (ie the evolution map degenerates the symplectic form and is therefore not able to preserve it).

Hence, the only general solution is α = 1 /2. Since it exactly preserves the symplectic form over the phase space, the Midpoint rule is called a symplectic integrator; meaning that the evolution map defined by the Midpoint rule is a symplectic transformation of the phase space.

As we shall highlight later in this section, the conservation of the symplectic form1 implies excellent conservation properties of the integrals of motion, even in the longtime dynamics (ie after a large number of time-steps). Symplectic integrators are the only viable way to correctly perform the integration of a Hamiltonian mechanical system over large times.

Let us now investigate the order of the error made due to the discretisation of the equation of motion. Starting from eq. (1.4.1b), one successively obtains

ω ab żb = ∂ a H , ie ω ac ω cb żb = ω ab ∂ b H -δ a b żb = ω ab ∂ b H ża = -ω ab ∂ b H =: χ a , then (1.4.26a) za = żb ∂ b χ a = χ b ∂ b χ a , and (1.4.26b) ... z a = żc ∂ c χ b ∂ b χ a + żc χ b ∂ c ∂ b χ a = χ c ∂ c χ b ∂ b χ a + χ c χ b ∂ c ∂ b χ a . (1.4.26c)
On the other hand, the exact evolution within a time-step δ can be obtained through the Taylor series as

Z a n+1 = z a n + δ ża (z n ) + δ 2 2 za (z n ) + δ 3 6 ... z a (z n ) + O δ 4 = z a n + δ χ a (z n ) + δ 2 2 χ b (z n ) ∂ b χ a (z n ) + δ 3 6 χ c (z n ) ∂ c χ b (z n ) ∂ b χ a (z n ) + χ c (z n ) χ b (z n ) ∂ c ∂ b χ a (z n ) + O δ 4 .
Then, noticing that z n+1z n = O(δ), and using eq. (1.4.14), one successively obtains

z a n+1 = z a n + δ χ a (z n + α (z n+1 -z n )) = z a n + δ χ a (z n ) + δ α z b n+1 -z b n ∂ b χ a (z n ) + δ α 2 2 z c n+1 -z c n z b n+1 -z b n ∂ c ∂ b χ a (z n ) + O δ 4 = z a n + δ χ a (z n ) + δ 2 α χ b (z n + α (z n+1 -z n )) ∂ b χ a (z n ) + δ 3 α 2 2 χ c (z n + α (z n+1 -z n )) χ b (z n + α (z n+1 -z n )) ∂ c ∂ b χ a (z n ) + O δ 4 = z a n + δ χ a (z n ) + δ 2 α χ b (z n ) + α z c n+1 -z c n ∂ c χ b (z n ) + O δ 2 ∂ b χ a (z n ) + δ 3 α 2 2 χ c (z n ) + O(δ) χ b (z n ) + O(δ) ∂ c ∂ b χ a (z n ) + O δ 4 = z a n + δ χ a (z n ) + δ 2 α χ b (z n ) ∂ b χ a (z n ) + δ 3 α 2 χ c (z n + α (z n+1 -z n )) ∂ c χ b (z n ) ∂ b χ a (z n ) + δ 3 α 2 2 χ c (z n ) χ b (z n ) ∂ c ∂ b χ a (z n ) + O δ 4 = z a n + δ χ a (z n ) + δ 2 α χ b (z n ) ∂ b χ a (z n ) + δ 3 α 2 [χ c (z n ) + O(δ)] ∂ c χ b (z n ) ∂ b χ a (z n ) + δ 3 α 2 2 χ c (z n ) χ b (z n ) ∂ c ∂ b χ a (z n ) + O δ 4 = z a n + δ χ a (z n ) + δ 2 α χ b (z n ) ∂ b χ a (z n ) + δ 3 α 2 χ c (z n ) ∂ c χ b (z n ) ∂ b χ a (z n ) + 1 2 χ c (z n ) χ b (z n ) ∂ c ∂ b χ a (z n ) + O δ 4 .
The error committed because of the discretisation of the equation of motion is therefore

z a n+1 -Z a n+1 = δ 2 α - 1 2 χ b (z n ) ∂ b χ a (z n ) + δ 3 α 2 - 1 6 χ c (z n ) ∂ c χ b (z n ) ∂ b χ a (z n ) + α 2 2 - 1 6 χ c (z n ) χ b (z n ) ∂ c ∂ b χ a (z n ) + O δ 4 , (1.4.27)
and since, for any α, the term of order δ 3 cannot be cancelled in general, one has

z a n+1 -Z a n+1 = δ 2 α - 1 2 χ b (z n ) ∂ b χ a (z n ) + O δ 3 . (1.4.28)
Thus, the Midpoint rule is a second-order integrator (z a n+1 -Z a n+1 = O δ 3 ) while all the other choices of α lead to first-order integrators (z a n+1 -Z a n+1 = O δ 2 ).

Introduction and preliminaries

Partitioned Euler

The partitioned Euler method is a variation of the method introduced above that uses the possible separability of the Hamiltonian. Here, for simplicity, we consider a system with only one degree of freedom and we also assume 1 that the dynamics are characterised by a separable Hamiltonian:

H(q, p) := T (p) + V (q) . (1.4.29)
The aim of a partitioned method is to separate the state vector into several sub-parts and then, to discretise each part on a different lattice, with different rules. In the present situation the state vector

z := q p , (1.4.30)
is partitioned into the two parts: q and p. Next, q(t) and p(t) are both discretised on a uniform uni-dimensional lattice but with different origins; they are shifted by a half-spacing:

t → t n := n δ , (1.4.31a) q(t) → q n := q(t n ) , and (1.4.31b) p(t) → p n := p t n + δ /2 . (1.4.31c)
Afterwards, we use the explicit Euler rule for q and the implicit one for p

q n+1 -q n δ = ∂T ∂p (p n ) , and (1.4.32a)
p n+1 -p n δ = - ∂V ∂q (q n+1 ) . (1.4.32b)
The resulting approximation is somehow a Midpoint rule in the sense that the finitedifference represents the derivative at the midway of the time-step. However, the interest of having separated the state vector results in the fact that the evolution map is not implicit (by contrast to the Midpoint approximation). Indeed, solving first the evolution of q and then the evolution of p, the application is explicit:

q n+1 = q n + δ ∂T ∂p (p n ) , and (1.4.33a
)

p n+1 = p n -δ ∂V ∂q (q n+1 ) = p n -δ ∂V ∂q q n + δ ∂T ∂p (p n ) . (1.4.33b)
Let us now have a look at the same questions as for the Euler methods: is the evolution map a symplectic transformation and what is the order of the integrator?

First, the symplectic form on the lattice is given by

ω n = 2 dq n ∧ dp n-1 /2 = 2 dq n ∧ dp n + dp n-1 2 = dq n ∧ (dp n + dp n-1 ) , and (1.4.34a) ω n-1 /2 = 2 dq n-1 /2 ∧ dp n-1 = 2 dq n + dq n-1 2 ∧ dp n-1 = (dq n + dq n-1 ) ∧ dp n-1 .
(1.4.34b)

1
this assumption is required; in contrast to the number of degree of freedom that we have reduced to the minimum in the aim to simplify the presentation.

By taking the exterior derivative on eqs. (1.4.33a) and (1.4.33b), we obtain the evolution equations for the forms

dq n+1 = dq n + δ ∂ 2 T ∂p 2 (p n ) dp n , and (1.4.35a
)

dp n+1 = dp n -δ ∂ 2 V ∂q 2 (q n+1 ) dq n+1 . (1.4.35b)
Then, one successively obtains

ω n+1 = dq n+1 ∧ (dp n+1 + dp n ) = dq n+1 ∧ dp n -δ ∂ 2 V ∂q 2 (q n+1 ) dq n+1 + dq n+1 ∧ dp n = 2 dq n+1 ∧ dp n = dq n+1 ∧ dp n + dq n + δ ∂ 2 T ∂p 2 (p n ) dp n ∧ dp n = (dq n+1 + dq n ) ∧ dp n = ω n+ 1 /2 (1.4.36) = dq n + δ ∂ 2 T ∂p 2 (p n ) dp n ∧ dp n + dq n ∧ dp n = 2 dq n ∧ dp n = dq n ∧ dp n + dq n ∧ dp n-1 -δ ∂ 2 V ∂q 2 (q n ) dq n = dq n ∧ (dp n + dp n-1 ) = ω n (1.4.37) = dq n ∧ dp n-1 -δ ∂ 2 V ∂q 2 (q n ) dq n + dq n ∧ dp n-1 = 2 dq n ∧ dp n-1 = dq n ∧ dp n-1 + dq n-1 + δ ∂ 2 T ∂p 2 (p n-1 ) dp n-1 ∧ dp n-1 = (dq n + dq n-1 ) ∧ dp n-1 = ω n-1 /2 . (1.4.38)
Hence,

ω n+1 = ω n+ 1 /2 = ω n = ω n-1 /2 . (1.4.39)
Therefore, the partitioned Euler approximation preserves the symplectic form. In fact, each sub-step is a symplectic transformation: the evolution map for q as well as the one for p independently preserve the symplectic form.

Let us now determine the order of the errors resulting from the discretisation process.

By Taylor expansion, we first obtain the exact evolution within a time-step δ

Q n+1 = Q n+1 + q n 2 + δ 2 qn+ 1 /2 + δ 2 8 qn+ 1 /2 + O δ 3 = q n + δ ∂T ∂p (p n ) - δ 2 4 ∂V ∂q q n+ 1 /2 ∂ 2 T ∂p 2 (p n ) + O δ 3
, and

P n+1 = P n+1 + p n 2 + δ 2 ṗn+ 1 /2 + δ 2 8 pn+ 1 /2 + O δ 3 = p n -δ ∂V ∂q (q n+1 ) - δ 2 4 ∂T ∂p p n+ 1 /2 ∂ 2 V ∂q 2 (q n+1 ) + O δ 3 .

Introduction and preliminaries

Recalling the equation of motion that governs the actual evolution of the system,

q n+1 = q n + δ ∂T ∂p (p n ) ,
and

p n+1 = p n -δ ∂V ∂q (q n+1 ) ,
we can estimate the errors produced by the numerical integrator at each time-step: and (1.4.40a)

q n+1 -Q n+1 = δ 2 4 ∂V ∂q q n+ 1 /2 ∂ 2 T ∂p 2 (p n ) + O δ 3 = O δ 2 ,
p n+1 -P n+1 = δ 2 4 ∂T ∂p p n+ 1 /2 ∂ 2 V ∂q 2 (q n+1 ) + O δ 3 = O δ 2 , (1.4.40b)
since the terms of O δ 2 do not vanish in general. Hence, the partitioned Euler method is a first-order (explicit) symplectic integrator.

The numerical behaviour of these integrators will be discussed below in section 1.4 page 36. In particular, we shall observe the importance of being symplectic or not.

We shall also observe that the higher the order of the integrator, the faster the errors will vanish when δ goes to zero (as a power law that corresponds to the order of the integrator). However, since our interest is in reaching long-time, it inhibits having a too small δ, and the order of the approximation will actually not be the most important deciding factor.

A selection of other well-known methods

As mentioned earlier, there exist many numerical methods to integrate odes [START_REF] Collatz | The Numerical Treatment of Differential Equations, Grundlehren der mathematischen Wissenschaften[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential -Algebraic Problems[END_REF][START_REF] Stoer | Introduction to Numerical Analysis[END_REF]. In the present section, we shall discuss two noteworthy classes of methods that can be constructed in such a way to be symplectic. On the one hand, the methods of the Runge -Kutta type and, on the other hand, the methods constructed from the Yoshida expansion of the evolution operator [START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF]. The aim of this section if only to give a sketch of these methods; it is absolutely not intended to be exhaustive. There are two reasons for that: firstly, because the class of the Runge -Kutta integrators is a wide subject, whose discussion should require an entire book. Secondly, because these two kinds of methods are not suitable for a generalization to Hamiltonian pdes.

In fact they can be generalised to Hamiltonian pdes: many papers treat of Runge -Kutta high-order multi-symplectic integrators, and it would be possible to generalise the Yoshida expansion as well. However, this generalization is at the cost of treating space and time in a different way, manifestly breaking the covariance of the theory. The importance of the covariance will be presented in the next sections but we shall carefully pay attention, and put a lot of efforts not to break this symmetry. Hence, treating space and time differently is not a satisfying discretisation approach.

Here we introduce these two classes of numerical approximations because they are notorious, and have interesting properties when applied on Hamiltonian odes. As earlier, the state vector is sampled on the nodes of the uniform uni-dimensional lattice that discretises time:

t → t n := n δ , and (1.4.42a) z(t) → z n := z(t n ) . (1.4.42b)
The exact evolution of the system in phase space is obtained by moving along the integral curves of the Hamiltonian vector field, thus

z a n+α = z a n + δ α 0 dβ χ a (z n+β ) , (1.4.43)
where α ∈ [0, 1]. The idea of the Runge -Kutta methods is to approximate this integral by a quadrature formula: the time interval, δ, is decomposed in M sub-intervals, and then the integral is approximated by the finite sum

z a n+ m /M = z a n + δ M M k=0 R mk χ a z n+ k /M , (1.4.44)
where m ∈ 1, M . This approximated evolution map is explicit as soon as R mk = 0 for all k ≥ m. The matrix R is the approximation pattern, and is chosen such that the first terms of the Taylor expansion vanish (up to O δ M +1 for an explicit integrator). For instance, explicit and implicit Euler methods are Runge -Kutta integrators of order M = 1 with, respectively,

R = 0 0 1 0 , and R = 0 0 0 1 . (1.4.45)
The choice of the coefficients of R has been widely studied, and it is possible to define high-order symplectic integrators in this way. We are not going to detail how to determine the R matrix, but we refer to the broad literature on the subject.

ii. The second class of methods introduced in this section is based on the Yoshida expansion of the evolution operator. These methods are based on a splitting of the Hamiltonian, and it is required to assume it as separable:

H(q, p) := T (p) + V (q) . (1.4.46)
The equation of motion (1.4.1b) can be rewritten in terms of Poisson bracket as q = {q, H} , and (1.4.47a)

ṗ = {p, H} . (1.4.47b)
The solutions of these equations can be abstractly written as q(t) = U (tt 0 ) q(t 0 ) , and (1.4.48a)

p(t) = U (t -t 0 ) p(t 0 ) , (1.4.48b)
where U is the evolution operator of the system:

U (t -t 0 ) := e (t -t 0 ) {•, H} = e (t -t 0 ) ({•, T } + {•, V }) . (1.4.49)
1

Introduction and preliminaries

For small enough evolution time, this evolution operator can be decomposed as

U (δ) = e δ {•, H} = e δ ({•, T } + {•, V }) = M m=1 e α m δ {•, T } • e β m δ {•, V } + O δ M +1 (1.4.50a) =: Ũ (δ) + O δ M +1 , (1.4.50b)
where the 2M coefficients {α m , β m } can be found using the Baker -Campbell -Hausdorff formula [START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF]. Ũ (δ) is, therefore, an approximate evolution operator of order M . Its interest is that it is explicit and preserves the symplectic structure of the phase space. Let us now prove this. We first notice that

{•, T } • 2 := {{•, T }, T } = 0 , and (1.4.51a) {•, V } • 2 := {{•, V }, V } = 0 . (1.4.51b) since {q, T } • 2 = {{q, T (p)}, T (p)} = ∂T ∂p (p), T (p) = 0 , {p, T } • 2 = {{p, T (p)}, T (p)} = {0, T (p)} = 0 , {q, V } • 2 = {{q, V (q)}, V (q)} = {0, V (q)} = 0 , and {p, V } • 2 = {{p, V (q)}, V (q)} = - ∂V ∂q (q), V (q) = 0 .
Hence, Ũ (δ) can be explicitly be written using

e α m δ {•, T } = id + α m δ {•, T } , and (1.4.52a) e β m δ {•, V } = id + β m δ {•, V } . (1.4.52b)
These operators act as transformations of the phase space: the first one acts as

Q(q, p) := e α m δ {•, T } q = q + α m δ ∂T ∂p (p) , and (1.4.53a)
P (q, p) := e α m δ {•, T } p = p , (1.4.53b)
and preserves the Poisson bracket:

{Q, P } qp = {q, p} qp + α m δ ∂T ∂p (p), p qp = {q, p} qp . (1.4.53c)
The second one acts as Q(q, p) := e β m δ {•, V } q = q , and (1.4.54a)

P (q, p) := e β m δ {•, V } p = p -β m δ ∂V ∂q (q) , (1.4.54b)
and therefore preserves the symplecticity of the phase space:

{Q, P } qp = {q, p} qp -β m δ q, ∂V ∂q (q) qp = {q, p} qp . (1.4.54c)
Finally, since Ũ (δ) is a composition of symplectic transformation, it is an (explicit, M -order) symplectic integrator for H.

Numerical integration of classical conservative field theories

Numerical results

In the present section we discuss the numerical results obtained using the four following methods:

i. the explicit Euler approximation as described in section 1.4 with α = 0, ii. the Midpoint rule (section 1.4, α = 1 /2), iii. the implicit Euler method (section 1.4, α = 1), iv. and finally, the partitioned Euler approach (see section 1.4). Two of them are symplectic (the Midpoint rule as well as the partitioned Euler method); the others are not. Two of them are explicit (the explicit and the partitioned Euler methods) while the two others are implicit. We thus have the four possible combinations of symplecticity and to be explicit or not. All the integrators are of the first order in δ, except the Midpoint rule which is second order. The integrators will be put to the test on two mechanical problems: on the one hand, the harmonic oscillator (a linear theory) and, on the other hand, the simple pendulum (a non-linear oscillator).

Once again, there is nothing original in this discussion; we do the same observations as described in the literature.

In all the following figures we use the symbolic operator ∆: it is an error estimator of the equation it is applied on, and has no dimension; it must be compared to 1.

The numerical evolution is presented in fig. I.2 for the harmonic oscillator and in fig. I.3 for the simple pendulum. First, we observe that all the methods have a tendency to produce a slower dynamics than the expected one (none of them is able to preserve the period of the motion). Secondly, we observe that the explicit and implicit Euler methods (the first row) manifestly do not preserve the volume of phase space; the geometry of the phase space is strongly broken. The explicit approximation amplifies the energy, leading to a divergence, while the implicit one decreases the energy of the system. The two symplectic approximations (Midpoint and partitioned Euler) preserve the volume in a much better way1 , even if they produce an unexpected deformation: they preserve the geometrical properties of the phase space but are still an approximation of the flow. The Midpoint rule follows the exact evolution in a closer way than the partitioned Euler approximation; we cannot attribute this observation to the fact they have a different approximation order since for such a δ the errors they commit are of the same magnitude (see fig. I.4 for δ = T /8). Hence, we assume (but we shall argue this later) that the implicit aspect of the approximation improves its behaviour. Moreover, on the phase space of the simple pendulum (fig. I.3) we observe that the two implicit methods (the right column) respect the Hamiltonian flow lines while the explicit methods (left column) produce a significant distortion of the latter; observe in particular the vicinity of the point q = 3π /2, p = 2 . Figures I.4 and I.5 represent the error committed by the different methods as a function of the time-step. The first remark we can make (looking at the graph showing the error committed on the solution) is that the methods present the expected convergence when δ goes to zero: all the methods are of the first-order, except the Midpoint rule which is a second-order integrator. Secondly, we observe that the explicit and implicit Euler methods behave in the same way except for large δ; the explicit approximation diverges, producing unbounded errors, while the implicit one rapidly loses all the initial energy,

Introduction and preliminaries

Explicit Euler

Implicit Euler

-1 0 1 2 -2 -1 0 1 q p Partitioned Euler Midpoint rule Figure I.
2 -Deformation of the phase space under the approximate evolution map for the four methods considered here when applied to the harmonic oscillator problem. The left column regroups the explicit integrators while the ones on the right column are implicit; the integrators on the second row are symplectic (those on the first row are not). The evolution flow is represented in light grey and the volumes drawn using the same colour represent the exact deformation of the phase space under the Hamiltonian flow. By volumes we point out the regions of phase space, initially drawn as squares, but that may be distorted by the integration process.

The filled volumes represent the successive deformation of the phase space obtained using the numerical integrators; a single time-step of length δ = T /8 is used to pass from a volume to the next one. We proceed as follow: we first discretise the volume in the phase space, then we independently perform a time-step of length δ = T /8 for each state constituting the initial volume; once done we reconstruct the transformed volume from these updated samples. The differences between the different integrators must not be attributed, here, to the fact they have a different order of approximation; indeed, for such a δ all the integrators produce errors of the same order of magnitude (the higher order of the Midpoint approximation distinguishes it from the other integrators only when δ goes to zero).

Numerical integration of classical conservative field theories

Explicit Euler Implicit Euler

-π 0 π 2π -4 -2 0 2 4 q p
Partitioned Euler Midpoint rule Figure I.3 -Consequences in phase space of numerically approximating the evolution of the simple pendulum using each of the four methods considered in this section. See fig. I.2 for a description of the content; here the process only differs in that the numerical evolution between two successive volumes is obtained using two timesteps of length δ = T (q 0 , p 0 )/16 (where (q 0 , p 0 ) is the centre of the initial square). The effective distance in between two successive volumes is still T (q 0 , p 0 )/8, however performing such an evolution in only one step produces too large deviations with respect to the exact evolution; the comparison becomes impossible, especially for the less precise integrators. and on the energy (bottom graph) as a function of the time-step for the harmonic oscillator problem. The worst error encountered in the course of the numerical integration of the system over a period of the motion (0 ≤ t ≤ T ) is represented here. Both graphs have the same horizontal axis. The exact integrator stands for the exact evolution map recursively applied; hence, the obtained solution is exact in the sense that there are no errors due to the discretisation but it is, however, subject to the accumulation of truncation errors. This exact integrator represents the maximum possible accuracy and is directly related to the precision of the representation employed by the computer to encode real numbers (here double floating-point). I.4 for a description of the content; here the exact integrator is not represented since its behaviour is unchanged with respect to the case of the harmonic oscillator. The exact integrator is not affected by the nonlinearities; only by the truncation errors that have no reason to be different in this situation. Moreover, it is extremely expensive in computation time since it requires to evaluate many Jacobi elliptic functions.

producing finite errors corresponding to the distance in the phase space between the origin (the zero energy point) and the exact solution. Regarding the two symplectic integrators, the same remark holds except that they have a different behaviour in the limit where δ goes to zero since they are of a different approximation order. Hence, we can conclude that the dynamics with a large time-step is most accurately described using implicit methods. In contrast, in the limit of small time-steps, the implicit aspect of the approximation is not decisive. Concerning the exact method, it is exact, and so it is only subject to the truncation errors due to the finite machine-precision (here almost 16 significant digits). These truncation errors increase by accumulation as the number of steps extends and, therefore, as δ decreases. All the above remarks hold for the conservation of the energy as well. However, in addition, there is a notable property: when applied to a linear problem, the Midpoint rule exactly preserves the energy (this property disappears with the non-linearity of the theory). This means that one can control how accurately the Midpoint rule conserves the energy by carefully analysing the non-linearity of the problem; this fact will be proven later.

Let us now discuss fig. I.6: it is a reduced-precision version of fig. I.4, where the computations were done with only 4 significant digits while all the other numerical results presented in this thesis were obtained in double floating-point precision (about 16 significant digits). In practice there are two sources of errors:

i. on the one hand, the errors due to the approximation of the equation of motion, and more precisely, the approximation of the derivative by a finite-difference, ii. on the other hand, the truncation errors due to the fact that real numbers are represented on the computer with a finite-precision.

In normal circumstances, the errors due to item i. decrease with δ. In contrast, the value of δ has almost no influence on the errors described in item ii.. Nonetheless, to reach a given final time, if δ decreases the number of steps increases and, with the latter, the number of times the truncation error will occur. Therefore, for large values of δ the errors almost exclusively come from item i. while, for small values of δ, the most significant errors come from item ii.. Thus, to reach a final given time, there is a threshold beyond which reducing δ will no longer improve the accuracy of the solution (see the 'V' shape of the errors plotted on fig. I.6). The higher the order of the integrator and the lower the precision of the numbers, the faster this threshold is reached.

Figures I.7 and I.8 represent the dynamic evolution of the errors in the course of a long integration process. Let us first observe the instantaneous errors; they have a very different behaviour depending on the symplecticity of the method: for such an integrator, the instantaneous errors are fluctuating (with some kind of periodicity) and are, therefore, not trustworthy. For the integrators that are not symplectic, the instantaneous errors are always the worst encountered ones since the initial time.

Hence, the relevant error estimator is the maximum encountered in the course of the integration process; not the instantaneous one, which might be much smaller than the correct estimator. As a second remark, no method can correctly describe the solution over long time-scales (especially for a non-linear theory). However, the reason depends on the method: the explicit Euler method diverges, the implicit one is subject to an energy loss; both approximations are not trustworthy. The symplectic integrators suffer the same disease but for a different reason: their dynamics are a bit slower than the expected ones (the periodicity of the motion is not exact). Depending on the situation this may be a lesser issue or not. Let us now observe how the energy is conserved: the explicit and implicit Euler methods are absolutely not able to describe Both graphs have the same horizontal axis. Beyond t = T the horizontal axis is shown in a different linear scale and the curves with closed symbols are not plotted since they vary too rapidly with respect to this new timescale (these represent instantaneous errors that in any case are not relevant on this time-scale). In order to fairly compare the different methods, δ is chosen independently for each of them such that they all produce the same error on the solution (10 -2 ) at t = T . In every case, δ is too small with respect to the time-scales considered here; there are only a few points represented on the graph, however, all the steps are taken into account for the worst encountered error.

Numerical integration of classical conservative field theories the evolution of the system since they strongly break the conservation of the energy. Concerning the symplectic integrators, they quite accurately preserve the energy, but the most important remark is that the worst errors occur in the very early evolution of the system: after only a few steps we can estimate the accuracy of the whole integration process; whatever its duration.

10 -5 10 -4 10 -3 10 -2 10 -1 ∆(z = z exact ) ∆(␣) max ≤t ∆(␣)
Considering all the elements presented here, the only viable approach to solve a Hamiltonian ode over long time-scales is a symplectic integrator; preferentially an implicit approximation. We shall show how to generalise this to pdes in section 3.

However, developing a multi-symplectic integrator requires first to introduce the correct framework to deal with Hamiltonian field theories. The aim of the next section, is to introduce a number of necessary concepts of the De Donder -Weyl Hamiltonian formalism.

The De Donder -Weyl Hamiltonian formulation of field theories

This section introduces the concept of multi-symplecticity through the framework of the De Donder -Weyl Hamiltonian formulation of field theories [START_REF] Donder | Théorie invariantive du calcul des variations[END_REF].

We first introduce the De Donder -Weyl formalism and the concept of Hamiltonian pdes. Next, we prove that the DW definition of the phase space is a multi-symplectic manifold, and we establish the conservation of the multi-symplectic structure under Hamiltonian flow. Afterwards, we discuss a possible issue in this fundamental structure and we present a way to address it. Finally, we shall discuss the definition and properties of the stress-energy tensor. This section is treated both in a general setup and with the example of the non-linear wave equation.

Most of the points presented here are just reminders except for two of them. Firstly, the link between the De Donder -Weyl formulation and multi-symplectic geometry is not so common (usually, the De Donder -Weyl formulation is treated through the formalism of the poly-symplectic geometry, see footnote 1 page 48). Secondly, as far as we know, the discussion of the degeneracy of the multi-symplectic structure, and especially its resolution in any dimension, is completely new. In fact, this issue was already discussed in [START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian {PDEs} that conserve symplecticity[END_REF] for the particular case of the non-linear wave equation in 1 + 1 dimensions. However, the argument they propose in this paper does not enable any generalisation to higher dimensions; the new argument we introduce here naturally extends in any dimension.

The De Donder -Weyl (DW) Hamiltonian formalism is a broad topic that is still not clearly understood; once again this presentation is not intended to be exhaustive. See [START_REF] Atanasiu | Distinguished torsion, curvature and deflection tensors in the multi-time Hamilton geometry[END_REF][START_REF] Donder | Théorie invariantive du calcul des variations[END_REF][START_REF] Echeverría-Enríquez | Multivector field formulation of Hamiltonian field theories: equations and symmetries[END_REF][START_REF] Forger | A general construction of poisson brackets on exact multusymplectic manifolds[END_REF][START_REF] Giachetta | Advanced classical field theory[END_REF][START_REF] Gotay | A multisymplectic framework for classical field theory and the calculus of variations II: Space+ time decomposition[END_REF][START_REF] Hélein | Covariant Hamiltonian formalism for the calculus of variations with several variables: Lepage-Dedecker versus De Donder-Weyl[END_REF][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]Basic structures of the covariant canonical formalism for fields based on the De Donder-Weyl theory[END_REF][START_REF]From the Poincaré-Cartan form to a Gerstenhaber algebra of Poisson brackets in field theory[END_REF][START_REF]On field theoretic generalizations of a Poisson algebra[END_REF][START_REF]Canonical structure of classical field theory in the polymomentum phase space[END_REF][START_REF]On the Duffin-Kemmer-Petiau formulation of the covariant Hamiltonian dynamics in field theory[END_REF][START_REF]On a generalization of the Dirac bracket in the De Donder-Weyl Hamiltonian formalism[END_REF][START_REF] Mclean | Covariant field theory on frame bundles of fibered manifolds[END_REF][START_REF] Paufler | De Donder-Weyl equations and multisymplectic geometry[END_REF][START_REF] Struckmeier | Covariant Hamiltonian field theory[END_REF][START_REF] Vankerschaver | The momentum map for nonholonomic field theories with symmetry[END_REF][START_REF] Vankerschaver | Geometric aspects of nonholonomic field theories[END_REF][START_REF] Rieth | The Hamilton-Jacobi theory of De Donder and Weyl applied to some relativistic field theories[END_REF] for further details.

From Lagrangian to DW Hamiltonian formulation

Generic field theory

Let us start with a space-time Lorentzian manifold (a manifold endowed with a Lorentzian metric), M, of dimension D = 1 + d. We assume M to be non-dynamic (ie the metric is not subject to an equation of motion) and flat1 with metric η :≡ diag (1, -1, • • • , -1). We parameterise M by the local coordinate system {x µ }, with ∂ µ := ∂ /∂x µ a basis of T M, and where µ ∈ 0, d . Next, we consider a field theory on M, described by the action S φ i , where φ i is a collection of dynamic fields with i ∈ 1, N . This action originates in a Lagrangian density, L, which is assumed to depend only on the field and its first derivatives:

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

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

The idea of the Hamiltonian formulation of classical mechanics is to substitute the generalised velocity ( 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, they define

ψ i µ := ∂L ∂(∂ µ φ i ) , (2.1.3)
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 , (2.1.4)
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 that parameterise a differential manifold, Ω, which is the De Donder -Weyl definition of phase space (a multisymplectic manifold as we shall prove later). The Hamilton equations generalise to

∂ µ ψ i µ = - ∂H ∂φ i , and (2.1.5a) ∂ µ φ i = ∂H ∂ψ i µ . (2.1.5b)
This is the DW Hamiltonian formulation of a classical field theory.

Equations (2.1.5a) and (2.1.5b) can be rewritten in a more symmetrical way following the same approach as in the case of the classical mechanics. Let us now parameterise the phase space, Ω, with the generalised coordinates {ζ a } :≡ φ i ∪ {ψ i µ } where a is an index conveniently chosen to sweep the collection. Let ∂ a := ∂ /∂ζ a be a basis of T Ω, and {d a := dζ a } be the dual basis in T * Ω. 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 state vector, ζ = ζ a ∂ a , that allows one to rewrite eqs. (2.1.5a) and (2.1.5b) in the abstract form

M µ • ∂ µ ζ = ∇H , (2.1.6)
where {M µ } is a set of constant skew-symmetric matrices of T Ω.

This abstract form, eq. (2.1.6), will be the starting point for the multi-symplectic integrator in the light-cone coordinates (msilcc) we shall develop in section 3. This numerical method is designed to accurately integrate any pde that can be written in such a form, and where the underlying multi-symplectic structure is non-degenerate (see below).

The multi-symplectic structure will directly arise from eq. (2.1.6), and the next step is now to exhibit it 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 example

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

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

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

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

H := ψ µ ∂ µ φ -L = 1 2 ψ µ ψ µ + V (φ) . (2.1.10)
Defining the state vector

ζ := φ ψ 0 ψ 1 • • • ψ d t , (2.1.11a)
the equation of motion is given by eq. (2.1.6) provided that

M µ ab := δ µ a-1 δ 0 b -δ µ b-1 δ 0 a , (2.1.11b)
where a, b ∈ 0, D .

In dimension D = 1 + 0 (d = 0), aliasing q :≡ φ, p :≡ ψ 0 and H :≡ H, we recover the expected Hamilton equations of mechanics

0 -1 1 0 • q ṗ =     ∂H ∂q ∂H ∂p     = V (q) p .
(2.1.12)

In dimension D = 1 + 1 (d = 1) the M matrices read

M 0 :=   0 -1 0 1 0 0 0 0 0   , and
(2.1.13a)

M 1 :=   0 0 -1 0 0 0 1 0 0   , (2.1.13b)
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 shall see in the next section, this fact causes some difficulties. For the moment, we are going to introduce the multi-symplectic structure of Ω. Then, we shall return to the vanishing eigenvalues, and we shall explain how to treat them in the particular case of the non-linear wave equation.

Multi-symplectic structure

The multi-symplectic structure This section will be devoted to the construction of the multi-symplectic structure on Ω, 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.1.6).

Let us first define what a multi-symplectic manifold is. First, we recall that a symplectic manifold is defined as being a differentiable manifold equipped with a closed non-degenerate 2-form. This particular 2-from is called the symplectic form. A possible generalization 1 of this, is the concept of multi-symplectic manifold: it is a differential manifold equipped with several independent symplectic forms (in the present case, it will be one per space-time direction).

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

ω µ := - 1 2 M µ ab d a ∧ d b . (2.2.1a)
They act on pair of vectors as

ω µ ( * , ) := M µ • * , = - * , M µ • , (2.2.1b)
for any * , ∈ T Ω, where ␣, ␣ is the scalar product on T Ω, defined through the metric

η = η ab d a ∨ d b .
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 shall come back to this at the end of this section) that they are all non-degenerate. Finally, since the {M µ } are 1 one can also encounter in the literature the concept of poly-symplectic (or n-plectic) manifold;

it is actually a different object: a n-plectic manifold being a differential manifold with a closed nondegenerate (n + 1)-form (poly-symplectic stands for any n > 1, while symplectic is equivalent to 1-plectic). A multi-symplectic manifold is necessarily poly-symplectic too since the wedge product of D independent 2-forms is a 2D-form. However, the reciprocal is not necessarily true. Finally, we want to stress that the meaning of poly-symplectic and multi-symplectic can be exchanged depending on authors.

linearly independent (as we shall show shortly) 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.1.6), 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 . (2.2.2)
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 multisymplectic only for d = 0 (that is to say, for the mechanical problem).

Degeneracy

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 shall 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 (ie 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 of space-time in fig. I.9. Now, to construct the new field theory explicitly, we first consider the collection of fields

φ ψ 0 ∂ 0 d = 0 d = 1 φ ψ 0 ψ 1 γ ∂ 0 ∂ 1 d = 2 φ ψ 0 ψ 1 ψ 2 γ 01 γ 02 γ 12 Γ ∂ 0 ∂ 1 ∂ 2
Φ (i) µ 1 •••µ i := Φ (i) [µ 1 •••µ i ] µ 1 ,••• ,µ i ∈ 0,d i∈ 0,D , (2.2.3)
where the pair of square brackets denote the anti-symmetrisation defined in section 1.2.

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.2.3) 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 , (2.2.4a)
where the ν indices are in 0, d , which is nothing else than the concatenation of 2 d independent theories 1 (note that all these theories need to belong in the class of the non-linear 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) -Φ (2j) µν 1 •••ν 2j-1 ∂ µ ∂ ν Φ νν 1 •••ν 2j-1 (2j)
, 1 the anti-symmetrisation of the indices in the potential is required since the fields are antisymmetric objects, and, in particular, V 01

(2) Φ

(2)

01

should coincide with V 10

(2) Φ

(2)

10 .

and then, using an integration by parts (and ignoring the boundary terms1 ), we claim that

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

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 (2.2.4c)
and this term does not affect the dynamics either.

We now consider the theory described by the Lagrangian density

L := L 0 + L + + L -, (2.2.5)
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.2.3), and for each direction of space-time, a conjugate momentum

Ψ µ 0 •••µ 2i (2i+1) := Ψ µ 0 [µ 1 •••µ 2i ] (2i+1) := ∂L ∂ ∂ µ 0 Φ (2i) µ 1 •••µ 2i (2.2.6a) = ∂ µ 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) , ( 2

.2.6b)

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)
(2.2.7a)

= (2i + 1) ∂ [µ 0 Φ µ 1 •••µ 2i ] (2i) + ∂ µ Φ µµ 0 •••µ 2i (2i+2) , ( 2

.2.7b)

we can prove that

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

Let us now prove eq. (2.2.8). By substituting the expression of

Φ µ 0 •••µ 2i (2i+1)
and

Φ µ 2 •••µ 2i (2i-1)
(eq. (2.2.7b)) in eq. (2.2.8), and after a straightforward identification with eq. (2.2.6b), 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, one successively obtains

(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)
.

The Euler -Lagrange equation reads

∂L ∂Φ (2i) µ 1 •••µ 2i = -V [µ 1 •••µ 2i ] (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) , ie Φ µ 1 •••µ 2i (2i) = -V [µ 1 •••µ 2i ] (2i)
.

(2.2.9)

Summarising, 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, 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 equations of motion 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
, and

(2.2.10a)

(2i + 1) ∂ [µ 0 Φ (2i) µ 1 •••µ 2i ] + ∂ µ Φ (2i+2) µµ 0 •••µ 2i = ∂H ∂Φ µ 0 •••µ 2i (2i+1) , (2.2.10b)
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 , (2.2.11)
is the Hamiltonian density associated to L.

In order to identify the form of eq. (2.1.6), we need to flatten all these indices (ie 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 , (2.2.12)
with µ 0 := -1, and where the stacking of sums can be re-expressed in term of generalised harmonic numbers as well.

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

ζ a = ζ a=fl(2i,α 1 •••α 2i ) := Φ (2i) α 1 •••α 2i (2.2.13a) = ζ a=fl(2i+1,α 1 •••α 2i+1 ) := Φ α 1 •••α 2i+1 (2i+1) , ( 2
M µ a b ∂ µ ζ b = ∂ a H , (2.2.13c) provided that M µ a=fl(i,α 1 •••α i ) b=fl(j,β 1 •••β j ) := (2.2.13d) (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; 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.2.3), and for notational convenience, we create aliases for these 2 2 fields as

Φ (0) ≡: φ , (2.2.14a)
Ψ 0

(1) = Φ 0 (1) ≡: ψ 0 , (2.2.14b)

Ψ 1 (1) = Φ 1 (1) ≡: ψ 1 , (2.2.14c) Φ (2) 01 = -Φ (2)
10 ≡: γ , (2.2.14d)

Ψ 001 (3) = -Ψ 010 (3) = Φ 1 (1)
≡: ψ 1 , and (2.2.14e)

Ψ 101 (3) = -Ψ 110 (3) = Φ 0 (1) ≡: ψ 0 . (2.2.14f)
Here, we have introduced one new dynamic field, γ, independent of φ. See fig. I.9 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 (φ) + Ṽ (γ) , (2.2.15
) 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 γ , and
(2.2.16a)

ψ 1 := ∂L ∂(∂ 1 φ) = ∂L ∂(∂ 0 γ) = -∂ 0 γ -∂ 1 φ . (2.2.16b)
From the action,

S := d 2 x L(φ, ∂ 0 φ, ∂ 1 φ, γ, ∂ 0 γ, ∂ 1 γ) , (2.2.17)
one obtains the two equations of motion for φ and γ respectively:

δS δφ = 0 = ∂L ∂φ -∂ 0 ∂L ∂(∂ 0 φ) -∂ 1 ∂L ∂(∂ 1 φ) = ∂L ∂φ -∂ 0 ψ 0 -∂ 1 ψ 1 = -V (φ) -∂ 0 2 φ -∂ 0 ∂ 1 γ + ∂ 1 ∂ 0 γ + ∂ 1 2 φ = -V (φ) -φ , ie φ = -V (φ)
, and (2.2.18)

δS δγ = 0 = ∂L ∂γ -∂ 0 ∂L ∂(∂ 0 γ) -∂ 1 ∂L ∂(∂ 1 γ) = ∂L ∂γ -∂ 0 ψ 1 -∂ 1 ψ 0 = Ṽ (γ) + ∂ 0 2 γ + ∂ 0 ∂ 1 φ -∂ 1 ∂ 0 φ -∂ 1 2 γ = Ṽ (γ) + γ , ie γ = -Ṽ (γ) . (2.2.19)
Hence, the equations of motion for φ and γ are indeed independent since the interaction terms in the Lagrangian (∂ 0 φ ∂ 1 γ and -∂ 0 γ ∂ 1 φ) only produce cross derivatives of the fields, that cancel thanks to Schwarz's theorem.

The Hamiltonian density is

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

ζ := φ ψ 0 ψ 1 γ t , (2.2.21a)
the dynamics are fully described by eq. (2.1.6) provided that

M 0 :=     0 -1 0 0 1 0 0 0 0 0 0 1 0 0 -1 0     , and
(2.2.21b)

M 1 :=     0 0 -1 0 0 0 0 1 1 0 0 0 0 -1 0 0     .
(2.2.21c)

These two matrices are no longer degenerate (their eigenvalues are +i, +i, -i and -i), and they define a multi-symplectic structure on the phase space. The block structure of these matrices highlight once again the independence of the two theories: there are no terms coupling the φ-theory to the γ-theory in the De Donder -Weyl -Hamilton equations.

Recapping, 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 setup, 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.

Returning to the particular setting of the non-linear wave equation in 1 + 1 dimensions, we want to stress that, posing Ṽ (γ) = 0, one obtains the same result as in [START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian {PDEs} that conserve symplecticity[END_REF].

In this paper, the authors consider γ as a Lagrange multiplier associated with the constraint ∂ 0 ψ1 + ∂ 1 ψ 0 = 0 without further explanations. This constraint is in fact the equation of motion for γ, but this does not provide any reason why this construction allows one to solve the problem of degeneracy affecting the multi-symplectic structure. The constraint imposed by γ can be initially seen as imposing equality of the cross derivatives 1 (which is not guaranteed anymore since the derivatives of a field are regarded as independent variables in the DW Hamiltonian formalism). However, introducing γ affects the definition of ψ 0 and ψ 1 , and the constraint it imposes (originally

∂ 0 ∂ 1 φ = ∂ 1 ∂ 0 φ
) is no longer exactly the equality of cross derivatives (it becomes

∂ 0 2 γ + ∂ 0 ∂ 1 φ = ∂ 1 2 γ + ∂ 1 ∂ 0 φ).
This does not matter, as if one tries to generalise this constraint in dimensions higher than 1 + 1, it does not work (and the multi-symplectic structure remains partially degenerated). The issue with this procedure is mainly that it does not add enough fields; it misses many of the necessary fields to completely remove the degeneracy of the multi-symplectic structure. This is only achieved (in any dimension) using the construction we have introduced in this section, and which deeply relies on the notion of canonical conjugate field.

Finally, we want to stress that from the numerical point of view, to integrate a Ddimensional 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 shall first prove the conservation of the multi-symplectic structure under the Hamiltonian flow. Then, we shall 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.1.6) which reads

M µ ab ∂ µ ζ b d a = ∂ a H d a , ie ω µ (∂ µ ζ, •) = dH . (2.2.22)
Now, taking the exterior derivative of it, one successively obtains

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, one has

∂ µ ω µ = 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.2.23)
By definition, a multi-symplectic integrator is a numerical method that exactly preserves the discrete version of eq. (2.2.23).

Stress-energy tensor and charges

The stress-energy tensor

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

T µν := 1 4 [ω µ (∂ ν ζ, ζ) + ω ν (∂ µ ζ, ζ) + 2 η µν ω κ (ζ, ∂ κ ζ)] + η µν H =: T νµ . (2.3.1)
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, . . . ), one has

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

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

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

∂ µ T µν = ∂ µ T νµ = 0 . (2.3.3)
Therefore, the stress-energy tensor is locally conserved (on-shell).

Furthermore, let us stress that definition (2.3.1) of the stress-energy tensor leads to the expected result that it is traceless for a bi-dimensional massless free field theory. Consider the Hamiltonian density of a massless free field theory in any dimension:

H 0 := 1 4 ω µ (∂ µ ζ, ζ) . (2.3.4)
(using the Euler -Lagrange equation and the fact that the global variation commutes with the space-time derivatives). Hence, the variation of the action, which by definition vanishes on-shell,

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

j µ := ∂L ∂(∂ µ φ i ) ∆φ i + L µ (2.3.12a) = ∂L ∂(∂ µ φ i ) δφ i - ∂L ∂(∂ µ φ i ) ∂ ν φ i -η µν L ν . (2.3.12b)
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 ν , (2.3.13)
since all the fields (whatever their nature) have no local variation (the Jacobian is the identity). Now, the stress-energy tensor is naturally defined as the collection of the Noether currents (each being associated with the global translational invariance along a direction of space-time). Its definition is then

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

∂ µ Θ µν = 0 . (2.3.15)
Note that this tensor is not necessarily symmetric but it can be symmetrised following the popular Belinfante procedure (see eg [START_REF] Peter | Primordial cosmology[END_REF]).

Using eqs. (2.1.3) and (2.1.4), the stress-energy tensor is expressed in terms of the Hamiltonian density as

-Θ µν = η µν H -η µν ψ i κ ∂ κ φ i + ψ i µ ∂ ν φ i , ie (2.3.16) = 1 2 [ω µ (∂ ν ζ, ζ) + η µν ω κ (ζ, ∂ κ ζ)] + η µν H , (2.3.17) since ω µ (∂ ν ζ, ζ)/2 produces terms of the form ψ i µ ∂ ν φ i -φ i ∂ ν ψ i µ /2
, which, after integration by parts, give ψ i µ ∂ ν φ i , and ω κ (∂ κ ζ, ζ)/2 produces terms of the form ψ i κ ∂ κ φ i (again, after integration by parts). Then, adding a term to restore its symmetry, the stress-energy tensor is indeed given (up to a sign) by eq. (2.3.1).

Charges

To conclude this section, let us recall the definition of 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µ .
(2.3.18)

The stress-energy tensor is locally conserved (eq. (2.3.3)) 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 , (2.3.19)
where j runs on 1, d .

Summary

We have first presented the De Donder -Weyl covariant Hamiltonian formulation of classical field theory, and we have shown that it leads to a multi-symplectic 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 multisymplecticity. 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.1.6) of the De Donder -Weyl -Hamilton equations.

Now we claim that, whenever a pde can be written in the form (2.1.6), and has a nondegenerate multi-symplectic structure, the multi-symplectic integrator in the light-cone coordinates (msilcc) can be applied; we are going to present how to implement it in the next section.

Multi-symplectic integrators

The method we have developed and we are now going to introduce, the multi-symplectic integrator in the light-cone coordinates (msilcc), is a mixing of ideas from [START_REF] Bridges | Multi-symplectic integrators: numerical schemes for Hamiltonian {PDEs} that conserve symplecticity[END_REF] and [START_REF] Boyanovsky | Approach to thermalization in the classical ϕ 4 theory in 1 + 1 dimensions: Energy cascades and universal scaling[END_REF]. The msilcc method is a centred box 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, inspired from 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], which has the advantage to restore the locality of the method (ie there is the same number of unknowns as equations in each cell).

After finishing writing this thesis we were informed that a similar work was already done in [START_REF] Mclachlan | The multisymplectic diamond scheme[END_REF]. In this paper, the authors introduce the same method but only in 1 + 1 dimensions; the difference is that they do so with a different discretisation step in space and time (regardless of preserving the covariance of the theory). Still only in 1 + 1 dimensions, they also introduce a high-order approximation (of the Runge -Kutta kind) that extends the centred box scheme. At the end of this section, we shall discuss the possible ways of improving our generalisation to higher dimensions by taking advantage of the extra ideas introduced for the 1 + 1 dimensional case in [START_REF] Mclachlan | The multisymplectic diamond scheme[END_REF].

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

Preliminaries

The (partitioned) Euler method

The objective of the method is to generalise the construction presented in section 1.4 to the case of a field theory.

The first step is to split the state vector into several components in such a way that, approximating the DW Hamilton equations by mixing explicit and implicit Euler rules, allows to write the discrete equations of motion at the Midpoint of all derivatives.

The solution is to split the state vector into its n components; all of them are discretised on a D-dimensional hypercubic lattice of spacing δ, but with different origins. The relative origin of these lattices follow the positioning (modulo a δ /2 factor) of the corresponding field (the component of the state vector) in the hypercube constructed when resolving the problem of the degeneracy of the multi-symplectic structure (fig. I.9). ie in each hypercubic cell of the space-time lattice, the fields are positioned as if the no-degeneracy hypercube where multiplied by δ /2 and placed in such a way that φ coincide with a node of the space-time lattice. Then, in the equations of motion, we mix the implicit and explicit Euler approximations such that all the derivatives are expressed by a finite difference that belongs in an hypercube of size δ, centred in φ.

It can be easily proven that this method collapses with the usual approximation of the Laplacian

∆φ ≈ φ n+1 -2φ n + φ n-1 δ 2 .
Later, we shall directly use this approximation without constructing it from the partitioning approach (see section 4.2). It is straightforward to generalise in higher dimensions.

This approximation is an explicit multi-symplectic integrator of first-order in spacetime; the proof can be easily obtained by generalising the results of section 1.4. However, as we shall observe in section 4.2, this approximation is not robust to the nonlinearity; from the von Neumann linear stability analysis, we know that an equal lattice spacing in all the directions of space-time (ie that preserves the covariance) respects the Courant -Friedrichs -Lewy (CFL) condition [START_REF] Courant | On the Partial Difference Equations of Mathematical Physics[END_REF]. It is actually the bound imposed by the CFL condition, but it is still in the stability range. Nevertheless, due to the non-linearity, the approximation rapidly becomes unstable; leading to a divergence of the solution. The stronger is the non-linearity, the faster the integrator diverges. In contrast, taking a different lattice spacing in space and time (ie explicitly breaking the covariance), allows one to move away from the bound of the CFL condition, and the approximation becomes much less sensitive to the non-linearity. Hence, the only way to make this approximation robust to the non-linearities requires to explicitly breaking the covariance of the theory; which is unsatisfactory.

The centred box scheme

The centred box scheme, introduced1 in [42], is a numerical approximation based on a concatenation of Midpoint rules. It is the starting point for the msilcc method. Let us now give a sketch of this approximation.

The idea is to discretise the state vector on an hypercubic lattice and to write down all the approximations at the centre of the elementary cell. It is a natural generalization of the Midpoint rule; it expresses all the quantities at the middle of a D-dimensional hypercube: the value of the field at the centre is obtained by averaging over its values on the 2 D vertices of the cell. The derivative along the µ-direction is the finite difference of two average values, each taken by averaging on the 2 d vertices of the d-hypercube at the corresponding end of the µ-direction.

This approximation is a second-order multi-symplectic integrator. The proof of multisymplecticity can be easily inferred from the demonstration we shall carry out for the msilcc method; as well as the proofs of all the others properties as well.

This integrator is very accurate and possesses a number of desirable properties (in particular, it preserves the covariance of the theory); however, it has a severe drawback: its non-locality. The centred box scheme is not well-defined locally: on each D-hypercube, one has 2 d n unknowns (the n fields at later time) but only n equations (the n equations of motion). Nevertheless, this integrator is well-defined at the global level: each unknown is shared in 2 d adjacent cells; hence, there are a total of n N unknowns (where N is the number of cells in the system) for the same amount of equations. However, solving at the global level requires much more memory and, if the problem is treated on a parallel computer, also requires an enormous amount of communications 1 between processes. Such an algorithm has a deplorable scalability. This method will not be treated in the section 4 since it is inapplicable in practice.

The aim of the msilcc method is to take advantage of the impressive properties of the centred box scheme (by using the same concatenation of Midpoint rules), but with a restored locality (by the use an optimised lattice).

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 , and (3.2.1a)

xj := √ 2 x j + δ j d x0 , (3.2.1b)
the inverse of which is

x 0 := 1 √ 2 d µ=0
xµ , and (3.2.2a)

x j := 1 √ 2 xj -δ j d x0 . ( 3

.2.2b)

Defining ∂µ := ∂ /∂ xµ , the associated vector basis is

∂0 = ∂ 0 -∂ d √ 2
, and (3.2.3a)

∂j = ∂ 0 + ∂ j √ 2 , (3.2.3b)
and its inverse

∂ 0 = ∂0 + ∂d √ 2
, and (3.2.4a)

∂ j = √ 2 ∂j -∂ 0 . ( 3

.2.4b)

1 communications are the bottleneck of supercomputers. One can remark that (3.2.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 in the usual definition, which mixes the time with only one of the dimensions of space), allowing one to treat all the directions of space-time in the same manner inside a cell (see below, eq. (3.2.6) and fig. I.10).

∂ 0 • • • • • D = 0 + 1 • • • • • • • • • ∂ 1 ∂ 0 • • • • D = 1 + 1 • • • • • • • • • • • • • • • ∂ 1 ∂ 2 ∂ 0 • • • • • • D = 2 + 1
Let us now sample the space-time manifold, M, using the lattice

M := n = o + δ ňµ ∂µ ňµ ∈ Z, n ∈ M , (3.2.5)
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 , (3.2.6)
if all the vertices of the cell belong in M (possibly using periodic boundary conditions).

Figure I.10 represents how the lattice looks like in low dimensions. On this lattice, the derivatives along the original directions of space-time (highlighted in colours) will simply be approximated using the Midpoint rule: they only involve two points of M , and this makes the method local.

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 , (3.2.7a) n R +| ∂µ := n + δ R +| ∂µ = n R -| ∂µ + δ ∂µ = n + δ ∂ 0 + σ ∂ ρ √ 2 , (3.2.7b) n L -| ∂µ := n + δ L -| ∂µ = n R +| ∂µ -σ δ √ 2 ∂ ρ = n + δ ∂ 0 -σ ∂ ρ √ 2
, and (3.2.7c)

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

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

The numerical approximation scheme

Definition

The approximation rules are constructed by applying the centred box scheme into each square of the elementary cell.

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
, and (3.3.1a)

∂µ f {ζ a } n C ≈ ϵ f {ζ a } n C , (3.3.1b)
with n C defined in eq. (3.2.8), and where

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

ϵ f {ζ a } n C := 1 δ f 1 2 ζ a n L +| ∂µ + ζ a n R +| ∂µ -f 1 2 ζ a n L -| ∂µ + ζ a n R -| ∂µ . (3.3.2b)
For the moment we do not know whether these approximation rules respect the rules of differential calculus (we shall 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 (ie 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, one obtains

∂ 0 ζ a n C = ∂0 + ∂d √ 2 ζ a n C ≈ [• • •] , and (3.3.3a) ∂ j ζ a n C = 2 ∂j -∂0 -∂d √ 2 ζ a n C ≈ [• • •] . (3.3.3b)
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 = (3.3.10a) - √ 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 , (3.3.10b) γ 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 , (3.3.10c
)

ψ 1 j n+1 -ψ 1 j n-1 -σ n ψ 0 j+σn n + σ n ψ 0 j n = 0 . (3.3.10d)
These are the equations used to integrate the λ φ 4 theory with the msilcc method in section 4. Again, they are not implicit but much too complicated to write down in an explicit form. Hence, we treat them as implicit equations and we solve them using the Levenberg -Marquardt algorithm for non-linear least-squares [START_REF] Madsen | Methods for non-linear least squares problems[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF].

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 , ie ≈ ζ a ∂µ ζ b n C + ζ b ∂µ ζ a n C , (3.4.1)
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.4.2)
First of all, eq. (3.4.1) 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 any approximation).

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

Preservation of Schwarz's theorem

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 limited to

ζ a , ζ a ∂µ , ϵ ζ a . (3.4.3)
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 . (3.4.4)
Following the same reasoning as when the field was ζ a defined on M , one has (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 . ( 3.4.5)
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) , (3.4.6)
proving the Schwarz's theorem 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) . (3.4.7)
1 when applied on objects of the collection z k (defined in section 3.4), 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
, and

z k (n) := 1 2D m∈cell(n-δ C ) z k (m) .

Local approximate conservation of the stress-energy tensor

Let us first obtain two preliminary results. On the one hand, one has the commutativity of the oriented average value (eq. (3.4.2)) 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) . (3.5.1)
On the other hand, one has 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) . (3.5.2)
Now we consider the non-symmetrised part of the stress-energy tensor,

T µν := 1 2 [ω µ (∂ ν ζ, ζ) -η µν ω κ (∂ κ ζ, ζ)] + η µν H(ζ) . (3.5.3)
Since T is a tensor (η is a tensor while ω and ∂ are vectors), one obtains

Ť µν := ∂ xµ ∂x ρ ∂ xν ∂x σ T ρσ (3.5.4a) = 1 2 ωµ ∂ν ζ, ζ -ηµν ωκ ∂κ ζ, ζ + ηµν H(ζ) (3.5.4b) ≈ Ť µν ,
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( ζ ) . ( 3
ϵ Ť µν = 1 2 0 -Ďν ωκ Ďκ ζ, ζ ∂κ + Ďν H( ζ ) = Ďν H( ζ ) - 1 2 ωκ Ďν Ďκ ζ, ζ ∂κ ∂ν + ωκ Ďκ ζ ∂ν , Ďν ζ ∂κ = Ďν H( ζ ) - 1 2 ωκ Ďκ Ďν ζ, ζ ∂ν ∂κ + ωκ Ďκ ζ ∂ν , Ďν ζ ∂κ = Ďν H( ζ ) -dH ζ Ďν ζ = Ďν H i ( ζ ) -dH i ζ Ďν ζ , (3.5.6a)
where

H i := H -H q , (3.5.6b)
is the non-quadratic part of H (H q being the quadratic part of the Hamiltonian density).

Accordingly, the msilcc scheme exactly preserves the local conservation of the stressenergy 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. (3.5.6a). Hence, let us approximate ϵ Ť µν by removing the lowest level average value:

ϵ Ť µν Ďν H i (ζ) -dH i ( ζ ) Ďν ζ . (3.5.7)
Obviously, this operation is strictly forbidden! Nevertheless, eq. (3.5.7) is in practice a very good estimator of eq. (3.5.6a). This can be understood if we remember that eq. (3.5.6a) mainly evaluates how much the chain rule is violated on the lattice for nonquadratic functions. 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.5.6a) (or explicitly eq. (3.5.8b)) and eq. (3.5.7) (explicitly eq. (3.5.9)) is negligible, and it is almost impossible to distinguish the two on the numerical results.

The substantial advantage of the estimator (3.5.7) 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 necessary communications.

The λ φ 4 theory in 1 + 1 dimensions

In the case of the λ φ 4 theory in 1 + 1 dimensions, and assuming that the extra field γ is free (ie used as a control parameter), eq. (3.5.6a) explicitly becomes

± := ϵ Ť µν (3.5.8a) = 1 4 δ φ U + φ ± 2 4 - φ D + φ ∓ 2 4 - 1 2 δ (φ U + φ ± -φ D -φ ∓ ) φ U + φ ± + φ D + φ ∓ 4 3 , (3.5.8b)
where ± selects ν = 0 or 1, and where φ U := 1 4 φ j n + φ j n+1 + φ j-σn n+1 + φ j n+2 , (3.5.8c)

φ ± := 1 4 φ j+ σn±1 2 n + φ j-σn∓1 2 n + φ j-σn∓1 2 n-1 + φ j-σn∓1 2 n+1
, and (3.5.8d)

φ D := 1 4 φ j n + φ j n-1 + φ j-σn n-1 + φ j n-2 , ( 3
.5.8e) while eq. (3.5.7) explicitly becomes

± 1 4 δ     φ j+ σn±1 2 n + φ j n+1 2   4 -   φ j n-1 + φ j+ σn∓1 2 n 2   4   (3.5.9) - 1 2 δ φ j+ σn±1 2 n + φ j n+1 -φ j n-1 -φ j+ σn∓1 2 n φ j n + φ j+σn n + φ j n-1 + φ j n+1 4 3 .

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

Let us briefly introduce the idea of 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]: the discretisation rules are applied on the energy instead of the equation of motion, then the constrains on the conservation of the energy are used as a pseudoequation of motion (see below, section 4.3, for a complete presentation).

The equation of motion and the conservation of the stress-energy tensor 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. One can imagine here to proceed in the same way as in the BDdV method, by applying the discretisation 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 ζ Ďν ζ (ie eq. (3.5.6a)).
Indeed, the dH i ζ Ďν ζ term being the source of the errors in the equation of motion, while Ďν H i ( ζ ) evaluates the errors produced in the derivatives of the stress-energy tensor. So, whether the discretisation 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.

Motivation to use the light-cone coordinates

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 centred 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 which is the Midpoint rule case). 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, one has 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 (ie 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 time-derivative of all the fields of the state vector. Hence the necessity to work with a formulation of the problem that will not lead to a degeneracy of the multi-symplectic structure. This is the main reason for requiring the multi-symplectic structure to be non-degenerate.

One can consider to perform the first step of the integration with another procedure, not requiring to remove the degeneracy of the multi-symplectic structure. However, this will be at the cost of breaking the local conservation of the multi-symplectic structure in the first step; hence, producing errors that might have a non-negligible impact on the rest of the integration process.

Alternative discretisation in dimensions higher than 1 + 1

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.5), is an attempt to generalise to higher dimensions the one introduced in [START_REF] Boyanovsky | Approach to thermalization in the classical ϕ 4 theory in 1 + 1 dimensions: Energy cascades and universal scaling[END_REF] for D = 1 + 1 (see fig. I.15 in section 4.3). However, we have experienced some instabilities of the method in dimension D > 1 + 1.

Actually, in the linear case, these instabilities can be demonstrated; using the von Neumann stability analysis, we have observed that the fastest modes of the linear wave equation in D = 2 + 1 are not stable under the msilcc approximation scheme (it is stable in D = 1 + 1). However, this change in the behaviour 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 (ie there are points in space-time that are not contained in any cell). ii. On the other hand, the oriented average, • • • ∂µ (see eq. One can imagine another generalization of the lattice introduced in [START_REF] Boyanovsky | Approach to thermalization in the classical ϕ 4 theory in 1 + 1 dimensions: Energy cascades and universal scaling[END_REF] 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 ν µ ∂ ν , (3.7.1)
where R is a rotation matrix (ie R ∈ SO(D)), such that the direction (1, 1, • • • , 1) is 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 this section. Now, first of all, let us introduce the model which will support the comparison.

4.1 The λ φ 4 theory in 1 + 1 dimensions

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 the second order non-linear pde (the equation of motion):

φ = ∂ 0 2 φ -∂ 1 2 φ = -V (φ) = -φ 1 + φ 2 , (4.1.1)
where x and t are space and time respectively, ∂ 0 := ∂ /∂ct, ∂ 1 := ∂ /∂x, c is a characteristic speed (eg 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; it has been set to one since one can always rescale the field to achieve this, with no loss of generality.

The other parameter in the potential, the one that accompanies the quadratic term, r, has also been set to one in such a way that the potential has only one absolute minimum at φ = 0. At the end of this section we shall quickly investigate the influence of changing this parameter, and especially when it becomes negative.

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]; they can be useful as a first check of the accuracy of a numerical integrator (see section 4.5).

Boundary and initial conditions

As previously mentioned, a finite-difference method can be decomposed in terms of two ingredients: the lattice and the discretisation 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 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 shall 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) , and (4.1.2a)

∂ 0 φ(x + L, t) = ∂ 0 φ(x, t) . ( 4

.1.2b)

1 which belong in the class of the non-linear wave equation, with the potential

V (φ) := r 2 φ 2 + λ 4 φ 4 .
We use an initial condition that complies with the pbc: φ(x, 0) = A sin 2π x L , and (4.1.3a)

∂ 0 φ(x, 0) = 0 . (4.1.3b)
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)) , ie (4.1.4a) = A 2 π 2 L + L 8 + 3A 2 32 . (4.1.4b)
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 behaviour of the solution turns out to be much more complex. 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 as well as 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.3 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 symmetric stress-energy tensor is Conversely, the numerical equivalents of these local conservation laws (eqs. (4.1.6a) and (4.1.6b)) 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 discretisation rules not always fulfil all the rules of differential calculus (Leibniz, . . . ).

T 00 := 1 2 (∂ 0 φ) 2 + 1 2 (∂ 1 φ) 2 + 1 2 φ 2 + 1 4 φ 4 , ( 4
T 11 := 1 2 (∂ 0 φ) 2 + 1 2 (∂ 1 φ) 2 - 1 2 φ 2 - 1 4 φ 4 . ( 4
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 , (4.1.8)
where ν is either 0 or 1. These are global quantities. Integrating over space, the local conservation of the stress-energy tensor leads to the conservation of the charges,

∂ 0 Q ν = 0 . (4.1.9)
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 shall 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 shall 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 time-scale decreases while the time-step remains fixed). So, we expect the errors to be lower at small A than at large A. Secondly, we shall fix A = 10 and observe how the errors behave as a function of time. We shall explore both short and long-time behaviour.

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 Euler method

The Euler method is probably the simplest finite-difference method one can develop. Its ease of use and its efficiency make it a classic. However, we shall 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 + , (4.2.1)
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 The lattice spacing in space and time are chosen to be equal to preserve covariance. This choice of an equal lattice spacing in space and time respects the CFL condition (computed from the von Neumann stability analysis) for the linear wave equation.

M := δ Z /N Z × δ N = δ Z N × δ N , ( 4
where n ∈ N and j ∈ 0, N .

By this sampling process, at a given time, one has 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 Euler scheme

In order to complete the construction of the finite-difference scheme, the second step is to provide the rules that will indicate how to combine the samples of the field (ie 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 Euler's rule

∂ 0 φ(x = j δ, t = n δ) ≈ D ± 0 ϕ j n := ± ϕ j n±1 -ϕ j n δ
, and (4.2.4a)

∂ 1 φ(x = j δ, t = n δ) ≈ D ± 1 ϕ j n := ± ϕ j±1 n -ϕ j n δ , ( 4

.2.4b)

where the + (respectively -) stands for the forward (respectively backward) approximation 1 . These two definitions (+ or -) are inequivalent. The centred Euler'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 .

1 or explicit (respectively implicit).

2 to show this fact, let us define the centred Euler'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 ,
The equation of motion can be approximated in two ways: using forward then backward rules or vice versa, ie ∂ 2 ≈ D -D + or ∂ 2 ≈ D + D -(either forward -forward or backward -backward leads to an inconsistent approximation of the second derivatives 1 ). 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 . (4.2.5)
This algebraic equation is explicit in ϕ j n+1 . The evolution of a given state can then be efficiently obtained using This concludes the definition of the Euler method. Let us now use these rules to obtain the discrete formulation of the stress-energy tensor.

ϕ j n+1 = ϕ j n-1 + ϕ j-1 n + ϕ j+1 n -δ 2 ϕ j n 1 + ϕ j n 2 . ( 4

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 , (4.2.7a)
T 01 ± = T 10 ± := -D ± 0 ϕ D ± 1 ϕ , and (4.2.7b)

T 11 ± := 1 2 (D ± 0 ϕ) 2 + 1 2 (D ± 1 ϕ) 2 - 1 2 ϕ 2 - 1 4 ϕ 4 , (4.2.7c)
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 ± . (4.2.8b)
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.

thus leads to two independent sub-lattices (the odd one and the even one) and is therefore not a valid approximation.

1 since

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

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

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 ± , (4.2.9)
where Q 0 ± is the charge defined as

Q 0 ± := δ N -1 j=0 T 00 ± . (4.2.10)
One can also envisage to combine these two definitions as We observe that both E+ and E-vary 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 discretisation 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).

E ave := E+ + E- 2 . ( 4
We also notice that there is no apparent change in the amplitude of the deviations as time elapses up to a time-scale at which the energies rapidly diverge. The divergence of 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 , (4.3.10)
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 pseudoequation 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 . ( 4

The stress-energy tensor

The 00 component of the stress-energy tensor was defined in eq. ( 4 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
, and (4.3.13a)

Ď∓ 0 T 01 ± -T 11 ± j n + Ď∓ 1 T 01 ± + T 11 ± j n =: √ 2 1 ± j n , (4.3.13b)
where

ϱ 0 T µν ∓ j n := ± T µν ∓ j-σn 2 n±1 -T µν ∓ j± 1 2 n δ
, and (4.3.14a)

ϱ 1 T µν ∓ j n := ± T µν ∓ j-σn 2 n±1 -T µν ∓ j∓ 1 2 n δ , (4.3.14b)
with both µ and ν being either 0 or 1. In practice, the two definitions of the residues behave in the same way and we shall 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 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 discretisation process; the difference decreases with a better sampling of the initial condition. We also observe that the total energy is exactly conserved, as expected.

E± n := Q 0 ± n , ( 4
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 4.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 properties

The multi-symplectic integrator in light-cone coordinates (msilcc) has been already exhaustively described in section 3. However, let us recall that it is designed such that the discretisation 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. Nevertheless, 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.

In fig. I.17, we show how the total energy (which is here uniquely defined) behaves in time. We first observe that there seems 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 discretisation 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-discretisation 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 in details how the different methods exposed here preserve the local conservation laws as well as the charges.

Numerical results

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 4.1.

Influence of the non-linearity

The first situation will explore the influence of the non-linearities (in coordination with the influence of the quality of the lattice spacing). Figures I. [START_REF] Atanasiu | Distinguished torsion, curvature and deflection tensors in the multi-time Hamilton geometry[END_REF] and I.19 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 (ie 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 4.1).

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

We emphasise that there are some missing data-points for the Euler method (beyond 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 Euler 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.

As already mentioned, taking a different lattice spacing in space and time, allows one to move away from the bound of the CFL condition, and the approximation becomes much less sensitive to the non-linearity. However, this is unsatisfactory since it explicitly breaks the covariance.

To pursue the remarks on the Euler 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.

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 have shown in section 3.5). 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. I. [START_REF] Donder | Théorie invariantive du calcul des variations[END_REF]. 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, here double floating-point precision) than the msilcc method which is itself better than the Euler 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 Euler method.

In conclusion, Euler's method presents a not so bad conservation of the stress-energy tensor as well as of the charges. However, when implemented such that it respects the covariance of the theory, it is unfortunately rapidly destabilised by the non-linearities. Concerning the BDdV method, it presents very poor local conservation properties that, quite surprisingly, lead to excellent conservation of the charges (due to a deceptive cancelation of the errors). Finally, the msilcc method behaves more like Euler's but with much better conservation features.

long-time behaviour

We now explore the long-time properties. Figures I.20 and I.21 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. I.18 and I.19 still hold and we shall only describe the time behaviour here.

Firstly, we observe that the Euler 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 behaviour: after a rapid evolution (of the order of the characteristic time-scale of the system, as we ), the error stabilises 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.

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 (ie 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 4.1 and for A = 10). In some situations the differences can become dramatically larger (up to 20%) as we show on fig. I.22 that represents the field after an integration time t /L = 1, still with periodic boundary conditions, but with an initial state that corresponds to 4.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 data-points represent the field evolved by the different numerical methods (diamonds for Euler, triangles for BDdV and circles for msilcc). On this scale we see no difference between the exact and the numerical solutions obtained with the Euler and msilcc methods.

∆ ∂ 0 Q 1 = 0 msilcc BDdV Euler max ≤t,x ∆(␣) max x ∆(␣)
a particular solution of eq. (4.1.1) in term of a Jacobi elliptic function [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF]:

φ(x, t) = 2 k 2 1 + k 2 sn p (x + ν t) k 2 , where (4.5.1a)
ν := 1 + 1 (1 + k 2 ) p 2 , and (4.5.1b)
p := 4 K k 2 L . ( 4.5.1c)

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 now want to show that the msilcc method behaves just as well for a potential in double well. Figure I.23 shows the error committed on the conservations of the charges and of 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 (r positive or negative). 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 already proven in section 3.5: the deviations from the conservation of the stress-energy tensor only arise with the nonlinear part of the Hamiltonian. In fact, this feature disappears in the large |r| limit for numerical reasons: multiplying only certain terms in the equations of motion by a large value has the tendency to increase the effects of the truncation errors.

-0.1 -1 -10 -100 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 settled to 1.

10 -16 10 -15 max t,x ∆ ∂ µ T µ0 = 0 max t,x ∆ ∂ µ T µ1 = 0 max t,x ∆ ∂ 0 Q 0 = 0 max t,x ∆ ∂ 0 Q 1 = 0 10 -5 10 -4 ∆(∂ µ T µν = 0), ∆(∂ 0 Q ν = 0) -0.1 0

Conclusion

The purpose of this work was to introduce a new numerical method to integrate partial differential equations stemming from the Hamiltonian dynamics of field theories. The method is a centred 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; especially in the problem we were initially interested in.

As mentioned earlier, we have recently learnt that a similar method has already been introduced by a different group [START_REF] Mclachlan | The multisymplectic diamond scheme[END_REF] with similar results. While they focus on the 1 + 1 dimensional case without any assumptions on the lattice spacing in space and time, our approach is more focused on a generalisation in any dimension. We also have imposed an equal lattice spacing in space and time to respect the covariance since it is an important symmetry of the theory as well as of the underlying multi-symplectic structure. Combining these two approaches seems to be promising to address the instability issues we encountered in dimensions higher than 1 + 1.

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.

Chapter II

Critical percolation in ferromagnetic Ising spin models 1 Introduction and preliminaries

Introduction

Statistical systems can be out-of-equilibrium in many different ways; as already stated, the non-equilibrium dynamics occurring after an instantaneous quench is particularly interesting, especially when the quench crosses a second order transition where the ordered phase spontaneously breaks an internal symmetry of the system. A natural question to address is how order emerges from the initially disordered state; it the so called coarsening dynamics, that progressively makes the short length scales to acquire the properties of the equilibrium target state.

In recent years, the interplay between percolation and coarsening [START_REF] Bray | Theory of phase-ordering kinetics[END_REF]121,[START_REF] Henkel | Non-equilibrium phase transitions[END_REF] bidimensional spin models was studied in quite some detail. A series of papers proved that the critical and sub-critical instantaneous quenches of the bi-dimensional ferromagnetic Ising model rather quickly approach a critical percolation state1 and later undergo the coarsening phenomenon. More precisely, in the quenches performed, the evolution starts from a totally random initial configuration mimicking equilibrium at infinite temperature and later evolve with different microscopic stochastic spin updates. This feature was demonstrated with extensive numerical simulations of the Glauber -Ising model for ferromagnetism [110,[START_REF] Blanchard | How soon after a zero-temperature quench is the fate of the Ising model sealed?[END_REF][START_REF] Blanchard | Critical percolation in the dynamics of the bidimensional ferromagnetic Ising model[END_REF][START_REF] Sicilia | Domain growth morphology in curvature-driven two-dimensional coarsening[END_REF] and the Kawasaki model for phase separation [START_REF] Sicilia | Geometry of phase separation[END_REF][START_REF]Phase separation and critical percolation in bidimensional spin-exchange models[END_REF], quenched into their symmetry broken phases. The effects of weak disorder were considered in [START_REF] Corberi | Coarsening and percolation in a disordered ferromagnet[END_REF]136]; the voter model dynamics was investigated in [START_REF] Dornic | Large deviations and nontrivial exponents in coarsening systems[END_REF][START_REF] Tartaglia | Percolation and coarsening in the bidimensional voter model[END_REF]; and, especially relevant for the present study, quenches to the critical point of the bi-dimensional ferromagnetic Ising model were considered in [START_REF] Blanchard | A morphological study of cluster dynamics between critical points[END_REF][START_REF]Response of non-equilibrium systems at criticality: ferromagnetic models in dimension two and above[END_REF]. The early approach to critical percolation also explained why zero temperature quenches of the bi-dimensional Ising model often get blocked in metastable states with infinitely long-lived flat interfaces [START_REF] Barros | Freezing into stripe states in two-dimensional ferromagnets and crossing probabilities in critical percolation[END_REF][START_REF]Zero-temperature relaxation of three-dimensional Ising ferromagnets[END_REF][START_REF] Olejarz | Zero-temperature freezing in the three-dimensional kinetic Ising model[END_REF][START_REF]Fate of 2D kinetic ferromagnets and critical percolation crossing probabilities[END_REF][START_REF]Zero-temperature coarsening in the 2d Potts model[END_REF][START_REF] Spirin | Fate of zero-temperature Ising ferromagnets[END_REF][START_REF]Freezing in Ising ferromagnets[END_REF]. Metastable states in quenches from the critical point to zero temperature were considered in [START_REF] Blanchard | Frozen into stripes: Fate of the critical Ising model after a quench[END_REF].

In statistical physics studies, quenches are assumed to be instantaneous. Indeed, the relevant time-scales in experimental realisations are such that the cooling time is much shorter than all other time-scales. Instead, in field theoretical models of cosmology, there was interest in determining the cooling rate dependencies induced by a very slow quench across a second order phase transition. The original Kibble arguments for the existence of spatial regions that are not causally connected long after going through the phase transition [START_REF] Kibble | Topology of cosmic domains and strings[END_REF] were complemented by a scaling proposal, by Zurek [START_REF] Zurek | Cosmological experiments in superfluid helium?[END_REF]158]. This argument allows one to estimate the correlation length reached when the system falls out-of-equilibrium when approaching a critical point from the symmetric phase with a weak finite speed. The interest in counting the number of topological defects left over after crossing the phase transition triggered by cosmology [START_REF] Kibble | Phase-transition dynamics in the lab and the universe[END_REF], prompted condensed-matter experimental physicists to try these measurements in the lab. This kind of experiments were first performed in Helium-3 [START_REF] Bäuerle | Laboratory simulation of cosmic string formation in the early Universe using superfluid 3 He[END_REF] and liquid crystals [START_REF] Chuang | Cosmology in the laboratory: Defect dynamics in liquid crystals[END_REF] more than twenty years ago. The subject was recently revived by the realisation of cold atom experiments in which the samples are taken across the critical region with a finite speed [START_REF] Chomaz | Emergence of coherence via transverse condensation in a uniform quasi-two-dimensional Bose gas[END_REF][START_REF] Donadello | Creation and counting of defects in a temperature-quenched Bose-Einstein condensate[END_REF][START_REF] Lamporesi | Spontaneous creation of Kibble-Zurek solitons in a Bose-Einstein condensate[END_REF][START_REF] Navon | Critical dynamics of spontaneous symmetry breaking in a homogeneous Bose gas[END_REF]. New studies in ion crystals [START_REF] Del Campo | Structural defects in ion chains by quenching the external potential: the inhomogeneous Kibble-Zurek mechanism[END_REF][START_REF] Silvi | Crossover from Classical to Quantum Kibble-Zurek Scaling[END_REF][START_REF] Ulm | Observation of the Kibble-Zurek scaling law for defect formation in ion crystals[END_REF], bi-dimensional colloidal suspensions [START_REF] Deutschländer | Kibble-Zurek mechanism in colloidal monolayers[END_REF][START_REF] Reichhardt | Multi-step ordering in kagome and square artificial spin ice[END_REF] have also been recently performed. Two recent reviews give a more complete summary of the status of this field [142,[START_REF] Del Campo | Universality of phase transition dynamics: Topological defects from symmetry breaking[END_REF].

Studies of cooling rate dependencies in statistical physics models were performed in a number of papers: for instance, the bi-dimensional Ising model with non-conserved order parameter dynamics was considered in [START_REF] Biroli | Kibble-Zurek mechanism and infinitely slow annealing through critical points[END_REF][START_REF]Nonequilibrium critical dynamics of ferromagnetic spin systems[END_REF], and the xy model in two dimensions (planar spins) in [START_REF] Jelić | Quench dynamics of the 2d XY model[END_REF] (the latter is relevant to discuss the recent experimental activity in Bose -Einstein condensates and colloidal suspensions). In the former model, the phase transition is a conventional second order one: from a symmetric to a symmetry broken phase; in the latter case, the transition is of Berezinskii -Kosterlitz -Thouless (BKT) kind and the target is a critical phase. The aim of these papers was to show that, contrary to what was usually claimed in the Kibble -Zurek literature, the dynamics are not frozen after the system falls out-of-equilibrium close to the critical point (be it second order or BKT). The critical or subcritical dynamics, at continuously changing control parameters, let the dynamic correlation length go on growing in time. Scaling arguments were used in these papers to derive the dependence of the growing correlation length, and hence the number of topological defects, as a function of time and cooling rate; they were favourably compared to the outcome of numerical simulations. Exact results for the one dimensional Ising chain and a variety of cooling procedures were derived in [START_REF] Krapivsky | Slow cooling of an Ising ferromagnet[END_REF]. The spherical ferromagnetic model with exponentially fast cooling was treated, also analytically, in [START_REF] Picone | Competition between dynamic and thermal relaxation in non-equilibrium critical spin systems above the critical point[END_REF]. A one-dimensional non-equilibrium lattice gas model with a phase transition was treated in [START_REF] Jain | Critical dynamics of classical systems under slow quench[END_REF]. Extensive numerical simulations of models for two dimensional atomic gases were very recently presented in [START_REF] Comaron | Dynamical critical exponents in drivendissipative quantum systems[END_REF][START_REF] Davis | Formation of Bose-Einstein condensates[END_REF][START_REF] Larcher | Critical Dynamics in Quenched 2D Atomic Gases[END_REF][START_REF] Proukakis | [END_REF]. The evolution of the order parameter in the finite dimensional Ising model slowly cooled to the critical point were studied with different microscopic stochastic rules in [START_REF] Liu | Dynamic scaling at classical phase transitions approached through nonequilibrium quenching[END_REF].

The aim of the work presented in this chapter is to revisit the slow cooling of the bi-dimensional Ising model [START_REF] Biroli | Kibble-Zurek mechanism and infinitely slow annealing through critical points[END_REF][START_REF] Liu | Dynamic scaling at classical phase transitions approached through nonequilibrium quenching[END_REF] paying now special attention to the geometric properties of the domain structures formed when approaching the critical point with a finite speed. The outline of the chapter is the following.

The remaining of this preliminary part, sections 1.2 and 1.3, will be devoted, on the one hand, to recalling some results of site percolation especially relevant for the present study, and, on the other hand, to briefly introduce the concept of stochastic Loewner evolution so that the tools used to characterise the geometry of the domains will make sense.

The second part of this chapter, section 2, will be devoted to introducing the model as well as the observables involved in this study. The model we focus on is the emblematic kinetic ferromagnetic Ising model in two dimensions and on the square lattice. Each site is occupied with probability p, and unoccupied with probability 1-p.

In the represented configuration, and assuming periodic boundary conditions, the largest cluster percolates only if the missing site is occupied (probability p).

In section 3, we establish a reference behaviour for the observables by studying their properties in equilibrium. Next, in section 4, we recall the essential steps to describe the dynamics following an instantaneous quench, both to the critical point and to zero temperature.

Finally, in section 5, we study the effects of a finite cooling rate. We first recall the Kibble -Zurek mechanism, and how to extend it so that the growth of the dynamic correlation length is correctly described in the out-of-equilibrium regime. Next, the fractal geometric properties of the domains walls will be characterised, both when the system reaches the critical point and in the course of the prior cooling process.

Some reminders of site percolation

Site percolation [START_REF] Aharony | Introduction to percolation theory[END_REF][START_REF] Christensen | Complexity and criticality[END_REF][START_REF]Fields, particles and universality in two dimensions[END_REF][START_REF] Saberi | Recent advances in percolation theory and its applications[END_REF] is a purely geometric problem in which particles are placed at the sites of a lattice with probability p ∈ [0, 1] (see fig. II.1). In particular, the properties of the site percolation problem on the bi-dimensional square lattice (Z 2 ) are especially relevant for the present study. The results presented in this section are just a selection of reminders from the literature.

Let us first recall the definition of the concept of cluster: it is a maximal subset of the lattice, constituted of occupied sites, pairwise connected by a path. A path being any sequence of displacements on the lattice from an occupied site to one of its nearest neighbours (occupied as well).

A representative question of percolation theory is to determine the probability, ϑ(p), of having at least one infinite cluster1 . Of course, ϑ(0) = 0 while ϑ(1) = 1, but it can also be proven that ϑ is an increasing function of p. In fact, there exists a percolation threshold, p c , below which there are no infinite clusters (with probability 1 ie ϑ(p < p c ) = 0). Above the threshold, the probability of having at least one infinite cluster becomes strictly positive (ϑ(p > p c ) > 0). Hence, for an infinite system, the site percolation threshold can be defined as

p c := sup p {ϑ(p) = 0} , (1.2.1)
and depends on the geometry and dimension of the lattice; here, for the Z 2 lattice, p c ≈ 0.593. This model undergoes a phase transition:

i. p < p c is the subcritical regime: the system is constituted of small clusters uniformly distributed. ii. p > p c is the supercritical phase: the majority of the system is populated, and it only remains small islands of unoccupied sites. iii. p = p c is the critical point: the behaviour is similar to the one at a thermodynamic second order critical point with universal critical exponents characterising various geometric quantities that one can define. It belongs to the same universality class as a stochastic Loewner evolution with κ = 6 (sle 6 , see section 1.3).

Finally, in the vicinity of the critical point, the typical size of the clusters 1 , ξ, diverges as

ξ ∼ |p -p c | -νp , (1.2.2)
where ν p = 4 /3 for the Z 2 site percolation problem.

To conclude this brief introduction, let us recall some finite size effects. Later, we shall only consider a L × L square lattice (with periodic boundary conditions). However, the notion of infinite cluster has no meaning on such a lattice, and the concept of percolation requires some adjustments. So, on a finite lattice, a cluster is now said to be percolating as long as it crosses the system (from the left to the right, or/and from the top to the bottom).

In this situation, ϑ(p < p c ) does not vanish anymore; however, in practice, it remains negligible as long as L is sufficiently large. Thus, the site percolation threshold is now in the region (hopefully, sufficiently thin) where ϑ(p) starts to be significant.

An interesting side effect of the finite size of the system is that it reveals a noteworthy property at the critical point: the largest percolating cluster is always much bigger than the second largest cluster, of almost one order of magnitude. This is a characteristic feature of percolation. Finally, Ising models (see section 2.1) can be thought of a site percolation problem after performing a one-to-one mapping between Ising spins (= ±1) and the fact that a site is occupied or not. For example, an infinite temperature configuration in which the spins take values ±1 with probability 1 /2 is a random percolation configuration with p = 1 /2, therefore, below the threshold on Z 2 for percolation of a cluster of occupied sites.

Criticality in Schramm -Loewner evolution

As stated, the present study mainly concerns the geometric properties of the interface between domains in kinetic Ising models (that is to say, the interface between clusters of connected spins pointing in the same direction); the out-of-equilibrium dynamics later considered are in fact, as we shall observe, influenced by two critical points. On the one hand, the Ising critical point, and, on the other hand, the one corresponding to the Z 2 site percolation problem at p = p c , and in which domains are still fractal objects but with a different criticality.

Schramm -Loewner evolutions [START_REF] Bauer | 2D growth processes: SLE and Loewner chains[END_REF][START_REF] Cardy | SLE for theoretical physicists[END_REF][START_REF] Duncan | An Introduction to the Loewner Equation and SLE[END_REF][101][START_REF] Henkel | Conformal invariance: an introduction to loops, interfaces and stochastic Loewner evolution[END_REF][START_REF] Kager | A guide to stochastic Löwner evolution and its applications[END_REF][START_REF] Kennedy | A fast algorithm for simulating the chordal Schramm-Loewner evolution[END_REF][105][START_REF] Lawler | An introduction to the stochastic Loewner evolution, Random walks and geometry[END_REF][107][START_REF] Saleur | Exact determination of the percolation hull exponent in two dimensions[END_REF][START_REF] Werner | Part II: Random planar curves and Schramm-Loewner evolutions[END_REF], or stochastic Loewner evolution, denoted sle κ , actually provide a common framework to deal with the interfaces between domains formed at both of these two critical points.

Roughly speaking, a sle κ is a bi-dimensional random path obtained by conformal transformation of a Brownian motion; it is a particular class of bi-dimensional conformal field theory. They are fractal objects whose criticality is characterised by a universal parameter, κ: the stochastic Loewner parameter.

Most of the statistical problems defined on a bi-dimensional lattice map to a sle κ at their critical point. Especially, ferromagnetic Ising spin models at the critical temperature map on sle κ=3 , and the interface of a cluster at the site percolation threshold (p = p c ) belongs to the universality class of sle κ=6 . If κ = 0, there are no criticality: the path is a straight line.

Two typical realisations of discrete sle κ are represented on fig. II.3 for κ = 3 (Ising) and κ = 6 (site percolation), highlighting the variety of possible behaviours depending on the value of the κ parameter.

The tools of Schramm -Loewner processes allow one to extract relevant geometric information. Among the multitude of interesting properties of sle κ , two are of a particular interest here for characterising the interfaces between domains. The first one is the variance of the winding angle, but we shall come back to this point later, in section 2.2. The second one is the Hausdorff dimension: for a closed sle κ , if ξ is the typical size of the enclosed area, the typical length of the interface scales as One can remark that some parts of these trajectories appear quite smooth, and therefore seem not to be fractal anymore. It is actually an issue of the method we have employed to generate these examples: the sampling is not uniform in the curvilinear abscissa. Hence, some samples are too distant, and it inevitably misses the details in between. In fact the fractality should be everywhere the same.

where

D κ := min 2, 1 + κ 8 (1.3.1b)
is the Hausdorff dimension.

In the next section, we introduce the kinetic Ising model; mentioned but not defined so far. We shall also introduce the relevant observables that allow one to characterise the geometry of the interfaces present in this model; in particular, the variance of the winding angle whose theoretical behaviour is obtained using the framework sketched here.

2 Ising models and observables

Ising models

The bi-dimensional ferromagnetic Ising model

In this work we focus on the emblematic ferromagnetic Ising model whose Hamiltonian is

H({σ i }) := -J i,j σ i σ j , (2.1.1)
with J > 0, and where spin variables take only two values: σ i := ±1; the symbol i, j means that the sum is over nearest neighbours only. In particular, we study its bi-dimensional realisation (d := 2) on the square lattice (Z 2 ).

The Hamiltonian tends toward aligning each spin with its nearest neighbours; disorder emerges from the contact with a bath at temperature T . The temperature of the bath will be our control parameter, and by equilibrium, we mean that the system is thermalised with respect to the bath.

The canonical equilibrium properties as a function of the parameter K := β J, with β the inverse temperature, are described by the partition function

Z(K) := {σ i =±1} exp K i,j σ i σ j . (2.1.2)
Hereafter we work with units such that J := 1 and β := 1 /T (k B := 1), ie K := 1 /T . This model undergoes a second order phase transition at

K c := 1 T c := log 1 + √ 2 2 ≈ 0.441 . (2.1.3)
T c is the so called Curie temperature; we shall also refer to it as the Ising critical point. At this critical point, and in equilibrium, the system is conformally invariant and the interfaces of the domains belong to the universality class of sle 3 . Moreover, the Hamiltonian presents a Z 2 symmetry spontaneously broken in the low temperature phase (K > K c ).

The total number of spins in the system is L × L, with L the linear length of the lattice measured in units of the lattice spacing a, and we assume periodic boundary conditions. All the numerical results presented in this chapter were obtained using L = 1024; the snapshots are taken in a system of size L = 128.

Glauber dynamics -kinetic Ising model

The model is endowed with microscopic Monte Carlo stochastic dynamics for the individual spins. The microscopic update rule is the one-spin flip Glauber dynamics which simulates the contact with a heat bath at temperature T .

With P(σ, t) the probability the system is in state σ := σ 0 • • • σ L 2 -1 t at time t, the process is described by the master equation

τ dP(σ, t) dt = σ W σ → σ P σ , t -W σ → σ P(σ, t) , (2.1.4)
where the transition rates for Glauber dynamics are

W σ → σ := L 2 -1 i=0 min 1, e -β δH(σ i →-σ i ) δ -σ i σ i L 2 -1 j=0 j =i δ σ j σ j , (2.1.5)
and τ is a characteristic time-scale.

In practice, the process is as follows: we randomly chose a spin i ∈ 0, L 2 in the system; the spin is flipped (σ i = -σ i ) with a probability p := min 1, e -β δH(σ i →-σ i ) , (2.1.6) where δH(σ i → -σ i ) is the difference produced in the Hamiltonian (2.1.1) by a potential flipping of the selected spin. On Z 2 , β δH can only takes five different values:

-8K, -4K, 0, 4K, or 8K. The process is controlled by the parameter K given by the external inverse temperature of the bath, β, times the exchange parameter, J. We have illustrated it on fig. II.4.

Repeating this process L 2 times constitutes one unit of time in the kinetic Ising model. Hereafter, the time appearing in dynamical studies is always in this unit. On the first line, the selected spin is flipped with probability 1 since that reduces the energy of the system. On the second line, the selected spin is now flipped with probability p = e -8K since that would produce four unfavourable links (-4K) instead of the four favourable ones already present (-4K). The smaller K (ie the higher the temperature of the bath), the more probably the selected spin will be flipped.

Observables

A small number of observables are enough to characterise the geometry of the domains present in the system both in equilibrium and through the out-of-equilibrium dynamics: the correlation length, the variance of the winding angle, and finally, the average occupancy rates of the largest clusters. We define these observables in this section.

Space time correlation function and correlation length

In equilibrium, the correlation of the spin fluctuations

C c (r = |i -j|) := σ i σ j -σ i σ j , (2.2.1)
where r ∈ 0, L /2 , allows one to extract the equilibrium correlation length, ξ eq , with different studies of its decaying properties over distance. For instance, one can extract it from the weighted integral

ξ eq :∼ Λ 0 r ζ C c (r) dr Λ 0 r ζ-1 C c (r) dr , (2.2.2)
with a convenient choice of the power ζ and the cut-off length Λ. In particular, we shall use ζ = 2 and Λ will be chosen as the largest possible distance such that C c (r) remains larger than its statistical fluctuations.

In dynamical studies, the time-dependent space correlation is defined just as in eq. (2.2.1), where the spins are now time-dependent variables; the average is taken over different histories of the dynamics (ie different realisations of the random noises) instead of the canonical statistical ensemble. The procedure in the right-hand-side of eq. (2.2.2) can then be applied to extract the dynamic growing length, ξ(t), that characterises the growth of equilibrium structures at T c or below T c .

Variance of the winding angle

The second observable of interest for characterising the interface between the domains is the variance of the winding angle (wav).

The winding angle, θ( ), is measured on any bi-dimensional curve as a function of the curvilinear abscissa, , as follows. We first chose an origin point for the curvilinear lengths. Then, we measure a reference angle, θ 0 , between a chosen fixed direction and the tangent to the curve at the origin of the curvilinear lengths. Now, for each point on the curve, we define η( ) ∈ [-π, π], the local angle between the same chosen fixed direction as earlier and the tangent to the curve at . Finally, the winding angle is obtained by integrating the variation of the local angle along the curve: For closed curves, after one turn (ie returning to the origin), we have ∆θ = 2nπ, where n ∈ Z is the number of loops. In particular, since the curve is an interface, it cannot cross itself and ∆θ = 0 or ∆θ = ±2π (where the sign changes whether the curve rotates clockwise or anticlockwise). The former (∆θ = 0) means that the interface spans the system from one border to another one, while the latter (∆θ = ±2π) means that the interface encloses a finite area.

θ( ) := θ 0 + 0 dη (2.
Note that on a square lattice η( ) (as well as θ 0 ) only takes four values.

The moments of the winding angle can then be computed by taking its desired power and performing the equilibrium or dynamic statistical averages.

Finally, for the universality class of sle κ (see section 1.3), the average of θ( ) vanishes at all length scales and its variance satisfies [START_REF] Saleur | Exact determination of the percolation hull exponent in two dimensions[END_REF]

θ 2 ( ) = C + 4 κ 8 + κ log , (2.2.4)
where is the curvilinear distance along the curve measured in units of the lattice spacing a, C is an irrelevant constant, and κ takes a universal value depending on the kind of criticality.

Equation (2.2.4) is of the utmost importance. This noteworthy form for the growth of the wav provides an accurate and efficient tool, widely used in the present study, to characterise the geometry (ie the kind of criticality) of the interfaces of the domains in the dynamic model.

Average occupancy rate

The last observable we consider in this work is the average occupancy rate; especially the one of the largest clusters.

As above, in dynamical studies, the average is taken over different realisations of the dynamics (ie of the random noises) while, in equilibrium, the average is taken over the canonical statistical ensemble. Then, in each statistical sample, we sort the clusters with respect to their surfaces; from the biggest to the smallest. Next, the clusters are labeled by their rank in the list established above. Finally, the average occupancy rate of the n th largest cluster is defined as the average of the occupancy rates of the clusters having label n. The occupancy rate obviously being the surface of the domain divided by the total surface of the system (L 2 ).

Note that this process is independently performed at every time in dynamical studies. For instance, in a given realisation, the largest cluster at a given time may be replaced by a completely different domain at the following step.

The number of clusters is not conserved, and if a label is not present in the sorted list (because the latter is too short) the surface of the corresponding cluster is considered to be zero.

This observable allows one to estimate how the system is distributed over all the domains (mainly between the three largest ones).

Equilibrium behaviour

In this section we review some properties of the equilibrium behaviour of the bidimensional Ising model at high temperature and at the critical point. Our objective is to establish a reference equilibrium behaviour for the observables relevant to our study.

Away from the critical point, correlations span a finite distance. The equilibrium correlation length diverges at the Curie critical temperature; in the close vicinity of the critical point, it does so as a power law,

ξ eq (τ ) ∼ τ -ν , where τ := T -T c T c > 0 (3.0.1)
is the distance to the critical point1 , and ν := 1 is the universal critical exponent of the Ising universality class associated to the correlation length. Equation (3.0.1) is only valid in a close vicinity of the critical temperature (τ 1); far from it, there are extra corrections to add, but we do not need them here. Another limitation of eq. (3.0.1) is that it is only valid for an infinite system; if the system size (L) is finite, it limits the growth of the correlation length to a saturation threshold that scales here with the system size as ξ eq (τ = 0) = ξeq ∼ L.

Let us now discuss the equilibrium behaviour of the variance of the winding angle (wav), ie the nature of the interfaces between domains; see fig. II.6. We observe that the wav increases logarithmically on short curvilinear length scales; the value of κ extracted from the slope of θ 2 (log ) is close to 6 at high temperature and close to 3 at T c . This means that, on short length scales, the interfaces of the domains are subject to a conformal invariance (with the criticality of percolation at high temperature and the one of Ising at T c ). There is nothing surprising here. Firstly, at the Ising critical point, the domains obviously have the criticality of the corresponding universality class. Secondly, at high temperature, the Ising model is a percolation problem (correlations are so short that one could argue that the spins are randomly chosen to point up or down with half probability, p = 1 /2). A typical configuration is, therefore, one of a site percolation problem away from its critical point [START_REF] Delfino | Field theory of Ising percolating clusters[END_REF] (recall that, on a square lattice, the critical percolation threshold is at p c ≈ 0.593 > 1 /2 = p). In consequence, on average, there are no percolating clusters in these configurations. This means that the conformal invariance disappears at sufficiently long length scales: ∼ |pp c | -νpDp ∼ 10 2 , where ν p = 4 /3 is the percolation correlation length critical exponent and D p = 7 /4 is the fractal dimension of the interface of a percolation cluster. This leads us to our second remark: at high temperature and long length scales, the wav does not grow logarithmically anymore; it increases much faster. This is, in fact, due to the finite size of the domains. Indeed, since we are far from the critical percolation threshold, the domains remain small, and the overall curvature necessary to close their interface is responsible for a faster growth of the wav. When the temperature decreases the domains swell (like the correlation length), and the wav stops its logarithmic growth at a longer and longer length scale. Obviously, when reaching T c , there is a true conformal invariance, and the wav increases logarithmically on all length scales. Considering only the short length scales, as the temperature decreases, the criticality smoothly evolves from the percolation universality class to the Ising one. This is most clearly shown in fig. II.7 (upper panel) where κ is plotted as a function of T . The slope is extracted from the wav by linear interpolation on short length scales; the longer length scales, where criticality disappears, are excluded from the interpolation set. The Ising criticality is only reached in a close vicinity of the critical point (< 1.1 T c ).

The fact that we observe critical percolation properties in the disordered phase is related to the presence of a critical curve in the temperature-field phase diagram of the bi-dimensional Ising model. It separates a phase with an infinite cluster of parallel spins (at sufficiently large external field) from one without (weak field). This critical curve joins the Ising critical point (Curie temperature and zero field) with the infinite temperature limit at non-vanishing value of the external field, while remaining close to the zero field axis [START_REF]Fields, particles and universality in two dimensions[END_REF]. The vicinity of this line at our working temperatures justifies the fact that we see (finite size) critical percolation geometric properties on the spins clusters.

The last quantity we want to discuss is the average occupancy rate of the largest clusters shown in fig. II.7 (lower panel). Firstly, at high temperature, all the clusters are more or less of the same size. Then, when temperature decreases, the bigger clusters start to grow by absorbing the smaller ones, up to a point (T ≈ 1.1 T c ) where only the two biggest prevail over all the others. Having two coexisting big clusters is a feature of percolation on Ising clusters1 . These two clusters will coexist up to a very close vicinity of the Curie temperature ( 1.01 T c ). In contrast, at the Ising critical point there is only one large cluster (much larger than all the others).

To summarise our observations, at T c , or in its very close vicinity, the system is occupied by only one large geometric cluster having the Ising criticality (κ ≈ 3) at all length scales. See the snapshot at T = T c in fig. II.8. At high temperature, the domains are much smaller. However, on short length scales, they have the geometric properties of critical percolation (κ ≈ 5.5, which is only 5% different from the slope expected with κ = 6). Finally, in between, the criticality smoothly changes from the

Instantaneous quenches

In this section we recall some features of the dynamics after instantaneous quenches to zero temperature and to the critical point, as interpreted from the geometric point of view that we adopt in this work.

Quench to T = 0

The second situation of interest, after the equilibrium, is the one of an instantaneous quench to zero temperature. We consider the following procedure: starting from an equilibrium state at T = 2 T c , at t = 0 we suddenly change the temperature of the bath to zero, ie

T (t) := 2 T c t ≤ 0 0 t > 0 , (4.1.1)
and we observe the further evolution of the system.

In such a procedure, the growing length is known to increase as a power law,

ξ(t) ∼ t 1 /z d , (4.1.2)
where z d := 2 is the dynamical exponent [START_REF] Allen | A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening[END_REF][START_REF] Bray | Renormalization-group approach to domain-growth scaling[END_REF]. Of course, this result holds only for t such that a ξ(t) L: ξ(t) cannot be smaller than the lattice spacing, and it is bounded by the finite system size. The growing length extracted from the correlation function after such a quench is compared to the theoretical expectation in fig. II.9 (upper panel). L is fulfilled. The lower panel still represents the wav, but with a different scaling: the horizontal axis is rescaled following eq. (4.1.4), and d (t) is evaluated through its theoretical expression (∼ t 1 /z d ).

Let us now discuss how the behaviour of the wav evolves in time, as displayed in fig. II.10 (upper panel). At the initial time, the system is in equilibrium, and the wav behaves as described in section 3 (see the yellow curve in fig. II.6). Then, the zero temperature dynamics start to smooth the interfaces: first, on short length scales, then, on longer length scales. This is the first part of the curve and the wav does not increase since a smooth interface has no criticality (κ = 0). In the meantime, the clusters swell, and since the system has not yet realised, at long length scales, that it is at zero temperature (and should have smooth interfaces), it develops the criticality of percolation. This is the second part of the curve; the wav restarts to grow logarithmically. The typical (curvilinear abscissa) length scale that separates these two behaviours is denoted d (t), and is related to the typical size of the domains: Starting from a high-temperature equilibrium state, all clusters are almost of the same size. Next, in the early dynamics, they all grow in the same way. As soon as the correlation length starts to grow, the larger clusters progressively swallow the smaller ones. Indeed, the smaller the clusters, the faster they disappear. This is the so-called coarsening dynamics. In particular, the second largest cluster lengthly coexists1 with the largest one. As already mentioned, this long coexistence of two large clusters having almost the same size (the third cluster is far smaller) is a typical feature of percolation of Ising clusters (see footnote 1 page 107).

d (t) ∼ ξ(t) D ∼ t D /z d = √ t , ( 4
In the course of this process, the quench protocol went through the Ising critical point, and there is no track of it. Now the question is: what happens if, like in real experiments, we cannot do the quench instantaneously? What is the influence of the time spent near the Curie temperature? Section 5 will address these questions. However, let us first explore the dynamics after an instantaneous quench to the Curie temperature.

Quench to T = T c

The process is the same as above, except that the temperature immediately after t = 0 is now the Curie temperature:

T (t) := 2 T c t ≤ 0 T c t > 0 . ( 4

.2.1)

In this situation, the correlation length still grows as a power law,

ξ(t) ∼ t 1 /zc , (4.2.2)
with z c :≈ 2.17 the critical dynamical exponent [START_REF] Blanchard | A morphological study of cluster dynamics between critical points[END_REF][START_REF] Duclut | Frequency regulators for the nonperturbative renormalization group: A general study and the model A as a benchmark[END_REF][START_REF] Nightingale | Monte Carlo computation of correlation times of independent relaxation modes at criticality[END_REF][START_REF] Hohenberg | Theory of dynamic critical phenomena[END_REF]; the growth of the correlation length is slightly slower than in the previous situation since Finally, the average sizes of the largest clusters evolve in a very similar way to the one found in the T = 0 quenches: the only perceptible differences are that the smallest clusters do not disappear (thanks to the thermal fluctuations), and the dynamics are slightly slower (since z c > z d ). See the lower panels in figs. II.9 and II.11.

Effects of a finite cooling rate

In the present section, we shall discuss how the time spent in the vicinity of the Ising critical point affects the dynamics.

Let us first describe the process considered in the remainder of this presentation. The system is initially placed in an equilibrium state at T = 2 T c (ie T (t) = 2 T c for all t ≤ 0). Next, at t = 0, the temperature of the bath is linearly cooled following In the present study, we only consider the dynamics above T c (ie t ∈ [0, τ Q ]). Studies of the cooling rate effects on the coarsening dynamics that is at work close and below the critical point, even after annealing, have been presented in [START_REF] Biroli | Kibble-Zurek mechanism and infinitely slow annealing through critical points[END_REF] for the bi-dimensional Ising model, in [START_REF] Jelić | Quench dynamics of the 2d XY model[END_REF] for the bi-dimensional xy model, in [START_REF] Jain | Critical dynamics of classical systems under slow quench[END_REF] for a one-dimensional non-equilibrium lattice gas model with a phase transition between a fluid phase with homogeneously distributed particles and a jammed phase with a macroscopic hole cluster, and in [START_REF] Comaron | Dynamical critical exponents in drivendissipative quantum systems[END_REF][START_REF] Davis | Formation of Bose-Einstein condensates[END_REF][START_REF] Larcher | Critical Dynamics in Quenched 2D Atomic Gases[END_REF][START_REF] Proukakis | [END_REF] for time-dependent dissipative and stochastic Gross -Pitaeviskii models relevant to describe cold boson gases.

T (t) T c := 2 - t τ Q , ( 5
-τ = t τ Q ξ eq (τ ), ξ(t) ξ eq (τ ) ξ(t), τ Q = 10 3 = 10 2 = 10 1 t /τ Q = 0 t /τ Q = 0.5 t /τ Q = 0.9 t /τ Q = 1 τ Q = 10 3 τ Q = 10 2 τ Q = 10 1

The Kibble -Zurek mechanism

Starting from a thermal state, the system will follow the equilibrium conditions dictated by the changing environment as long as it can: ie up to a time, called t, when the time needed to thermalise becomes too long with respect to the relative rate of variation of temperature. Next, the system falls out-of-equilibrium and its further evolution will be discussed later. Obviously, the slower the cooling, the later the system will fall out-of-equilibrium. For an infinite system size, the time required to thermalise at the Ising critical point is infinite; it actually scales as L zc , and, unless cooling rates are scaled with the system size in a convenient way, the system will necessarily fall out-of-equilibrium at a certain point. Conversely, for finite-size systems, there exists a sufficiently slow cooling rate such that the system never goes out-of-equilibrium; we shall discuss this point in section 5.3.

We suppose the cooling to be sufficiently slow so that the system falls out-of-equilibrium only in a close vicinity of the critical point. On the one hand, in equilibrium, the correlation length depends on the distance from the critical point as τ -ν . On the other hand, close to T c , the dynamic correlation length is assumed to grow in time as

ξ(t) ∼ (t + ξ 0 zc ) 1 /zc
, where ξ 0 is the initial correlation length and some constant factor. After an instantaneous quench from a state with correlation length ξ(t = 0) = ξ 0 to a temperature at a distance τ > 0 from the critical point, the thermal state is reached after a time τ th (τ ) such that ξ t = τ th (τ ) ∼ ξ eq (τ ) ∼ τ -ν . Assuming that the instantaneous quench is performed from 2 T c and neglecting the microscopic correlation length at this temperature, ξ 0 = 0, we have τ th (τ ) ∼ τ -νzc . Now, following the argument proposed by Zurek [START_REF] Zurek | Cosmological experiments in superfluid helium?[END_REF], the system falls out-of-equilibrium at a time t, when the remaining time needed to reach T c , τ Q -t in the linear cooling procedure, becomes smaller than the time needed to thermalise at the current temperature T (the standard notation is such that the temperature and time at which the system falls out-of-equilibrium are noted by hats). Hence, we have

τ Q -t ∼ τ th (τ ) ∼ τ -νzc ∼ 1 - t τ Q -νzc ∼ τ Q - t τ Q -νzc
, where τ is the distance from the critical temperature at t. Therefore, the system falls out-of-equilibrium at

t = τ Q -τ Q νzc /(1+νzc) , (5.1.1)
where νz c/(1 + νz c ) ≈ 0.685 and another constant factor.

In many papers from the KZ literature dealing with the slow cooling of atomic systems the assumption is that, after t, the system remains frozen and correlations do not grow beyond the correlation length present when it has fallen out-of-equilibrium:

ξ := ξ t ∼ τ -ν ∼ τ Q νzc /(1+νzc) τ Q -ν ∼ τ Q ν /(1+νzc) , (5.1.2)
where ν /(1 + νz c ) ≈ 0.315. This, however, is not correct in coarsening systems as already discussed in [START_REF] Biroli | Kibble-Zurek mechanism and infinitely slow annealing through critical points[END_REF][START_REF] Comaron | Dynamical critical exponents in drivendissipative quantum systems[END_REF][START_REF] Jain | Critical dynamics of classical systems under slow quench[END_REF][START_REF] Jelić | Quench dynamics of the 2d XY model[END_REF], for example.

The out-of-equilibrium dynamics

At early times, such that t < t, the system evolves in equilibrium and the correlation length grows as the equilibrium one at the temperature reached at the measuring time Consequently, when reaching the critical point, the system is (at least) thermalised up to a scale s, as soon as the cooling is slower than τ ths Q which is such that

Q -νDc /(1+νzc) θ 2 ( , t = τ Q ) - 4 κ c 8 + κ c log τ Q = 10 1 ≈ 10
s ∼ ξ ∼ τ ths Q ν /(1+νzc) ie τ ths Q ∼ s zc+ 1 /ν . (5.3.2)
We recall that for an infinitely fast quench to T = T c , the scale s is thermalised after a time τ ths ∼ s zc .

(5.3.3)

Since z c + 1 /ν ≈ 3.17 > 2.17 ≈ z c , an instantaneous quench is more efficient than a linear cooling to create the structures of the Ising critical point; the time spent far from T c is not useful to develop the Ising criticality, the system develops, instead, the percolation one.

This feature can also be observed by looking at the average sizes of the largest clusters by comparing the lower panels in figs. II.11 and II.14. Indeed, on fig. II.11, at t = 10 5 ∼ τ th L=1024 /10, the second largest cluster has already started to be swallowed by the first one, and the third and fourth have almost reached their equilibrium average sizes. In contrast, on fig. II.14, at τ Q = 10 5 , all the largest clusters are still far from their equilibrium average sizes. Moreover, the first and second are still of the same order of magnitude.

These results confirm that the dynamics are affected by the Ising critical point only in its close vicinity, and the time spent far from it is not helpful to get closer to equilibrium at T c .

Dynamics before reaching the critical point

In the previous section, we have shown that the behaviour when reaching the critical point does not really rely on the exact out-of-equilibrium dynamics in the range t ∈ t, τ Q ; whether the system remains frozen or evolves like in a post-quench dynamics.

In this last section we try to clarify the situation.

Let us consider that the system's typical length continues to grow as after an instantaneous quench after t; the correlation length should then grow as

ξ(t) ∼ t -τ Q + τ Q νzc /(1+νzc) 1 /zc ∼ τ Q 1 /zc τ Q -1 /(1+νzc) -τ 1 /zc ,
and forgetting the dependence in the cooling rate, as ξ(t) ∼ (τ ) 1 /zc , (5.4.1)

where the factor has changed but is still positive. Thus, for τ (or τ small enough), the correlation length is almost constant, and the system seems to be frozen. Moreover, the shape described by eq. ( 5 Figure II.17 -Approach to the critical point; dependency in the cooling rate. On the figure we have represented the evolution of the slopes of the wav when approaching the critical point together with the equilibrium one. κ S is extracted from the slope of the wav at short curvilinear length scales and is expected to have the Ising criticality when reaching the Curie temperature. κ L is extracted from the slope of the wav at long curvilinear length scales and corresponds to the percolation criticality.

Let us now recall that the correlation length at the time or temperature at which the system falls out-of-equilibrium scales as

ξ ∼ τ Q ν /(1+νzc) .
This is only valid while τ Q is such that ξ ≤ ξ eq (τ = 0). Beyond this point, the correlation length saturates to ξ = ξeq := ξ eq (τ = 0) ∼ L; especially for an infinitely slow cooling (equilibrium). Doing the rescaling ξ(t) → ξ(t) ξ and τ → τ ξν (5.4.2)

(where ξ is taken as its saturation value for the equilibrium curve), the lower panel in fig. II.16 shows that the correlation length has a universal behaviour. This is in agreement with a power law growth of the correlation length in the out-of-equilibrium phase of the dynamics, as assumed by eq. (5.4.1). Nonetheless, universality disappears far from the critical point since the equilibrium correlation length is subject to nonalgebraic corrections in this region.

Let us finally discuss the "change in criticality", as measured by the evolution of the parameter κ in the course of the cooling process and compare it to the equilibrium (κ eq ). As done before, κ is extracted from the wav. Figure II.17 represents the evolution of the criticality on short (κ S ) and long (κ L ) curvilinear abscissa length scales. The long length scales have almost always the criticality of percolation (κ = 6) except for very slow cooling rates and in the vicinity of the Curie temperature where the criticality starts to be affected by the Ising critical point. On short length scales, the system is able to achieve the equilibration process, and the criticality corresponds to the equilibrium one discussed in section 3 and represented by κ eq . However, for the fastest cooling rates, eg τ Q = 10 1 , even the short scales cannot follow the equilibrium.

Conclusions

The purpose of this work was to study the influence of a slow annealing of the control parameter on the dynamics close to a second order critical point (between a symmetric and symmetry broken phase; here, for Ising models, the Z 2 symmetry). More precisely, we analysed the evolution of the geometric and scaling properties of the interfaces between domains, close and at the critical point.

In order to set the stage, we first studied the fractal properties of the interfaces in equilibrium in the disordered phase. The analysis of the wav allowed us to reach our first conclusion: in equilibrium at T > T c , the properties of the interfaces are the ones of critical percolation until a crossover length-scale that decreases with increasing temperature. A temperature dependent crossover towards critical Ising fractality at short-length scales arises in a close vicinity of the critical point, visible only below, say, T = 1.1 T c .

So far, the influence of critical percolation on the dynamics of the bi-dimensional Ising -Glauber model after instantaneous quenches from infinite temperature to the critical point [START_REF] Blanchard | A morphological study of cluster dynamics between critical points[END_REF] and below it [110,[START_REF] Blanchard | How soon after a zero-temperature quench is the fate of the Ising model sealed?[END_REF][START_REF] Blanchard | Critical percolation in the dynamics of the bidimensional ferromagnetic Ising model[END_REF][START_REF] Sicilia | Domain growth morphology in curvature-driven two-dimensional coarsening[END_REF] was studied. The equilibrium result just mentioned indicates that this critical percolation geometry is actually already encoded in high temperature equilibrium configurations.

Next, we recalled some basic features of the coarsening dynamics following an instantaneous quench from equilibrium at T = 2 T c both to zero and the Curie temperatures. On the one hand, we confirmed that correlation functions scale with a growing length that increases algebraically with time. On the other hand, we highlighted the nontrivial evolution of the geometry of the domains of parallel spins. The critical percolation geometry of the interfaces present in the initial state is progressively transformed, starting from the short scales, towards the one of the target temperature: smooth at zero temperature and critical Ising at T c .

We then explained the Kibble -Zurek mechanism [START_REF] Kibble | Topology of cosmic domains and strings[END_REF][START_REF] Zurek | Cosmological experiments in superfluid helium?[END_REF] allowing one to estimate when the system falls out-of-equilibrium while approaching a critical point from the symmetric phase with a finite speed. While Zurek assumes that the system freezes after falling out-of-equilibrium, following [START_REF] Biroli | Kibble-Zurek mechanism and infinitely slow annealing through critical points[END_REF][START_REF] Jelić | Quench dynamics of the 2d XY model[END_REF][START_REF] Proukakis | [END_REF] we argued that the correlation length continues to grow in this regime as if the system were instantaneously quenched to the critical point. Our argument does not affect the scaling in the cooling rate predicted by Zurek, but offers a more accurate description of the growing of the correlation length after the system has fallen out-of-equilibrium. We examined this scaling numerically and we found that after a slow linear cooling with rate τ Q to T c , the dynamic growing length extracted from the analysis of the space-time correlation function scales as ξ(τ Q ) ∼ τ ν /(1+νzc) Q

.

During the slow cooling process, while far from the critical point, the interfaces keep the fractal properties of critical percolation over a wide range of length scales, up to a temperature dependent crossover length. However, when approaching temperatures close enough to the critical point, we observe that the variance of the winding angle satisfies a scaling with respect to ξ(τ Q ), and besides, the interfaces with critical Ising properties span the length scales that are shorter than τ Q νDc /(1+νzc) .

Finally, our results prove that the Ising critical point influences the dynamics only in its close vicinity. Therefore, an instantaneous quench procedure is much more efficient to create the structures of the Curie temperature than a slow annealing. As a matter of fact, the time spent far from the Ising critical point does not contribute to the thermalisation of the system; instead, the dynamics of the system are governed by critical percolation.

This study is also a complement to works that try to elucidate the role played by the initial conditions on the post-quench dynamics of the Ginzburg -Landau scalar field theory [START_REF] Bray | Kinetics of ordering for correlated initial conditions[END_REF] and, more recently, the kinetic Ising model [START_REF] Chakraborty | 119, Fractality in persistence decay and domain growth during ferromagnetic ordering: Dependence upon initial correlation[END_REF]119,[START_REF] Corberi | Role of initial state and final quench temperature on aging properties in phase-ordering kinetics[END_REF] as well as the influence of a non-vanishing cooling rate on the scaling properties of discrete models close to their phase transition [START_REF] Liu | Dynamic scaling at classical phase transitions approached through nonequilibrium quenching[END_REF]. In the latter paper the emphasis was set on the scaling properties of the order parameter and how these depend, or not, on the microscopic stochastic updates. We focus instead on the geometrical and scaling properties of the structures when slowly approaching the critical point.

Bibliography

Chapter I -Numerical integration of classical conservative field theories

The content of this chapter has resulted in the following publication

Chapter II -Critical percolation in ferromagnetic Ising spin models

The content of this chapter has resulted in the following publication

List of Figures

I.1 Illustration of the tangent space to a 2-dimensional manifold. T z Ω is the tangent space to Ω at z. The basis vectors, ∂ 1 and ∂ 2 , can be seen as velocity vectors for a point moving respectively along the curvilinear coordinate z 1 or z 2 . z 1 , z 2 is a local coordinate system that parametrises Z, an open neighborhood of z on Ω. In this example, the coordinate system z 1 , z 2 is in fact globally defined on Ω.

• • • • • • • • • • • • • • • • • • I.2 Explicit Euler Implicit Euler -1 0 1 2 -2 -1 0 1 q p Partitioned Euler Midpoint rule
Deformation of the phase space under the approximate evolution map for the four methods considered here when applied to the harmonic oscillator problem. The left column regroups the explicit integrators while the ones on the right column are implicit; the integrators on the second row are symplectic (those on the first row are not). The evolution flow is represented in light grey and the volumes drawn using the same colour represent the exact deformation of the phase space under the Hamiltonian flow. By volumes we point out the regions of phase space, initially drawn as squares, but that may be distorted by the integration process. The filled volumes represent the successive deformation of the phase space obtained using the numerical integrators; a single time-step of length δ = T /8 is used to pass from a volume to the next one. We proceed as follow: we first discretise the volume in the phase space, then we independently perform a time-step of length δ = T /8 for each state constituting the initial volume; once done we reconstruct the transformed volume from these updated samples. The differences between the different integrators must not be attributed, here, to the fact they have a different order of approximation; indeed, for such a δ all the integrators produce errors of the same order of magnitude (the higher order of the Midpoint approximation distinguishes it from the other integrators only when δ goes to zero). Consequences in phase space of numerically approximating the evolution of the simple pendulum using each of the four methods considered in this section. See fig. I.2 for a description of the content; here the process only differs in that the numerical evolution between two successive volumes is obtained using two time-steps of length δ = T (q 0 , p 0 )/16 (where (q 0 , p 0 ) is the centre of the initial square). The effective distance in between two successive volumes is still T (q 0 , p 0 )/8, however performing such an evolution in only one step produces too large deviations with respect to the exact evolution; the comparison becomes impossible, especially for the less precise integrators.

• • • • • • • • • • • • • • • • • • • • I.3
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • I.8
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 81 I.15 x = j + 1+σn √ 2 δ t = n δ/ √ 2 δ δ j n
Lattice description of the space-time manifold used in the BDdV method.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 85
List of Figures I.17 Error committed by the different methods on the local conservation of the stress-energy tensor as a function of the initial amplitude (discretisation is still on 128 points). The red line is the exact solution of eq. (4.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 data-points represent the field evolved by the different numerical methods (diamonds for Euler, triangles for BDdV and circles for msilcc). On this scale we see no difference between the exact and the numerical solutions obtained with the Euler and msilcc methods.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • I.19
• • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Résumé

Cette thèse est constituée de deux parties indépendantes.

Dans le premier chapitre, nous introduisons une méthode numérique permettant d'intégrer des équations aux dérivées partielles représentant la dynamique Hamiltonienne de théories des champs. Cette méthode est un intégrateur multi-symplectique qui préserve localement le tenseur énergie-impulsion sur de très longues périodes de temps et avec précision. Son principal avantage est d'être extrêmement simple (c'est essentiellement un schéma à point médian) tout en restant bien définie localement. Nous la mettons à l'épreuve sur le cas particulier de l'équation d'ondes non-linéaire en une dimension d'espace (en considérant un terme potentiel quartique) ; nous expliquons également comment l'implémenter en dimensions supérieures. De plus, nous faisons une présentation géométrique du problème considéré, mettant ainsi en évidence la structure multi-symplectique sous-jacente. En outre, nous introduisons une construction permettant de résoudre le problème de dégénérescence pouvant affecter cette structure géométrique essentielle.

Le second chapitre traite d'aspects hors équilibre dans les systèmes statistiques : nous nous intéressons en particulier à la question de l'impact d'un taux de refroidissement fini lors d'une trempe à travers une transition de phase du second ordre. Pour décrire plus fidèlement le régime hors équilibre qui se produit avant la transition de phase, nous étendons le mécanisme dit de Kibble -Zurek qui décrit la dynamique d'un système statistique lors d'une telle procédure. Nous décrivons comment la taille caractéristique des objets géométriques présents dans le système dépend du temps, mais aussi comment elle dépend du taux de refroidissement ; ceci, au point critique et avant que celui-ci ne soit atteint. Ces prédictions théoriques sont mises en relief sur un exemple concret, via l'étude numérique d'un modèle emblématique : la version cinétique du modèle d'Ising ferromagnétique sur réseau carré. Nous faisons également une description des propriétés géométriques des domaines qui apparaissent dans le système au cours de la dynamique de refroidissement ; nous caractérisons également ces propriétés géométriques au moment particulier où l'on atteint la transition.

Figure I. 1

 1 illustrates this construction.

  2.5b) that become d b (∂ a ) := δ b a , and (1.2.8a) ∂ b (d a ) := δ a b . (1.2.8b) Let us now show how co-vectors behave under reparameterisation of the manifold. Again, let {z a } and {z a } be two curvilinear coordinate system on the manifold; let ∂ a := ∂ /∂z a and ∂a := ∂ /∂ za be the associated two basis in T Ω; and let {d a := dz a } and { da := dz a } be the corresponding two dual basis in T * Ω. The basis co-vectors generally transform as d a = T a b db , and (1.2.9a) da = T a b d b .

i.

  Let us start with the methods of the Runge -Kutta type. Let us first rewrite the equation of motion (1.4.1b) as ża = -ω ab ∂ b H =: χ a .(1.4.41)

Figure I. 4 -

 4 Figure I.4 -Error committed by the different methods on the solution (top graph)and on the energy (bottom graph) as a function of the time-step for the harmonic oscillator problem. The worst error encountered in the course of the numerical integration of the system over a period of the motion (0 ≤ t ≤ T ) is represented here. Both graphs have the same horizontal axis. The exact integrator stands for the exact evolution map recursively applied; hence, the obtained solution is exact in the sense that there are no errors due to the discretisation but it is, however, subject to the accumulation of truncation errors. This exact integrator represents the maximum possible accuracy and is directly related to the precision of the representation employed by the computer to encode real numbers (here double floating-point).

Figure I. 5 -

 5 Figure I.5 -Error committed by the different methods on the solution (top graph) and on the energy (bottom graph) as a function of the time-step for the simple pendulum problem. See fig.I.4 for a description of the content; here the exact integrator is not represented since its behaviour is unchanged with respect to the case of the harmonic oscillator. The exact integrator is not affected by the nonlinearities; only by the truncation errors that have no reason to be different in this situation. Moreover, it is extremely expensive in computation time since it requires to evaluate many Jacobi elliptic functions.

  Figure I.6 -Error committed in reduced precision by the different methods, on the solution as well as on the energy, as a function of the time-step for the harmonic oscillator problem. See fig.I.4 for a description of the content; here the machine precision has been artificially reduced to only 4 significant digits with the aim to feign the possibly significative accumulation of truncation errors that may occur with the very large number of steps necessary to reach long-times. For the sake of reducing the computation cost, this aspect is simulated here on much shorter times by amplifying the effects of the truncation. All the other numerical results presented in this thesis were obtained in double floating-point precision.

1 ∆

 1 Figure I.8 -long-time behaviour of the error committed by the different methods when applied to the simple pendulum problem. See fig. I.7 for a description of the content.The error committed by the explicit Euler method on the solution suddenly goes out of range after t/T (q 0 , p 0 ) ∼ 100.

Figure I. 9 -

 9 Figure I.9 -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.2.14a) to (2.2.14d).

Figure I. 10 -

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

.5. 5 )

 5 Now, we use the exact conservation of the multi-symplectic structure ϵ ωµ Ďν ζ, ζ ∂ν = 0 = ωµ ϵ Ďν ζ, ζ ∂ν ∂µ + ωµ Ďν ζ ∂µ , ϵ ζ ∂ν , and the dual of the equation of motion (2.1.6) to prove that ωµ ϵ Ďν ζ, ζ ∂ν ∂µ = ωµ ϵ ζ ∂ν , Ďν ζ ∂µ = dH ζ Ďν ζ .

  (3.4.2)), no longer coincides with the full average • • • (eq. (3.3.2a)).

Figure I. 11

 11 Figure I.11 represents the short-time behaviour 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. I.11, A = 3, the small circles at the centre 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.

  .1.5a) T 01 = T 10 := -∂ 0 φ ∂ 1 φ , and (4.1.5b)

Figure I. 11 -Figure I. 12 -

 1112 Figure I.11 -Space-time plots of the solutions of eq. (4.1.1) with initial and boundary conditions given by eq. (4.1.3) and eq. (4.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 colour is constant in between. Figure I.12 represents a cross-section of these space-time plots for the smallest values of A.

.2. 2 )Figure I. 13 -

 213 Figure I.13 -Lattice description of the space-time manifold in the Euler method.The lattice spacing in space and time are chosen to be equal to preserve covariance. This choice of an equal lattice spacing in space and time respects the CFL condition (computed from the von Neumann stability analysis) for the linear wave equation.

.2. 6 )

 6 The nodes of the lattice that appear in this equation are highlighted in fig. I.13: 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 .

Figure I. 14 -

 14 Figure I.14 -Time evolution of the total energy for A = 10 with the Euler integrator and L /δ = 128.

.2. 11 )

 11 These three definitions are represented as a function of time for A = 10 in fig. I.14.

.3. 11 )

 11 The nodes of the lattice that are involved in this equation are highlighted in fig. I.15: the vertex in the left hand side of the equation is represented as , and the ones in the right hand side are represented as .

Figure I. 16 -√ 2 δ = 128 .

 162128 Figure I.16 -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.

  .3.15) and we follow their time evolution independently. The numerical outcome is shown in fig. I.16.

Figure I. 17 -√ 2 δ = 128 .

 172128 Figure I.17 -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. I.16. 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].

1 Figure I. 18 -

 118 Figure I.18 -Error committed by the different methods on the local conservation of the stress-energy tensor as a function of the initial amplitude (discretisation is still on 128 points).

1 Figure I. 19 -

 119 Figure I.19 -Error committed by the different methods on the conservation of the charges as a function of the initial amplitude (discretisation is still on 128 points).

10 ∆

 10 Figure I.20 -long-time behaviour of the error committed by the different methods on the local conservation of the stress-energy tensor for A = 10 (discretisation is still on 128 points). The upper pair of curves are for the BDdV method (triangles), the intermediate ones for the Euler 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 I. 21 -Figure I. 22 -

 2122 Figure I.21 -long-time behaviour of the error committed by the different methods on the conservation of the charges for A = 10 (discretisation is still on 128 points). Same symbol convention and time axis as in fig. I.20.

Figure I. 23 -

 23 Figure I.23 -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. (4.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

Figure II. 1 -

 1 Figure II.1 -Drawing of a random configuration for the site percolation problem on Z 2 .Each site is occupied with probability p, and unoccupied with probability 1-p.In the represented configuration, and assuming periodic boundary conditions, the largest cluster percolates only if the missing site is occupied (probability p).

Figure II. 2

 2 Figure II.2 shows typical snapshots of site percolation for different values of the occupancy probability, and on a finite Z 2 lattice of size L = 128.

5 p 9 Figure II. 2 -

 592 Figure II.2 -Typical snapshots of the site percolation problem on a L × L = 128 2 square lattice with periodic boundary conditions; different values of the occupancy probability, p, are represented.

cFigure II. 3 -

 3 Figure II.3 -Two typical discrete realisations of a Schramm -Loewner evolution with κ = 3 on the left, and κ = 6 on the right. The two trajectories on the right were obtained with the same noise realisations as for the trajectories on the left.One can remark that some parts of these trajectories appear quite smooth, and therefore seem not to be fractal anymore. It is actually an issue of the method we have employed to generate these examples: the sampling is not uniform in the curvilinear abscissa. Hence, some samples are too distant, and it inevitably misses the details in between. In fact the fractality should be everywhere the same.

Figure II. 4 -

 4 Figure II.4 -Illustration of the Glauber dynamics at work in the kinetic Ising model.On the first line, the selected spin is flipped with probability 1 since that reduces the energy of the system. On the second line, the selected spin is now flipped with probability p = e -8K since that would produce four unfavourable links (-4K) instead of the four favourable ones already present (-4K). The smaller K (ie the higher the temperature of the bath), the more probably the selected spin will be flipped.

2 . 3 )Figure II. 5 -

 235 Figure II.5 -Illustration of the procedure employed to measure the winding angle along a closed bi-dimensional curve. The curve, C, in black thick line, is parameterised by x( ) = sin 2π and y( ) = cos 3π cos π + sin 2π /2, with the curvilinear abscissa ∈ [0, 1] (the curve rotates clockwise). The value of the local angle, η( ), is represented in yellow thick line on a graph which is bent such that it follows C. Here, we clearly see that η( = 1) = η( = 0) while θ( = 1) = θ 0 -2π = θ 0 = η( = 0) ie there is one clockwise loop per turn. All the angles are measured with respect to the horizontal axis.

Figure II. 6 -

 6 Figure II.6 -Equilibrium behaviour above the Curie point (T ≥ T c ). The figure shows the wav as a function of the curvilinear length on the interfaces, at different temperatures. The two straight lines, κ = 3 and κ = 6, are the expected slopes for the Ising and percolation universality classes, respectively.

Figure II. 7 - 1 τ

 71 Figure II.7 -Equilibrium behaviour above the Curie point (T ≥ T c ). The upper panel displays, as a function of T , the value of κ extracted from the slope of θ 2 (log ) at short length . The horizontal axis, the same as on the graphic below, is a logarithmic scale where we added the two extreme points, 0 and ∞. The values of κ corresponding to the two universality classes (Ising and percolation) are labeled on the graph. The lower panel shows the average occupancy rates of the first largest clusters when approaching the critical temperature.

Figure II. 9 -

 9 Figure II.9 -Out-of-equilibrium evolution in post-quench dynamics (from T = 2 T c to T = 0). The upper panel shows the evolution over time of the correlation length, ξ(t), extracted from the space-time two point correlation function. Its theoretical time-dependence is shown with a dashed line; the range of validity of this prediction is highlighted by the grey shading (a ξ(t) L). The lower panel represents, as a function of time, the average occupancy rates of the first largest clusters.

Figure II. 10 -

 10 Figure II.10 -Out-of-equilibrium evolution in post-quench dynamics (from T = 2 T cto T = 0). The upper panel represents the variance of the winding angle (wav) at different times following the quench. These times are reported in the upper panel of fig. II.9, and are chosen such that the constraint a ξ(t) L is fulfilled. The lower panel still represents the wav, but with a different scaling: the horizontal axis is rescaled following eq. (4.1.4), and d (t) is evaluated through its theoretical expression (∼ t 1 /z d ).

Figure II. 9 (

 9 Figure II.9 (lower panel) shows the evolution of the average size of the largest clusters.Starting from a high-temperature equilibrium state, all clusters are almost of the same size. Next, in the early dynamics, they all grow in the same way. As soon as the correlation length starts to grow, the larger clusters progressively swallow the smaller ones. Indeed, the smaller the clusters, the faster they disappear. This is the so-called coarsening dynamics. In particular, the second largest cluster lengthly coexists 1 with the largest one. As already mentioned, this long coexistence of two large clusters having almost the same size (the third cluster is far smaller) is a typical feature of percolation of Ising clusters (see footnote 1 page 107).

.0. 1 )

 1 where τ Q is the cooling time up to the Curie temperature (see fig. II.13).

Figure II. 11 -Figure II. 12 -

 1112 Figure II.11 -Out-of-equilibrium evolution in critical post-quench dynamics (temperature is instantaneously taken from T = 2 T c to T = T c ). Graphics are organised in the same manner as in fig. II.9. However, in the lower panel, we added with dashed lines the equilibrium values of the average occupancy rates of the first largest clusters at T c .

Figure II. 13 -

 13 Figure II.13 -The left panel in the upper row describes the cooling process. Temperature is linearly decreases from T = 2 T c at t = 0 to T = T c at t = τ Q . τ Q controls the cooling rate, and the larger the values it takes, the slower the cooling. The lower panel shows typical snapshots of the system (L = 128) in the course of the cooling process, and for different values of the cooling rate. The right panel in the upper row displays the evolution of the correlation length extracted from the space-time correlation function in the course of cooling in a system with L = 1024 (note that the maximum value of ξ is close to 10, much shorter than the system size). We have also represented the equilibrium correlation length at the corresponding temperatures.

Figure II. 14 -

 14 Figure II.14 -Dependence on the cooling rate after a linear cooling to the critical point. The upper panel shows, as a function of the cooling rate, the measured correlation length when reaching the critical point; the dashed line is its predicted evolution (see eq. (5.2.5a)). The second panel represents the average occupancy rates of the largest clusters, and the dashed lines highlight the expected values for an infinitely slow annealing (ie the values in equilibrium at T c ).

θ 2 (

 2 , t = τ Q )

5 Figure II. 15 -

 515 Figure II.15 -Dependence on the cooling rate after a linear cooling to the critical point. The upper panel represents the wav for different cooling rates, and the lower one shows the same quantity, after the rescaling proposed in eq. (5.3.1); κ c = 3 and D c = 1.375 are the same as in section 4.2.

  .4.1) is in a quite good agreement with the numerical results presented in fig. II.16 (upper panel).

Figure II. 16 -κ

 16 Figure II.16 -Approach to the critical point; dependency in the cooling rate. The upper panel shows the increase of the correlation length during the cooling for different cooling rates; the equilibrium correlation length is also shown. The lower panel represents the same quantities, but with a different scaling (following eq. (5.4.2)).

1 A 11 A√ 2 δ 77 I. 13 x = j δ t = n δ δδ

 11127713 H = H0) long-time behaviour of the error committed by the different methods when applied to the simple pendulum problem. See fig. I.7 for a description of the content. The error committed by the explicit Euler method on the solution suddenly goes out of range after t/T (q 0 , p 0 ) ∼ 100. 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.2.14a) to (2.2.14d). • • • • • • • • • • • • • • • • portion of the lattice M in D = 0 + 1,1 + 1 and 2 + 1 dimensions where we have highlighted the elementary cell.• • • • • • • • • • • • • • • • • • • 63 I., t)/ASpace-time plots of the solutions of eq. (4.1.1) with initial and boundary conditions given by eq. (4.1.3) and eq. (4.1.2) respectively. Different panels show data for different values of A, obtained with the msilcc method and L / = 1024. Lines are iso-levels of the field while colour is constant in between. Figure I.12 represents a cross-section of these space-time plots for the smallest values of A.• • • • • • • • • • • • • • • • • • • • • • • • • • •76 of the space-time plots in fig. I.11 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 . • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •Lattice description of the space-time manifold in the Euler method. The lattice spacing in space and time are chosen to be equal to preserve covariance. This choice of an equal lattice spacing in space and time respects the CFL condition (computed from the von Neumann stability analysis) for the linear wave equation.• • • • • • • • • • • • • • • • • • • • • • • • • • • 79Time evolution of the total energy for A = 10 with the Euler integrator and L /δ = 128.

√ 2 δ = 128 .

 2128 Time evolution of the total energy for A = 10 using the BDdV method with L / 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.

√ 2 δ = 128 .

 2128 Time evolution of the total energy for A = 10 using the msilcc method with L / The horizontal axis is the same as in fig. I.[START_REF] Orioli | Universal self-similar dynamics of relativistic and nonrelativistic field theories near nonthermal fixed points[END_REF]. 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].• • • • • • • • • •

1

 1 Error committed by the different methods on the conservation of the charges as a function of the initial amplitude (discretisation is still on 128 points). behaviour of the error committed by the different methods on the local conservation of the stress-energy tensor for A = 10 (discretisation is still on 128 points). The upper pair of curves are for the BDdV method (triangles), the intermediate ones for the Euler 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).• • • • • • • • • • • • • • • • I.21

  long-time behaviour of the error committed by the different methods on the conservation of the charges for A = 10 (discretisation is still on 128 points). Same symbol convention and time axis as in fig. I.20.

  ••| (for contravariant indices) and |•••| (for covariant indices), that exclude the enclosed indices from an outer symmetrisation/anti-symmetrisation operation.

  i , t that depends on the Darboux coordinates and explicitly on time. Its differential,

	df =	∂f ∂t	dt +	∂f ∂q i dq i +	∂f ∂p i	dp i ,	(1.3.21)
	df dt	=	∂f ∂t	+	∂f ∂q i	dq i dt	+	∂f ∂p i	dp i dt
		=	∂f ∂t	+	∂f ∂q i	qi +	∂f ∂p i	ṗi
		=	∂f ∂t	+	∂f ∂q i	∂H ∂p i	-	∂f ∂p i	∂H ∂q i
		=	∂f ∂t	+ {f, H} .
	Since this holds for any smooth function, one finally obtains
				d dt	=	∂ ∂t	+ {•, H} .	(1.3.22)
	Thus, the Poisson bracket allows to rewrite the set of Hamilton's eqs. (1.3.9a) to (1.3.9c) in a more symmetrical way
				qi = q i , H ,	(1.3.23a)
				ṗi = {p i , H} , and Ḣ = 0 .	(1.3.23b) (1.3.23c)

is hence composed of the exterior derivative on Ω (that takes into account the dependence on the Darboux coordinates) plus the part corresponding to the time dependence. Dividing by dt and identifying the Hamilton's eqs.

(1.3.9a

) and (1.3.9b) and then the definition of the Poisson bracket, one successively obtains Another consequence of eq. (1.3.22) is that any function whose Poisson bracket with H vanishes and that does not depend explicitly on time is a conserved quantity under Hamiltonian flow. Such quantities are called integrals of motion. Any linear combination of integrals of motion in an integral of motion too and, following the Jacobi identity

(1.3.20)

, the Poisson bracket of two integrals of motion is an integral of motion too. Hence, the set of all the integrals of motion is a Lie algebra.

  2 arccos dn t + τ k 2

	-2 sgn sn τ k 2 arccos dn τ k 2 p(t) = 2k cn t + τ k 2 , where	, and	(1.4.9a) (1.4.9b)
	τ = cn -1 p 0 2k	k 2 .		(1.4.9c)

  .1.3) since D := 1 is the Hausdorff dimension of the smooth interfaces on short length scales. The wav has a universal behaviour in time that is highlighted by the rescaling

	θ 2 ( , t) → θ 2 ( , t) and	→	d (t)	,	(4.1.4)

ie θ 2 ( , t) is a function of / d (t) only, once again, while ξ(t) is in the range a ξ(t) L. See the lower panel in fig. II.10.

  1 /z c ≈ 0.461 < 1 /2. See the upper panel in fig. II.11. Again, this result is only true for ξ(t) in the range a ξ(t) L.Concerning the wav, it behaves exactly as in the zero temperature quench except that, instead of the smooth zero temperature thermal state, it is the Ising criticality (κ = 3) that develops over short length-scales, see fig. II.12 (upper panel). The typical (curvilinear abscissa) length scale that separates the Ising criticality from the percolation one, c (t), scales now differently with the correlation length: := 3 is the same universal parameter as in the pre-factor in front of the logarithmic growth of the wav. Moreover, the Ising criticality grows faster than how the smoothness was growing in the zero temperature quench since D c/z c ≈ 0.634 > 1 /2. The wav still has a universal behaviour, now highlighted by the rescaling

	c (t) ∼ ξ(t) Dc ∼ t Dc /zc .			(4.2.3a)
	Since the interfaces on short length scales are not smooth anymore, their fractal Haus-dorff dimension is given by
	D c := 1 +	κ c 8	= 1.375 ,			(4.2.3b)
	where κ 4 κ c 8 + κ c	log	and	→	c (t)	,	(4.2.4)

c θ 2 ( , t) → θ 2 ( , t)where 4 κ c/(8 + κ c ) ≈ 1.

09, and still while ξ(t) is in the range a ξ(t) L; see the lower panel in fig. II.12.

in this chapter, we directly consider a local field theory; not an underlying model on a lattice whose field theory would be an approximate description.

most of the times, these indices are omitted; here we only wrote them to highlight on what the permutation apply.

the product and the bracket notations will be indifferently used.

or simply a metric.

here we assume that Ω is such that there exists a (0, 2)-tensor that is never degenerated (tensors can be defined on any differentiable manifold, but not all the differentiable maniflods admit a metric).

sometimes known as Clairaut's theorem.

a vector field is a continuous set of vectors from the tangent space for all the points on the manifold ie it is a section of the tangent bundle (again, the tangent bundle being the total space obtained by fibration of the manifold with its tangent space ie z∈Ω TzΩ).

one for each degree of freedom.

using the inverse of the Legendre transformation that were used to define the Hamiltonian (see eq. (1.3.6)).

therefore the Liouville's theorem is sometimes improperly presented as the conservation of the volume.

ie the Lie derivative with respect to the Hamiltonian vector field.

that is the fundamental geometric object over the phase space as we have shown earlier.

not perfectly since the preservation of the symplectic structure implies the exact conservation of the volume form, not of its integral.

or, at least, to be reducible to a flat manifold by a bijective transformation of space-time (at the cost of affecting the equations of motion).

we ignore the boundary terms since they produce constant terms in the action (ie independent of the variations of the fields), which do not affect the "position" of least action (ie the solutions of the Euler -Lagrange equations).

ie Schwarz's theorem.

in the 1 + 1-dimensional case, but the generalization to higher dimensions is straightforward.

in a time-scale that scales as, typically, a small power of the system size.

called a percolating cluster.

ie the correlation length.

here we are only interested in the behaviour above the Curie temperature (ie T > Tc).

in site percolation, at p = pc, the largest cluster (the percolating one) is much larger than the second one (of approximately one order of magnitude). In Ising models two percolating clusters are in competition: the up spins one and the down spins one; generally, in a q-state Potts model, the q largest clusters are of the same order of magnitude while the q + 1-th will be much smaller.

their size are of the same order of magnitude while the third largest cluster is far smaller.

Remerciements

For such a theory, the trace of the stress-energy tensor reduces to

and hence, the stress-energy tensor of a massless free field is indeed traceless (T = 0) in two dimensions (D = 1 + 1).

Let us now sketch why we have defined the stress-energy tensor as eq. (2.3.1). The following discussion is taken from [START_REF] Peter | Primordial cosmology[END_REF]. Suppose the action of our theory,

be invariant under the local transformation

where 1 is constant, while may depend on x. The local variation of the field is

while the global variation is

The local variation only depends on the nature of the field (scalar, vector, tensor, . . . ), while the global variation includes the effect of the local transformation on the functional form of the fields. They are related by

(2.3.9)

On the other hand, the measure, affected by a Jacobian, reads

Then, the variation of the action is

The global variation of the Lagrangian density reads

After some straightforward algebraic manipulations, we obtain

Equation (3.3.4) 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.1.6) in cell(n) at n C is

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 the time lattice M :

q n := q(t = n δ) , and

p n := p(t = n δ) .

(3.3.6b)

Then, applying the msilcc scheme, we get the discrete version of the equations of motion

, and (3.3.7a)

This is the Midpoint rule (see section 1.4). 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.

The λ φ 4 theory in 1 + 1 dimensions

The lattice is now the same as for the Boyanovsky -Destride Vega (BDdV) method (see [START_REF] Boyanovsky | Approach to thermalization in the classical ϕ 4 theory in 1 + 1 dimensions: Energy cascades and universal scaling[END_REF] or section 4.3).

Defining

we sample the state vector through the space-time lattice as

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.1.6):

and we obtain the set of skew-symmetric matrices in the new coordinate system

M j , and (3.4.9a)

From M µ , we define ωµ that actually behaves as a component of a D-vector in M:

The set { ωµ } defines the multi-symplectic structure in the light-cone coordinate system and is subject to the same conservation law Then, the local conservation of multi-symplecticity is, numerically,

since the contraction of the symmetric object ∂ a ∂ b with the skew-symmetry of the wedge product vanishes.

So, the msilcc scheme is a multi-symplectic integrator since it exactly preserves the discrete version of the conservation law of the multi-symplectic structure.

Conservation of the stress-energy tensor

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

mapped to (1, 0, • • • , 0) (in that way, each cell will only contain one unknown). So,

Then, it remains to orthogonalize the remaining rows of R, which can be achieved by defining

where µ ∈ 1, d and ν ∈ 0, d .

The space-time is now discretised using a hypercubic lattice rotated by R, and the discretisation rules are simple concatenations of Midpoint rules as introduced in section 3.1 (a centred box scheme). All the proofs exposed throughout this section are still valid for such a lattice (with some adjustments), whatever its orientation. This alternative approach has therefore the same remarkable properties.

At the initial time, the first D -1 integration steps pose the same problem as in the previous version of the lattice: some nodes of the cell are at negative times. This issue is addressed by progressively increasing the dimension: by applying the centred box scheme in one dimension for the first step; in two dimensions for the second one; and so on and so forth, up to the D th step, for which the entire cell is now accessible. This still requires having a non-degenerate multi-symplectic structure.

Unfortunately, this modification of the method does not address the stability issues; the fastest modes are still unstable. The problem being that the truncation errors produce a noise, with no spacial correlation; inevitably leading to populate the fastest modes and, eventually, to a divergence of the solution.

This instabilities on short distances are not yet understood and an in depth investigation must be done to identify the origin of the problem. Nonetheless, the msilcc method remains functional in D ≤ 1+1 and provides accurate solutions as we highlight in the next section through an applicative example.

A possible way of improvement would be to consider the first generalisation of the lattice we proposed, but now with a different lattice spacing in space and time, as proposed in [START_REF] Mclachlan | The multisymplectic diamond scheme[END_REF] for the 1 + 1 dimensional case. Of course, one have to carefully analyse the impact that explicitly breaking the covariance will have on the quality of the method. If successful, a promising second way of improvement could be to generalise on this lattice the high-order integrators introduced in [START_REF] Mclachlan | The multisymplectic diamond scheme[END_REF].

4 Application: the λ φ 4 theory in 1 + 1 dimensions

This last 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 partitioned Euler approximation, 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

15 -Lattice description of the space-time manifold used in the BDdV method.

the energies is directly due to the divergence of the solution which has been destabilised by the integration method.

We conclude that, although Euler'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, Euler'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.

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 stress-energy 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

where

M is represented in fig. I.15 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 Euler method, the field, φ : M → R, can be sampled through the lattice as ϕ : M → R, such that

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

, and (4.3.4a)

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 previous one since the discretisation rules are not applied to the equation of motion but to the energy; 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.

The local energy density (the 00 component of the stress-energy tensor) is approximated in two ways (following the same principle as for the two possible Euler approximations of the derivatives):

where the n and j indices can no longer be omitted since they are not obvious. The difference between these two possible definitions reads

where

On the other hand, the total energy is given by

It can be shown (using periodic boundary conditions) that these two definitions are equivalent,

(τ (t)):

When t exceeds t, the correlation length does not grow fast enough, and the system falls out-of-equilibrium. Now, since the equilibrium correlation length (the expected one at the slowly varying temperature) soon becomes much longer than the dynamic growing length, we can make the simple assumption that the system behaves as if it were in contact with a bath right at the Curie temperature. This proposal amounts to treating the problem as after an instantaneous quench at t = t from τ = τ to τ = 0. Hence, the correlation length is assumed to continue to grow after t, but now as

where the term in ξzc takes into account the non-vanishing correlation length at t = t.

Then, imposing the consistency of the correlation length before and after t,

we have

where the second line is obtained from eq. (5.2.2), where we have replaced

In particular, when reaching the critical point,

.2.5b)

While Zurek assumes that the system is frozen immediately after falling out-ofequilibrium, here we claim that the dynamic growing length, ξ(t), continues to grow after t. However, the growth between t and τ Q , when the cooling reaches T c , only affects the pre-factor and not the scaling with τ Q that is not modified (see eq. (5.2.5b)). Therefore, if the interest is set upon the scaling properties of the system at the critical point (and not far below it) one can assume that the dynamic correlation length takes the form it had at t.

The next section will be devoted to putting this prediction to the test.

behaviour at the Ising critical point

We are now going to describe the state of the system when reaching the critical point, and how it depends on the cooling rate. First of all, we can easily check that eq. (5. Error committed by the different methods on the solution (top graph) and on the energy (bottom graph) as a function of the time-step for the harmonic oscillator problem. The worst error encountered in the course of the numerical integration of the system over a period of the motion (0 ≤ t ≤ T ) is represented here. Both graphs have the same horizontal axis. The exact integrator stands for the exact evolution map recursively applied; hence, the obtained solution is exact in the sense that there are no errors due to the discretisation but it is, however, subject to the accumulation of truncation errors. This exact integrator represents the maximum possible accuracy and is directly related to the precision of the representation employed by the computer to encode real numbers (here double floating-point). Error committed by the different methods on the solution (top graph) and on the energy (bottom graph) as a function of the time-step for the simple pendulum problem. See fig. I.4 for a description of the content; here the exact integrator is not represented since its behaviour is unchanged with respect to the case of the harmonic oscillator. The exact integrator is not affected by the nonlinearities; only by the truncation errors that have no reason to be different in this situation. Moreover, it is extremely expensive in computation time since it requires to evaluate many Jacobi elliptic functions. Error committed in reduced precision by the different methods, on the solution as well as on the energy, as a function of the time-step for the harmonic oscillator problem. See fig. I.4 for a description of the content; here the machine precision has been artificially reduced to only 4 significant digits with the aim to feign the possibly significative accumulation of truncation errors that may occur with the very large number of steps necessary to reach long-times. For the sake of reducing the computation cost, this aspect is simulated here on much shorter times by amplifying the effects of the truncation. All the other numerical results presented in this thesis were obtained in double floating-point precision. and the instantaneous error (closed symbols). Both graphs have the same horizontal axis. Beyond t = T the horizontal axis is shown in a different linear scale and the curves with closed 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). In order to fairly compare the different methods, δ is chosen independently for each of them such that they all produce the same error on the solution (10 -2 ) at t = T . In every case, δ is too small with respect to the time-scales considered here; there are only a few points represented on the graph, however, all the steps are taken into account for the worst encountered error.

List of Figures I.23 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. ( 4.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 settled to 1.

Drawing of a random configuration for the site percolation problem on Z 2 . Each site is occupied with probability p, and unoccupied with probability 1p. In the represented configuration, and assuming periodic boundary conditions, the largest cluster percolates only if the missing site is occupied (probability p).

Two typical discrete realisations of a Schramm -Loewner evolution with κ = 3 on the left, and κ = 6 on the right. The two trajectories on the right were obtained with the same noise realisations as for the trajectories on the left. One can remark that some parts of these trajectories appear quite smooth, and therefore seem not to be fractal anymore. It is actually an issue of the method we have employed to generate these examples: the sampling is not uniform in the curvilinear abscissa. Hence, some samples are too distant, and it inevitably misses the details in between. In fact the fractality should be everywhere the same.

Illustration of the Glauber dynamics at work in the kinetic Ising model. On the first line, the selected spin is flipped with probability 1 since that reduces the energy of the system. On the second line, the selected spin is now flipped with probability p = e -8K since that would produce four unfavourable links (-4K) instead of the four favourable ones already present (-4K). The smaller K (ie the higher the temperature of the bath), the more probably the selected spin will be flipped.

Illustration of the procedure employed to measure the winding angle along a closed bi-dimensional curve. The curve, C, in black thick line, is parameterised by x( ) = sin 2π and y( ) = cos 3π cos π + sin 2π /2, with the curvilinear abscissa ∈ [0, 1] (the curve rotates clockwise). The value of the local angle, η( ), is represented in yellow thick line on a graph which is bent such that it follows C. Here, we clearly see that η( = 1) = η( = 0) while θ( = 1) = θ 0 -2π = θ 0 = η( = 0) ie there is one clockwise loop per turn. All the angles are measured with respect to the horizontal axis.

Equilibrium behaviour above the Curie point (T ≥ T c ). The figure shows typical snapshots of the equilibrium state of the system (L = 128), at different temperatures above Out-of-equilibrium evolution in post-quench dynamics (from T = 2 T c to T = 0). The upper panel shows the evolution over time of the correlation length, ξ(t), extracted from the space-time two point correlation function. Its theoretical time-dependence is shown with a dashed line; the range of validity of this prediction is highlighted by the grey shading (a ξ(t) L). The lower panel represents, as a function of time, the average occupancy rates of the first largest clusters. Out-of-equilibrium evolution in post-quench dynamics (from T = 2 T c to T = 0). The upper panel represents the variance of the winding angle (wav) at different times following the quench. These times are reported in the upper panel of fig. II.9, and are chosen such that the constraint a ξ(t) L is fulfilled. The lower panel still represents the wav, but with a different scaling: the horizontal axis is rescaled following eq. (4.1.4), and d (t) is evaluated through its theoretical expression (∼ t 1 /z d ). Out-of-equilibrium evolution in critical post-quench dynamics (temperature is instantaneously taken from T = 2 T c to T = T c ). Graphics are organised in the same manner as in fig. II.9. However, in the lower panel, we added with dashed lines the equilibrium values of the average occupancy rates of the first largest clusters at T c . Out-of-equilibrium evolution in critical post-quench dynamics (temperature is instantaneously taken from T = 2 T c to T = T c ). Graphics are organised in the same manner as in fig. II.10. However, the upper panel is now scaled following eq. (4.2.4), where c (t) is evaluated through its theoretical expression (∼ t Dc /zc ).

List Dependence on the cooling rate after a linear cooling to the critical point.

The upper panel shows, as a function of the cooling rate, the measured correlation length when reaching the critical point; the dashed line is its predicted evolution (see eq. ( 5.2.5a)). The second panel represents the average occupancy rates of the largest clusters, and the dashed lines highlight the expected values for an infinitely slow annealing (ie the values in equilibrium at T c ). Approach to the critical point; dependency in the cooling rate. The upper panel shows the increase of the correlation length during the cooling for different cooling rates; the equilibrium correlation length is also shown. The lower panel represents the same quantities, but with a different scaling (following eq. (5.4.2)).

II.17 Approach to the critical point; dependency in the cooling rate. On the figure we have represented the evolution of the slopes of the wav when approaching the critical point together with the equilibrium one. κ S is extracted from the slope of the wav at short curvilinear length scales and is expected to have the Ising criticality when reaching the Curie temperature. κ L is extracted from the slope of the wav at long curvilinear length scales and corresponds to the percolation criticality.