Iain Henderson
email: henderso@insa-toulouse.fr

Pascal Noble

Olivier Roustant

Characterization of the second order random fields subject to linear distributional PDE constraints

Keywords: Generalized functions, Linear constraints, Linear Partial Differential Equations, Second order random fields

Let L be a linear differential operator acting on functions defined over an open set D Ă R d . In this article, we characterize the measurable second order random fields U " pU pxqqxPD whose sample paths all verify the partial differential equation (PDE) Lpuq " 0, solely in terms of their first two moments. When compared to previous similar results, the novelty lies in that the equality Lpuq " 0 is understood in the sense of distributions, which is a powerful functional analysis framework mostly designed to study linear PDEs. This framework enables to reduce to the minimum the required differentiability assumptions over the first two moments of pU pxqqxPD as well as over its sample paths in order to make sense of the PDE LpUωq " 0. In view of Gaussian process regression (GPR) applications, we show that when pU pxqqxPD is a Gaussian process (GP), the sample paths of pU pxqqxPD conditioned on pointwise observations still verify the constraint Lpuq " 0 in the distributional sense. We finish by deriving a simple but instructive example, a GP model for the 3D linear wave equation, for which our theorem is applicable and where the previous results from the literature do not apply in general.

Introduction

When dealing with an unknown function of interest u : D Ñ R where say D Ă R d , it is common (as e.g. in Bayesian inverse problems) to assume that it is a sample path of a random field U " pU pxqq xPD . Incorporating prior knowledge over u, such as smoothness, is then achieved by constraining the law of U accordingly. Sometimes, this prior knowledge comes from physical considerations. If u describes a positive quantity such as mass or energy, then the random variables U pxq should all be positive almost surely (a.s.). In many cases, this physical constraint can be more precisely translated as a partial differential equation (PDE). Such equations are a pivotal tool for modelling, understanding and predicting real-life phenomena such as those arising from fluid mechanics, electromagnetics or biology to name a few. The most simple (yet central) PDEs are those that are linear. In this article, we will only consider homogeneous linear PDEs, which take the form Lpuq :"

ÿ |α|ďn a α pxqB α u " 0. (1)
Above, u is the unknown function of interest, defined over an open set D Ă R d , and L is a linear partial differential operator. In [START_REF] Adler | The Geometry of Random Fields[END_REF], for a multi-index α " pα 1 , ..., α d q T P N d , we used the notations |α| " α 1 `... `αd and B α " pB x1 q α1 ...pB x d q α d . Homogeneous PDEs, i.e. PDEs with a null term on the right-hand side of [START_REF] Adler | The Geometry of Random Fields[END_REF], are often encountered to describe conservation laws, such as conservation of mass, energy or momentum in closed systems [START_REF] Serre | Systems of conservation laws[END_REF].

In order to incorporate the knowledge that Lpuq " 0 in the prior U , a natural question is whether one can characterize, in terms of their law, the random fields U whose sample paths are all solutions to the PDE [START_REF] Adler | The Geometry of Random Fields[END_REF]. Let U be a centered second order random field with covariance function k: under the assumption that U is a Gaussian process (GP) whose sample paths are n times differentiable, [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF] proved for some classes of differential operators L of order n that ( [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF], Sections 3.3 and 4.1)

PpLpU q " 0q " 1 ðñ @x P D, Lpkpx, ¨qq " 0.

(

) 2
This property provides a simple characterization of the GPs that incorporate the PDE constraint (1) sample path-wise. Such GPs would fall in the category of "physics-informed" GPs in the machine learning community. In the proof of this property, the fact that the sample paths are n times differentiable, i.e. that the PDE (1) can be understood pointwise, is central. These functions are then strong solutions of the PDE (1) (see Definition 3.1).

In the standard PDE approach though, equation ( 1) is reinterpreted by weakening the definition of the derivatives of u, thereby weakening the required regularity assumptions over u. It can indeed happen in practice that the sought solutions of the PDE Lpuq " 0 are not n times differentiable or even continuous (see e.g. [START_REF] Evans | Partial differential equations[END_REF], Section 2.1), and they are only solutions of some weakened formulation of equation [START_REF] Adler | The Geometry of Random Fields[END_REF]. This is typically the case for hyperbolic PDEs such as the wave equation presented in Section 4. We introduce here the distributional formulation of the PDE [START_REF] Adler | The Geometry of Random Fields[END_REF], where the regularity assumptions over u are relaxed to the maximum. As such, this formulation enables working with potentially singular solutions of equation [START_REF] Adler | The Geometry of Random Fields[END_REF], solutions which are not allowed to appear in more restrictive functional frameworks (see also the upcoming Remark 3.4). Another advantage of the distributional formulation is that it provides a unifying framework for dealing with linear PDEs, independently of their nature. In contrast, traditional weak or variational formulations vary greatly depending on the nature of the PDE. As an illustration, in [START_REF] Evans | Partial differential equations[END_REF], one can compare the different function spaces for weak solutions associated to elliptic PDEs (Section 6.1.2), parabolic PDEs (Section 7.1.1(b)) and hyperbolic PDEs (Section 7.2.1(b)). The distributional formulation will be our main object of interest in this article, and can be seen as a weakened form of weak formulations of PDEs. Consider equation [START_REF] Adler | The Geometry of Random Fields[END_REF], and "test it locally": that is, multiply it by a compactly supported, infinitely differentiable test function φ (i.e. φ P C 8 c pDq) and integrate over D:

@φ P C 8 c pDq, ÿ |α|ďn ż D φpxqa α pxqB α upxqdx " 0. (3)
For each integral term above, perform |α| successive integrations by parts to transfer the derivatives from u to φ. Since φ is identically null on a neighbourhood of the boundary of D, the boundary terms of each integration by parts vanish and we obtain that

@φ P C 8 c pDq, ż D upxq ÿ |α|ďn p´1q |α| B α pa α φqpxqdx " 0. ( 4
)
To make sense of (4), one only requires u to be locally integrable, i.e. ş K |upxq|dx ă `8 for all compact set K Ă D. We then say that a locally integrable function u is a solution to Lpuq " 0 in the sense of distributions, or distributional sense, if u verifies [START_REF] Álvarez | Linear latent force models using Gaussian processes[END_REF]. In this case, u is a solution to equation [START_REF] Adler | The Geometry of Random Fields[END_REF] in the sense of "all smooth local averages" (i.e. for all φ P C 8 c pDq), though not pointwise in general: taking φpxq " δ 0 px ´x0 q is not allowed without additional assumptions over u.

The distributional formulation of the PDE Lpuq " 0 is "compliant with physics" too, as pointed out by W. Rudin ([46], p. 150): most of the sensors we use in practice are only capable of computing local averages of the physical quantity they are measuring. Suppose one wishes to check experimentally that a temperature field obeys the heat equation, by using a set of thermometers: then one will actually only deal with the distributional formulation of the heat equation.

The natural question that follows from this new definition is whether one can characterize, in terms of their law, the random fields whose sample paths are solutions to the PDE Lpuq " 0 in the distributional sense. The answer is yes, and is the main content of this article. Under the assumptions that U is a measurable centered second order random field and that its standard deviation function σ : x Þ ÝÑ a kpx, xq is locally integrable, we show in Proposition 3.5 that PpLpU q " 0 in the distrib. senseq " 1 ðñ @x P D, Lpkpx, ¨qq " 0 in the distrib. sense. ( 5)

Related literature

It is known, at least since the fifties, that some covariance functions are naturally linked to certain stochastic partial derivative equations (SPDEs), i.e. PDEs where the source term is random. For example, it was already observed in 1954 by [START_REF] Whittle | On stationary processes in the plane[END_REF] that the covariance function of a stationary GP U verifying the two dimensional SPDE pα 2 ´B2 xx ´B2 yy q 3{4 U " W S , where W S is a spatial white noise process, is exponential, i.e. of the form CovpU px `hq, U pxqq " C expp´α|h|q. Already for this SPDE, the differentiation has to be understood in a weakened sense as white noise processes are not random fields in the usual sense. In [START_REF] Roques | Spatial statistics and stochastic partial differential equations: A mechanistic viewpoint[END_REF], SPDEs describing the random motion of micro-particles are introduced to link certain covariance functions, Matérn in particular, with an underlying physical model. We also refer to [START_REF] Lindgren | The spde approach for gaussian and non-gaussian fields: 10 years and still running[END_REF] for a large overview of the possible applications and recent developments pertaining to random fields defined by SPDEs. A general framework for the study of SPDEs was recently reintroduced in [START_REF] Carrizo-Vergara | A general framework for SPDE-based stationary random fields[END_REF], which was then used to classify the stationary generalized random fields that are solutions of a wide class of linear SPDEs. In particular, [START_REF] Carrizo-Vergara | A general framework for SPDE-based stationary random fields[END_REF] provides a description of all the second order stationary generalized random fields that are solutions to certain homogeneous PDEs, and the 3D wave equation in particular (which we also study in Section 4), in terms of their covariance operator. Loosely speaking, generalized random fields are function-indexed random fields where the covariance function is replaced by a covariance operator. From a functional analysis point of view, this is actually very close to the tools we use here, although in this article we constrain ourselves to work with (standard) random fields with well-defined sample paths, as these are the objects that arise the most in the random function models met in practice. The two other key differences between this work and [START_REF] Carrizo-Vergara | A general framework for SPDE-based stationary random fields[END_REF] are that piq we do not focus on stationary random field models for u and piiq we focus on the homogeneous case for PDE Lpuq " 0.

The literature concerning random fields that are PDE-constrained at the level of the sample paths is rather sparse. In [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF], general theorems are exposed for many different classes of linear operators acting on suitable spaces of functions. These theorems take the form of equation ( 2), and can in turn be applied to certain differential operators (see [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF], Sections 3.3 and 4.1). [START_REF] Scheuerer | Covariance models for divergence-free and curl-free random vector fields[END_REF] builds covariance functions that ensure that the sample path of a given two or three dimensional random vector field are either divergence or curl free. This result is notable because "any" three dimensional vector field can be decomposed as a sum of divergence and curl free vector fields through the Helmholtz-Hodge decomposition theorem. Moreover, divergence or curl free vector fields are commonly encountered in fluid mechanics. [START_REF] Fiedler | Distances, Gegenbauer expansions, curls, and dimples: On dependence measures for random fields[END_REF] extends the results of [START_REF] Scheuerer | Covariance models for divergence-free and curl-free random vector fields[END_REF] to random fields on the sphere of R 3 , which has been rediscovered later in [START_REF] Fan | Modeling tangential vector fields on a sphere[END_REF]. In [START_REF] Estrade | Anisotropic gaussian wave models[END_REF], stationary GPs are represented in terms of a random wavevector. [START_REF] Estrade | Anisotropic gaussian wave models[END_REF] then characterizes the stationary GPs whose sample paths verify a homogeneous linear PDE, in terms of the spectral measure of the GP and in terms of its random wavevector. [START_REF] Estrade | Anisotropic gaussian wave models[END_REF] additionally requires that the sample paths be infinitely differentiable, that the PDE's coefficients be constant and that only even orders of differentiation appear in the PDE. This is then applied to a few wave models. A simple algorithm for building linearly constrained GPs is proposed in [START_REF] Jidling | Linearly constrained Gaussian processes[END_REF], based on formal GPR derivations upon (1); however, partly because the assumed regularity of u is not fully addressed, the claim that the sample paths of the underlying GP are indeed linearly constrained is left unproved. This is clarified in [START_REF] Lange-Hegermann | Algorithmic linearly constrained Gaussian processes[END_REF], where the requirement that u P C 8 pR d q is made explicit and the enforcement of the PDE on the sample paths is proved for GPs whose sample paths are smooth. The algorithm from [START_REF] Jidling | Linearly constrained Gaussian processes[END_REF] is then supplemented in [START_REF] Lange-Hegermann | Algorithmic linearly constrained Gaussian processes[END_REF], where parametrizations of the solution spaces of (1) thanks to Gröbner bases are proposed. In [START_REF] Lange-Hegermann | Linearly constrained Gaussian processes with boundary conditions[END_REF], the same author completes the approach from [START_REF] Lange-Hegermann | Algorithmic linearly constrained Gaussian processes[END_REF] by incorporating boundary conditions on hypersurfaces in the Gröbner basis parametrization. With the idea to apply GPR to rigid body dynamics, [START_REF] Geist | Learning constrained dynamics with gauss' principle adhering gaussian processes[END_REF] enforces Gauss' principle of least constraint on the sample paths of a GP.

One can understand our main result (Proposition 3.5) as a characterization of the "physicsinformed" random fields that incorporate the distributional PDE constraint Lpuq " 0 at the level of the sample paths. It turns out that the design of similar "physics-informed priors" has received a lot a attention from the machine learning community since the early 2000' ( [START_REF] Graepel | Solving noisy linear operator equations by gaussian processes: Application to ordinary and partial differential equations[END_REF]), in the context of Gaussian process regression (GPR); see Section 3.3.1 for a description of this technique. GPR is a Bayesian framework for function regression and interpolation which is well suited for handling linear constraints, partly because GPs are "stable under linear combinations", see Section 2.1.3. The recurring idea is to assume that the function u in equation ( 1) is a sample path of a (centered) GP U and to draw the consequences of equation ( 1) on the covariance function of U . The covariance function of U is then expected to incorporate the constraint Lpuq " 0 in some sense. The majority of these works (except those mentioned above) do not aim at analysing whether the obtained covariance function indeed yields sample path PDE constraints over U : this is justified by the fact that they are only concerned with imposing the constraints on the function provided by GPR to approximate u. This approximation of u, which we denote by m, is called the Kriging mean in the GPR context; see equation [START_REF] Lange-Hegermann | Algorithmic linearly constrained Gaussian processes[END_REF] for a definition.

While they do not primarily focus on investigating sample path PDE constraints (contrarily to this article), the works coming from the GPR community are still very connected to this article. Indeed, they are concerned with designing explicit covariance functions that verify constraints of the form Lpkp¨, xqq " 0 for all x P D (the PDE is understood in the strong sense in these works). Indeed, this constraint ensures that all the possible regression functions m provided by the corresponding GPR model verify the constraint Lp mq " 0 (as seen in equation ( 33)). Note that "Lpkp¨, xqq " 0 @x P D" is the right-hand side of equation ( 2): actually knowing covariance functions that verify this constraint is a necessary complement to the condition we prove in this article (Proposition 3.5 is otherwise useless in practice). Explicit PDE constrained covariance functions were designed for a number of classical PDEs, namely: divergence-free vector fields [START_REF] Narcowich | Generalized Hermite interpolation via matrix-valued conditionally positive definite functions[END_REF][START_REF] Scheuerer | Covariance models for divergence-free and curl-free random vector fields[END_REF], curl-free vector fields [START_REF] Fuselier | Refined error estimates for matrix-valued radial basis functions[END_REF][START_REF] Scheuerer | Covariance models for divergence-free and curl-free random vector fields[END_REF], the Laplace equation [START_REF] Schaback | Solving the Laplace equation by meshless collocation using harmonic kernels[END_REF][START_REF] Mendes | Bayesian inference in the numerical solution of Laplace's equation[END_REF][START_REF] Albert | Gaussian process regression for data fulfilling linear differential equations with localized sources[END_REF]], Maxwell's equations [START_REF] Wahlstrom | Modeling magnetic fields using Gaussian processes[END_REF][START_REF] Jidling | Linearly constrained Gaussian processes[END_REF][START_REF] Lange-Hegermann | Algorithmic linearly constrained Gaussian processes[END_REF] (although [START_REF] Wahlstrom | Modeling magnetic fields using Gaussian processes[END_REF][START_REF] Jidling | Linearly constrained Gaussian processes[END_REF] only exploit curl/divergence free constraints), the 1D heat equation [START_REF] Albert | Gaussian process regression for data fulfilling linear differential equations with localized sources[END_REF], Helmholtz' 2D equations [START_REF] Albert | Gaussian process regression for data fulfilling linear differential equations with localized sources[END_REF], and linear solid mechanics [START_REF] Jidling | Probabilistic modelling and reconstruction of strain[END_REF]. [START_REF] Lange-Hegermann | Linearly constrained Gaussian processes with boundary conditions[END_REF] and [START_REF] Gulian | Gaussian process regression constrained by boundary value problems[END_REF] enforce homogeneous boundary conditions on the covariance function.

We finish with a brief overview of the alternative "physics-informed" GPR models. Contrarily to the equation ( 1) considered here, one may put a random source term f in the PDE and study instead the SPDE Lpuq " f : see [START_REF] Raissi | Machine learning of linear differential equations using Gaussian processes[END_REF], [START_REF] Álvarez | Linear latent force models using Gaussian processes[END_REF] and [START_REF] Owhadi | Bayesian numerical homogenization[END_REF] for entry points on the related literature. A recent article [START_REF] Chen | Solving and learning nonlinear PDEs with Gaussian processes[END_REF] extended the use of GPR to nonlinear PDEs by imposing the nonlinear interpolation constraints on the collocation points, setting the way forward for many possible applications of GPR to nonlinear realistic PDE models, as found e.g. in fluid mechanics. In [START_REF] Nguyen | Gaussian functional regression for linear partial differential equations[END_REF], the variational formulation (see [START_REF] Evans | Partial differential equations[END_REF], Section 6.1.2 for a definition) of certain linear PDEs has been incorporated into a GPR framework. This approach requires the use of Gaussian generalized random fields (see [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF], Section 2.2.1.1), or "functional Gaussian processes" following [START_REF] Nguyen | Gaussian functional regression for linear partial differential equations[END_REF]. The variational formulation of a PDE differs from its distributional formulation in the choice of the space of test functions.

Contribution and organisation of the paper. Consider the PDE in equation [START_REF] Adler | The Geometry of Random Fields[END_REF], where the coefficients of the differential operator L have possibly limited smoothness. Consider also a centered second order measurable stochastic process U " pU pxqq xPD with covariance function kpx, x 1 q (see Sections 2.1.2 and 2.1.3). Under the assumption that its standard deviation function σ : x Þ ÝÑ kpx, xq 1{2 is locally integrable, we show in Proposition 3.5 that the announced equation ( 5) holds. The result is then compared to a previous result from [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF], which ensures pointwise linear differential degeneracy of the sample paths of U under stronger assumptions. We then provide a simple corollary which states that linear distributional differential constraints are preserved when a GP U is conditioned on pointwise observations, in view of GPR applications.

As an application example, we derive a general Gaussian process model for the homogeneous 3D free space wave equation, for which the solutions are not smooth in general. This equation is central for describing finite speed propagation phenomena as found e.g. in acoustics. Plugging this model in a GPR framework yields potential applications in different inverse problems related to this PDE, such as thermoacoustic tomography (i.e. initial condition reconstruction, [START_REF] Kuchment | Mathematics of photoacoustic and thermoacoustic tomography[END_REF], Section 19.3.1.1), source localization or propagation speed estimation, following e.g. the GPR methodology from [START_REF] Raissi | Machine learning of linear differential equations using Gaussian processes[END_REF] or [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF], Section 4.2.

This model is derived by putting GP priors over the initial conditions of the wave equation and in Proposition 4.1, we obtain "explicit" formulas for the covariance function of the solution process, in the form of convolutions. From Propositions 3.5 and 4.1, we obtain that the sample paths of the corresponding (nonstationary) GP all verify the wave equation in the distributional sense. When the covariance functions of the initial conditions are not smooth enough, the result from [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF] cannot be applied. Explicitly, for this PDE, choosing the commonly used 3/2-Matérn covariance functions for the initial position is enough to land outside the scope of the result from [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF] (Section 4.2.1).

We emphasize that the covariance expressions exposed in Proposition 4.1 are original and interesting in themselves, as they can be used for efficient GPR for the wave equation. Specifically, the key difference with the wave equation covariance functions presented in [START_REF] Carrizo-Vergara | A general framework for SPDE-based stationary random fields[END_REF] is that here, no stationarity assumptions are made on the solution stochastic process U . In particular the spectral measure provided by Bochner's theorem [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], which is the key tool used in [START_REF] Carrizo-Vergara | A general framework for SPDE-based stationary random fields[END_REF], is not available anymore. We thus resort to more standard integration techniques to prove Proposition 4.1.

The paper is organized as follow. For self-containment, Section 2 is dedicated to reminders on random fields and generalized functions. This Section and all the proofs are detailed enough so that this article is accessible to the analyst, the probability theorist and the statistician. In Section 3, we state and prove our new necessary and sufficient condition on random fields that are subject to linear distributional differential constraints. Section 4 is dedicated to the study of a GP model for the wave equation. We conclude in Section 5.

Background 2.1 Random fields

Let pΩ, A, Pq be a probability space. For convenience, we will assume that it is complete, i.e. that A contains the subsets of sets A P A such that PpAq " 0.

2.1.1 Basic definitions. Let D Ă R d be an open set. In this article, a random field U " pU pxqq xPD is a collection of real random variables defined on Ω. We define its sample path at point ω P Ω to be the deterministic function x Þ ÝÑ U pxqpωq, and we denote it by U ω . Given an operator acting on the sample paths of U , an event of the form tLpU q P Au will always be understood sample path wise: that is, by definition, tLpU q P Au :" tω P Ω : LpU ω q P Au. Such sets are not automatically measurable; still, they are measurable as soon as they contain an event of probability 1 (as the ones in Propositions 3.2 and 3.5), since pΩ, A, Pq is a complete probability space.

2.1.2 Measurable random fields. In view of our main theorem, a necessary notion is that of the measurability of the random field U . U is said to be measurable ( [START_REF] Doob | Stochastic processes[END_REF], p. 60 or [START_REF] Gall | Mouvement brownien, martingales et calcul stochastique[END_REF], p. 34) if it is measurable seen as a bivariate map

U : pΩ ˆD, A b BpDqq ÝÑ pR, BpRqq, pω, xq Þ Ñ U pxqpωq.
Here, BpSq denotes the Borel σ-algebra of S and A b BpDq denotes the product σ-algebra of A and BpDq.

To work with measurable random fields, one will often consider random fields U which are continuous in probability, i.e. for all x P D and ε ą 0, Pp|U pxq ´U px `hq| ą εq Ñ 0 when h Ñ 0. Indeed, continuity in probability implies the existence of a measurable modification of U , i.e. a measurable random field Ũ such that Pp Ũ pxq " U pxqq " 1 for all x P D ([13], Theorem 2.6, p. 61). One then implicitly works with Ũ . In this article, we will directly assume that we deal with measurable random fields instead of assuming any continuity regularity on the sample paths of the said stochastic process. This is because pointwise continuity is not really relevant when working with PDEs in a weak sense; actually, one of the main points of working with weakened formulations is to avoid strong (i.e. pointwise) formulations. Note however that ensuring measurability outside of the above mentioned theorem, though possible, rapidly becomes tedious (see e.g. [START_REF] Doob | Stochastic processes depending on a continuous parameter[END_REF], Theorem 2.3). A famous theorem from Kolmogorov ([10], Theorem 3.3 p. 73 and Theorem 3.4 p. 74) provides sufficient conditions for almost sure continuity of the sample paths, which in turn implies continuity in probability of the random field. This condition is phrased in terms of a sufficient Hölder control of the expectation of the increments of the process. Refinements in the case of Gaussian processes exist: see e.g. [START_REF] Adler | Random Fields and Geometry[END_REF], Theorem 1.4.1, p. 20. On a final note, the measurability assumption is discussed in [START_REF] Steinwart | Convergence types and rates in generic Karhunen-Loeve expansions with applications to sample path properties[END_REF] (Theorem 3.3), where it is shown to be a necessary condition for the existence of Karhunen-Loève expansions of second order random fields.

2.1.3 Second order random fields, Gaussian processes.. Note L 2 pPq the Hilbert space of real valued random variables X such that ErX 2 s ă `8. A stochastic process pU pxqq xPD is said to be second order if for all x P D, U pxq P L 2 pPq. One can then define its mean and covariance functions by mpxq " ErU pxqs and kpx, x 1 q " ErpU pxq ´mpxqqpU px 1 q ´mpx 1 qqs respectively. One can then also define its standard deviation function

σ : x Þ Ñ a kpx, xq. (6)
A Gaussian process pU pxqq xPD over D is a random field over D such that for any px 1 , ..., x n q P D n and any pa 1 , ..., a n q P R n , ř i a i U px i q is a Gaussian random variable; that is, the law of pU px 1 q, ..., U px n qq T is a multivariate normal distribution. The law of a GP is characterized by its mean and covariance functions ( [START_REF] Janson | Gaussian Hilbert spaces[END_REF], Section 8). We write pU pxqq xPD " GP pm, kq. Given a GP pU pxqq xPD , we will sometimes use the space LpU q " SpanpU pxq, x P Dq, i.e. the Hilbert subspace of L 2 pPq induced by U . Since L 2 pPq-limits of Gaussian random variables drawn from the same GP remain Gaussian ( [START_REF] Janson | Gaussian Hilbert spaces[END_REF], Section 1.3), LpU q only encompasses Gaussian random variables. Whereas m can be any function, the covariance function k has to be symmetric and positive definite: for all px 1 , ..., x n q in D n , the matrix pkpx i , x j qq 1ďi,jďn is symmetric and nonnegative definite.

Symmetric positive definite functions verify the Cauchy-Schwarz inequality [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] :

@x, x 1 P D, |kpx, x 1 q| ď a kpx, xq a kpx 1 , x 1 q. ( 7
)
Note that there is a one-to-one correspondence between positive definite functions and the laws of centered GPs ( [START_REF] Doob | Stochastic processes[END_REF], Theorem 3.1). We provide below two examples of radial Matérn covariance functions ( [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], pp. 84-85), which will be useful in Section 4. Set r " ||x ´x1 ||, the Euclidean distance between x and x 1 , then the following two functions are valid covariance functions, given any l ą 0:

k 1{2 px, x 1 q " expp´r{lq, k 3{2 px, x 1 q " p1 `r{lq expp´r{lq. ( 8
)
These covariance functions are widely used in machine learning, especially k 3{2 . Almost surely, the sample paths of a GP with a Matérn covariance function k ν with ν " m `1{2, m P N, are of differentiability class C m and not C m`1 . They are thus commonly used to model functions with finite smoothness.

Tools from functional analysis

We refer to [START_REF] Rudin | Functional analysis[END_REF] and [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF] for further details on generalized functions and Radon measures. In this whole subsection, D is an open set of R d .

2.2.1

Class C m functions, test functions, locally integrable functions. Given m P N, C m pDq denotes the space of real-valued functions defined over D of class C m , and C m c pDq denotes the subspace of C m pDq of functions φ whose support Supppφq is compact. Recall that Supppφq is the closure of the set tx : φpxq ‰ 0u. The space C 8 c pDq, which we will rather denote DpDq, is the space of compactly supported infinitely differentiable functions supported on D, also known as test functions. L 1 loc pDq denotes the space of measurable scalar functions f defined on D that are locally integrable, i.e. such that ş K |f | ă `8 for all compact sets K Ă D. Two locally integrable functions are equal in L 1 loc pDq when they are equal almost everywhere (a.e.) in the sense of the Lebesgue measure over R d . L 1 loc pDq is a very large space which contains the space of piecewise continuous functions, but also all the local Lebesgue spaces L p loc pDq, p ě 1 and thus all the Sobolev spaces of nonnegative exponent. It is in fact the largest space of functions that can be alternatively viewed as continuous linear forms over DpDq (see Section 2.2.4 below).

Generalized functions.

We endow DpDq with its usual LF-space topology, defined for example in [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], Chapter 13. LF stands for "strict inductive limit of Fréchet spaces". As it will appear in several places later on, we briefly describe the LF topology following [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], although this is not necessary for understanding the article. Assume that a vector space E can be written as E "

Ť n E n where pE n q is an increasing sequence of Fréchet spaces (i.e. metrizable complete locally convex topological vector spaces), such that the natural injection E n Ñ E n`1 is a linear homeomorphism over its range. The LF topology over E is defined as follow: a convex set V Ă E is a neighborhood of 0 if and only if V X E n is a neighborhood of 0 for all n. It is remarkable that LF topologies are not metrizable except if for some n 0 , E n " E n0 for all n ě n 0 ([54], Remark 13.1). In return, this allows for some other very nice topological properties to hold, e.g., LF spaces are complete ( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], Theorem 13.1).

For DpDq, the LF topology is the one corresponding to the decomposition DpDq " Ť i D Ki pDq, where D Ki pDq :" tφ P C 8 pDq : Supppφq Ă K i u, and pK i q iPN is an increasing sequence of compact subsets of D such that [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], pp. 131-133). This LF topology does not depend on the choice of pK i q iPN . An example of metric inducing the Fréchet topology of D Ki pDq is the following:

Ť i K i " D ([
d i pφ, ψq :" sup N PN 2 ´N p N,i pφ ´ψq 1 `pN,i pφ ´ψq , p N,i pφq :" max |α|ďN sup xPKi |B α φpxq|. (9)
It is given in [START_REF] Rudin | Functional analysis[END_REF], Section 1.46 p. 34 and Remark 1.38pcq p. 29. Explicitly, a sequence pφ n q Ă DpDq converges to φ P DpDq if there exists a compact set K Ă D such that Supppφ n q Ă K for all n P N and for all α P N d , ||B α φ n ´Bα φ|| 8 Ñ 0 ( [START_REF] Rudin | Functional analysis[END_REF], Theorem 6.5pf q and the remark following p. 154).

We call generalized function any continuous linear form on DpDq, i.e. any element of DpDq 1 , the topological dual of DpDq. We will rather denote it by D 1 pDq as in [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], Notation 21.1. The topology of DpDq is such that T P D 1 pDq if and only if for all compact set K Ă D, there exists C K ą 0 and a nonnegative integer n K such that

@φ P DpDq such that Supppφq Ă K, |T pφq| ď C K ÿ |α|ďn K ||B α φ|| 8 . (10)
We recall that we use the following notations: for a multi-index α " pα 1 , ..., α d q P N d , we denote |α| " α 1 `... `αd and B α :" pB x1 q α1 ...pB x d q α d where B αi xi is the α th i derivative with reference to the i th coordinate x i . Generalized functions are also called "distributions", a terminology we will only use when there is no risk of confusion with probability distributions. The duality bracket will be denoted x, y: for all φ P DpDq and T P D 1 pDq, we have xT, φy :" T pφq.

Generalized functions and differentiation.

Any generalized function T can be infinitely differentiated ( [START_REF] Rudin | Functional analysis[END_REF], Section 6.12, p. 158 or [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], pp. 248-250) according to the following definition

B α T : φ Þ ÝÑ xT, p´1q |α| B α φy. (11)
The derivative B α T is then also a continuous linear form over DpDq, i.e. B α T P D 1 pDq.

2.2.4

Regular generalized functions. Any function f P L 1 loc pDq can be injectively identified to a generalized function T f ([54], p. 224 or [START_REF] Rudin | Functional analysis[END_REF], Section 6.11, p. 157) defined as follow

@φ P DpDq, xT f , φy :" ż D f pxqφpxqdx. ( 12
)
The map L 1 loc pDq Q f Þ ÝÑ T f is linear and injective; any generalized function T that is of the form T f for some f P 1 loc pDq is said to be regular. Throughout this article, we will use the abusive notation xT f , φy " xf, φy, as if x, y were the L 2 inner product. Observe that equations ( 11) and ( 12) combined provide a flexible definition of the derivatives of any function f P L 1 loc pDq up to any order. One also sees that weak derivatives, as encountered in Sobolev spaces ( [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF], Section 9.1) and weak formulations of PDEs, are particular cases of distributional derivatives: given α P N d and two locally integrable functions f and f α , f admits f α for its α th weak derivative if and only if B α T f " T fα . One then conveniently writes B α f " f α (B α f is unique in L 1 loc pDq from the injectivity of the mapping ( 12)).

Radon measures.

This subsection and the ones that follow are only necessary for dealing with the wave equation in Section 4. In this article, we call positive Radon measure any positive measure over D that is Borel regular ( [START_REF] Evans | Measure theory and fine properties of functions[END_REF], Definition 1.9) and that has finite mass over any compact subset of D. A Radon measure is a linear combination of positive Radon measures. In [START_REF] Lang | Real and functional analysis[END_REF], Chapter 9, it is proved that the space of Radon measures over D is isomorphic to the space of continuous linear forms over C c pDq, the space of compactly supported continuous functions on D endowed with its usual LF-space topology described e.g. in [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], pp. 131-133. The corresponding isomorphism is given by

µ Þ ÝÑ T µ : # C c pDq ÝÑ R f Þ ÝÑ ş D f pxqµpdxq. (13)
We have the following facts. piq Any signed measure that admits a density f with reference to the Lebesgue measure such that f P L 1 loc pDq is a Radon measure ( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], p. 217). piiq Any Radon measure can be injectively identified to a generalized function T µ by replacing C c pDq by DpDq in equation [START_REF] Doob | Stochastic processes[END_REF]. In particular, Radon measures can be differentiated up to any order through equation [START_REF] Dieudonné | Sur les espaces de Montel métrisables[END_REF]. piiiq Any Radon measure µ, can be uniquely written as µ " µ `´µ ´where µ ànd µ ´are positive Radon measures ( [START_REF] Lang | Real and functional analysis[END_REF], Chapter 9). We then define its total variation by |µ| :" µ ``µ ´.

Finite order generalized functions.

A generalized function T is said to be of finite order if there exists a nonnegative integer m such that one can take n K " m, independently of K, in the definition of the continuity of T , i.e. equation [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF]. The order of T is then the smallest of those integers m. The space of generalized functions of order m is isomorphic to C m c pDq 1 , the space of continuous linear forms over C m c pDq, when C m c pDq is endowed with its usual LF-space topology ( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], pp. 131-133). The key property for us is that such generalized functions can be represented thanks to Radon measures. If L is of order m, there exists a family of Radon measures tµ p u |p|ďm over D such that

T " ÿ |p|ďm B p µ p , (14)
where the equality in equation ( 14) holds in D 1 pDq and C m c pDq 1 ([54], p 259). Among the finite order generalized functions are those that are compactly supported, i.e. those for which the measures µ p such that T " ř |p|ďm B p µ p all have compact support.

Convolution with generalized functions.

As above, we consider C m c pR d q endowed with its LF-space topology. Let f P C m c pR d q and T P C m c pR d q 1 . Note τ x f the function y Þ ÝÑ f py ´xq and f the function y Þ ÝÑ f p´yq. Then ( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], p. 287, Section 27) one may define the convolution between T and f by

T ˚f : x Þ ÝÑ xT, τ ´x f y, (15)
and T ˚f is a function in the classical sense, i.e. defined pointwise. When T is a regular generalized function, equation [START_REF] Estrade | Anisotropic gaussian wave models[END_REF] reduces to the usual convolution of functions through the identification defined in equation [START_REF] Doob | Stochastic processes depending on a continuous parameter[END_REF]. Similarly if T is in fact a Radon measure µ:

pT ˚f qpxq " ż R d f px ´yqµpdyq. ( 16
)
More general definitions of generalized function convolution are available ( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], Chapter 27) but this one is sufficient for our use.

2.2.8

Tensor product of generalized functions. For two generalized functions T 1 P D 1 pD 1 q and T 2 P D 1 pD 2 q, T 1 b T 2 P D 1 pD 1 ˆD2 q denotes their tensor product( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], pp. 416-417), which is uniquely determined by the following tensor property:

@φ 1 P DpD 1 q, @φ 2 P DpD 2 q, xT 1 b T 2 , φ 1 b φ 2 y " xT 1 , φ 1 y ˆxT 2 , φ 2 y. (17)
T 1 b T 2 reduces to the tensor product of functions (respectively, measures) when T 1 and T 2 are functions (respectively, measures) through the identification of equation ( 12) (respectively, equation ( 13)).

Random fields under linear differential constraints

The results in this section state that under suitable assumptions over the first two moments of a given second order random field U " pU pxqq xPD , sample path degeneracy properties with reference to differential constraints can be read on the first two moments of U , namely the mean function and the functions k x : y Þ ÝÑ kpx, yq, where k is the covariance function of U . This is remarkable because the space induced by the sample paths of U is a priori much larger than the space spanned by the functions k x , x P D. Moreover, the functions k x are "accessible", i.e. checking that these functions indeed verify the linear constraint can usually be done with direct computations.

We begin by recalling a result from [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF] in the case of pointwise defined derivatives. We next state and prove a result similar to that of [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF], where we interpret the derivatives in the distributional sense.

The case of classical derivatives

We start by properly defining the notion of strong solutions of a PDE. Definition 3.1 (Strong/classical solutions). Let L be a differential operator defined as in equation [START_REF] Adler | The Geometry of Random Fields[END_REF], with continuous coefficients. We say that a function u is a classical or strong solution to the PDE Lpuq " 0 if u is n times differentiable and u verifies the PDE pointwise:

@x P D, Lpuqpxq " ÿ |α|ďn a α pxqB α upxq " 0. ( 18
)
Note that the space of n times differentiable functions does not have the nice topological properties of C n pDq and in most cases met in practice, one rather requires that strong solutions lie in C n pDq. It is however in the sense of the definition 3.1 that the theorem from [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF] is best understood. This theorem, which we remind in Proposition 3.2, is the one proved and used in [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF] to build a Gaussian process whose sample paths are all strong solutions to the Laplace equation on a 2D circular domain.

We first introduce some notations. Let pU pxqq xPD be a centered Gaussian process with covariance function k. Denote FpD, Rq the space of real-valued pointwise-defined functions on D (often alternatively denoted R D ). We will only use FpD, Rq as a set, therefore we do not consider any topology over it. We refer to [START_REF] Schwartz | Sous-espaces Hilbertiens d'espaces vectoriels topologiques et noyaux associés (noyaux reproduisants)[END_REF], Section 9, for details on FpD, Rq seen as a topological vector space. Denote H k the reproducing kernel Hilbert space (RKHS) associated to k (see [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF], Definition 1 p. 7 and Theorem 3, p. [START_REF] Fasshauer | Meshfree approximation methods with MATLAB[END_REF]). H k is a Hilbert space of pointwise-defined functions (i.e. H k Ă FpD, Rq as sets), such that the pointwise evaluation maps l x : f Þ Ñ f pxq are continuous functionals. Although belonging to FpD, Rq does not seem very restrictive at first glance, this clashes with the usual L p and Sobolev spaces encountered in PDE theory, which are sets of functions defined up to a set of null Lebesgue measure. Proposition 3.2 (sample paths of GPs under linear constraints [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF]). Let `U pxq ˘xPD " GP p0, kq. Note for all x P D the function k x : y Þ ÝÑ kpx, yq. Let E be a real vector space of functions defined on D that contains the sample paths of U almost surely and L : E ÝÑ FpD, Rq be a linear operator. Assume that for all x P D, LpU qpxq P LpU q, where LpU q is the closure of SpantU pxq, x P Du in L 2 pPq. Then there exists a unique linear operator L : H k ÝÑ FpD, Rq such that @x, y P D, ErLpU qpxqU pyqs " L pk y qpxq, and such that for all x P D, h P H k and sequence ph n q Ă H k such that h n Ý Ñ h for the topology of H k , we have L ph n qpxq ÝÑ L phqpxq. Finally, the following statements are equivalent: (i) PpLpU q " 0q " 1.

(ii) @x P D, L pk x q " 0.

A sufficient condition ensuring that the sample paths of a GP lie in C n pDq is found in [START_REF] Adler | Random Fields and Geometry[END_REF], Theorem 1.4.2. More broadly, both necessary and sufficient conditions over the first two moments of a GP for its sample path to be (Hölder) continuous are well-known: see e.g. [START_REF] Adler | The Geometry of Random Fields[END_REF], Theorems 3.3.3 and 8.3.2.

The proof of Proposition 3.2 heavily relies on the Loève isometry ([6], Theorem 35, p. 65) between the two Hilbert spaces H k and LpU q (see Section 2.1.3 for details on LpU q). This theorem can be applied when L is a differential operator as discussed in [START_REF] Ginsbourger | On degeneracy and invariances of random fields paths with applications in Gaussian process modelling[END_REF]. However, in Proposition 3.2, the differential operator L of order n has to be valued in FpD, Rq; in particular for u P E, the function Lpuq has to be defined pointwise in order to use the Loève isometry. To summarize, in all generality the derivatives in L have to be understood in a classical sense and E has to be contained in D n pDq, the space of n times differentiable functions on D. Requiring that E Ă D n pDq is a very strong assumption with reference to the sample paths of U ; furthermore, this is not compliant with the usual way of studying PDEs where derivatives are understood in a weaker sense. We present in Proposition 3.5 an adaptation of Proposition 3.2 where the derivatives are understood in the distributional sense. By transferring all the derivatives on the test function, we will be liberated from any differentiability assumptions over the sample paths of U , effectively replacing D n pDq with L 1 loc pDq. Finally, the random field U will not be assumed Gaussian and will only be required to be measurable second order.

The case of distributional derivatives

3.2.1 Distributional solutions of PDEs. In this section, we elaborate a bit more on the notion of distributional solutions to a given PDE. Let L " ř |α|ďn a α pxqB α be a linear differential operator, and assume for the moment that its coefficients are infinitely differentiable. We briefly recall the steps described in the introduction that lead to the definition of distributional solutions presented in equation ( 4). Start from a strong solution u of class C n of Lpuq " 0, multiply this PDE by a test function φ P DpDq, integrate over D and perform |α| integration by parts to transfer all derivatives from u to φ. Since the support of φ is a compact subset of the open set D, the boundary terms of each integration by parts vanish, leading to

@φ P DpDq, ż D upxq ÿ |α|ďn p´1q |α| B α pa α φqpxqdx " 0. ( 19
)
Following equation [START_REF] Fasshauer | Meshfree approximation methods with MATLAB[END_REF] we introduce L ˚, the formal adjoint of L, acting on DpDq, defined by the following formula ( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], pp. 247-249)

L ˚: φ Þ ÝÑ ÿ |α|ďn p´1q |α| B α pa α φq. (20)
Note that for equation [START_REF] Fasshauer | Meshfree approximation methods with MATLAB[END_REF] to be well defined, the assumptions that u P L 1 loc pDq and a α P C |α| pDq are sufficient. More precisely, these assumptions are enough to show that the map Lpuq defined by duality Lpuq :

# DpDq ÝÑ R φ Þ ÝÑ ş D L ˚pφqpxqupxqdx (21)
defines a continuous linear form over DpDq, i.e. Lpuq P D 1 pDq (see equations ( 10) and ( 26) for a rigorous proof of this statement). This definition extends the definition of distributional derivatives from Section 2.2.3 to differential operators. By construction, L and L ˚verify a duality identity: given φ P DpDq and u P L 1 loc pDq, xLpuq, φy " xu, L ˚pφqy. As in Section 2.2.4, the assumption that u P L 1 loc pDq is in fact a continuity assumption over the associated linear form Lpuq (a more general and theoretical analysis of such observations can be found in [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], pp. 247-251). This finally leads to the following definition, following e.g. [START_REF] Duistermaat | Distributions[END_REF], p. 10:

Definition 3.

(Distributional solutions). A function u P L 1

loc pDq is said to be a solution to the PDE Lpuq " 0 in the sense of distributions if Lpuq " 0 in D 1 pDq, i.e. when Lpuq is seen as en element of D 1 pDq through equation ( 21) and 0 is the null linear form over DpDq.

As weak derivatives are a particular case of distributional derivatives (Section 2.2.4), one expects that the distributional solutions of a PDE that admit some weak derivatives are in fact weak solutions, i.e. solutions of some weak formulation of that PDE. Rigorous statements of this general fact have to be checked on a case-by-case basis, depending on the weak formulation at hand (a more in-depth discussion falls outside of the scope of this article). As an example, this is the case for the weak formulation of elliptic PDEs in H 1 0 pDq (see e.g. [START_REF] Evans | Partial differential equations[END_REF], Section 6.2), where H 1 0 pDq is the closure of DpDq in the Sobolev space H 1 pDq :" tu P L 2 pDq : ∇u exists as a weak derivative and ∇u P L 2 pDq d u.

Remark 3.4 (Measure-valued solutions of PDEs). Although it is not the main focus of the paper, we can even allow u in Definition 3.3 to be a Radon measure by replacing upxqdx with µpdxq in equation ( 21). This will be useful in Section 4.1, where we will encounter a measure-valued PDE solution which is central from a physical viewpoint, with the wave equation's Green's function (it is not actually a function!). Notice that weak formulations in Sobolev spaces, say H 1 pDq, are not well-equipped to work with such solutions, and our distributional framework becomes needed.

Random fields under distributional differential constraints.

We can now state the following proposition, based on Definition 3.3.

Proposition 3.5 (sample paths of random fields under linear differential constraints, distributional derivatives). Let D Ă R d be an open set and let L "

ř a α pxqB α , |α| ď n, be a linear differential operator of order n with coefficients a α pxq P C |α| pDq. Let U " `U pxq ˘xPD be a measurable second order random field with mean function mpxq and covariance function kpx, x 1 q. For all x P D, note k x : y Þ ÝÑ kpx, yq. Suppose that m P L 1 loc pDq and σ P L 1 loc pDq, where σ : x Þ Ñ kpx, xq 1{2 . 1) Then PpU P L 1 loc pDqq " 1 and for all x P D, k x P L 1 loc pDq. 2) Suppose that Lpmq " 0 in the sense of distributions. Then the following statements are equivalent:

(i) PpLpU q " 0 in the sense of distributionsq " 1.

(ii) @x P D, Lpk x q " 0 in the sense of distributions.

Explicitly, by piq we mean that there exists a set A P A with PpAq " 1 such that for all ω P A,

@φ P DpDq, xU ω , L ˚φy " ż D U ω pxqL ˚φpxqdx " 0. ( 22
)
The fact that the functions x Þ ÝÑ U ω pxq and y Þ ÝÑ k x pyq lie in L 1 loc pDq ensure the existence of the integrals in equations [START_REF] Gapaillard | Sur les processus linéaires définis sur un espace nucléaire[END_REF] (see Point 2 of the proof of Proposition 3.5) as well as the continuity of the associated linear forms over DpDq, following the definition of equation [START_REF] Fuselier | Refined error estimates for matrix-valued radial basis functions[END_REF]. The assumption that a α P C |α| pDq is not very strong, in the sense that it is the minimal assumption to ensure that the adjoint L ˚is well-defined (equation ( 20)), and thus that Definition 3.3 even makes sense. Likewise, requiring that σ P L 1 loc pDq is not very restrictive (see Section 2.2.1). However, ensuring the measurability of the random process U is more demanding in practice, because it is difficult to ensure this property outside of having continuity in probability (see Section 2.1.2).

The following lemma will be crucial for the proof of Proposition 3.5:

Lemma 3.6. DpDq is sequentially separable, i.e. there exists a countable subset F Ă DpDq such that for all φ P DpDq, there exists a sequence pφ n q Ă F such that φ n Ñ φ in DpDq for its LF topology.

Recall that a topological space E is separable if there exists a countable subset F Ă E such that its closure in E is equal to E. If the topology of E is metrizable (as e.g. for Fréchet spaces), sequential separability and separability are equivalent. If this topology is not metrizable (as e.g. for LF spaces), then sequential separability implies separability but the converse need not hold. Below, we provide a short proof of Lemma 3.6, as we could not find it in the literature. The weaker property that DpDq is separable is already difficult to track down, see e.g. [START_REF] Gapaillard | Sur les processus linéaires définis sur un espace nucléaire[END_REF], Corollaire (1).2, p. 78 or [START_REF] Gel | Applications of harmonic analysis[END_REF], p. 73, (3).

Proof. We first show that the spaces D Ki pDq introduced in Section 2.2.2 are separable Fréchet spaces. The Fréchet topology of D Ki pDq is the one induced by the usual Fréchet topology of C 8 pDq when D Ki pDq is seen as a subspace of C 8 pDq ( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], pp. 131-132). As a Fréchet space, C 8 pDq is metrizable ( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], p. 85). But C 8 pDq is also a Montel space ( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], Proposition 34.4, p. 357): as a metrizable Montel space, it is automatically separable ( [START_REF] Schaefer | of Graduate Texts in Mathematics[END_REF], p. 195 or [START_REF] Dieudonné | Sur les espaces de Montel métrisables[END_REF]). Thus D Ki pDq is also separable as a subset of the separable metrizable space C 8 pDq ([7], Proposition 3.25, p. 73).

Denote now F i a countable dense subset of D Ki pDq and consider F :" Ť iPN F i . Let φ P DpDq and i P N such that Supppφq Ă K i , where pK i q iPN is the sequence of compact sets from Section 2.2.2. Then φ P D Ki pDq and there exists a sequence pφ n q Ă F i Ă F such that φ n Ñ φ in the sense of the Fréchet topology of D Ki pDq, i.e. the metric d i in equation [START_REF] Chen | Solving and learning nonlinear PDEs with Gaussian processes[END_REF]. From equation ( 9), ||B α φ n ´Bα φ|| 8 Ñ 0 for all α P N d . Since Supppφ n q Ă K i for all n P N, we have that φ n Ñ φ in DpDq (see Section 2.2.2).

We are now able to prove Proposition 3.5.

Proof. Suppose first that U is centered, i.e. m " 0. 1) We begin by showing that the sample paths of U almost surely lie in L 1 loc pDq. Note first that thanks to the Cauchy-Schwarz inequality, Er|U pxq|s ď σpxq. Now, let pK i q iPN be an increasing sequence of compact subsets of D such that Ť iPN K i " D. Using Tonelli's theorem, we have that for any n P N,

E " ż Ki |U pxq|dx ȷ " ż Ki Er|U pxq|sdx ď ż Ki σpxqdx ă `8, (23)
since σ P L 1 loc pDq. Note that in order for the integrals above to be well defined, imposing that U is a measurable random field cannot be circumvented. Equation ( 23) yields a set B n Ă Ω of probability 1 over which the random variable ω Þ ÝÑ ş Ki |U ω pxq|dx is finite (from Fubini's theorem again, the map ω Þ ÝÑ ş Ki |U ω pxq|dx is measurable). Consider now the set B " Ş nPN B n which remains of probability 1. For all compact subset K Ă D, there exists an integer n K such that K Ă K n K and thus for all ω P B,

ż K |U ω pxq|dx ď ż Kn K |U ω pxq|dx ă `8, (24)
which shows that the sample paths of U lie in L 1 loc pDq almost surely. Similarly, we check that for all x P D, k x lies in L 1 loc pDq: for any compact set K, since σ P L 1 loc pDq and because of equation ( 7),

ż K |k x pyq|dy " ż K |kpx, yq|dy ď σpxq ż K σpyqdy ă 8.
2) Let us check in advance that whatever f P L 1 loc pDq, the map T pf q : φ Þ ÝÑ xf, L ˚φy is a continuous linear form over DpDq. Since a α P C |α| pDq, we can apply Leibniz' rule on L ˚φ " ř |α|ďn p´1q |α| B α pa α φq. This yields a family tf α u |α|ďn of continuous functions over D such that @φ P DpDq, @x P D, L ˚φpxq "

ÿ |α|ďn f α pxqB α φpxq. (25)
For all f P L 1 loc pDq, for all compact set K Ă D and for all φ P DpDq such that Supppφq Ă K, we have SupppL ˚φq Ă K and equation ( 25

This proves that T pf q : φ Þ ÝÑ xf, L ˚φy is a continuous linear form over DpDq (see equation ( 10)). piq ùñ piiq: Suppose piq. Let φ P DpDq. There exists a set A Ă Ω such that PpAq " 1 and such that

@ω P A, xU ω , L ˚φy " ż D U ω pxqL ˚ϕpxqdx " 0.
Multiplying equation above with the random variable U px 1 q, taking the expectation and formally permuting (for now) the integral and the expectation, we obtain

0 " E " U px 1 q ż D U pxqL ˚φpxqdx ȷ " ż D L ˚φpxqErU pxqU px 1 qsdx " ż D L ˚φpxqkpx, x 1 qdx " xk x 1 , L ˚φy.
The integral-expectation permutation is justified by writing down the expectation as an integral and using Fubini's theorem, checking that the below quantity is finite. We use Tonelli's theorem and the Cauchy-Schwarz inequality:

E " ż D |U px 1 qU pxqL ˚φpxq|dx ȷ " ż D |L ˚φpxq|Er|U pxqU px 1 q|sdx ď ż D |L ˚φpxq|ErU pxq 2 s 1{2 ErU px 1 q 2 s 1{2 dx ď σpx 1 q ż D |L ˚φpxq|σpxqdx ă `8.
Indeed, since σ P L 1 loc pDq, setting f " σ in equation [START_REF] Graepel | Solving noisy linear operator equations by gaussian processes: Application to ordinary and partial differential equations[END_REF] shows that the last integral is indeed finite. Thus, @x P D, @φ P DpDq, xk x , L ˚φy " 0 which proves that piq ùñ piiq. piiq ùñ piq: Suppose piiq. Let φ P DpDq, we have xk x 1 , L ˚φy " 0. Multiplying this with L ˚φpx 1 q and integrating with reference to x 1 yields 0 "

ż D L ˚φpx 1 q ż D L ˚φpxqkpx, x 1 qdxdx 1 " ż D ż D L ˚φpxqL ˚φpx 1 qErU pxqU px 1 qsdxdx 1 .
Permuting formally the expectation and the integrals (justified in equation ( 27)) yields 0 "

ż D ż D L ˚φpxqL ˚φpx 1 qErU pxqU px 1 qsdxdx 1 " E «˜ż D L ˚φpxqU pxqdx ¯2ff " E " xU, L ˚φy 2 ‰ ,
and thus xU, L ˚φy " 0 a.s. : there exists A φ P A with PpA φ q " 1 such that @ω P A φ , xU ω , L ˚φy " 0. We justify the expectation-integral permutation with the computation below

ż D ż D |L ˚φpxqL ˚φpx 1 q|Er|U pxqU px 1 q|sdxdx 1 ď ż D ż D |L ˚φpxqL ˚φpx 1 q|σpxqσpx 1 qdxdx 1 ď ˜żD |L ˚φpxq|σpxqdx ¸2 ă `8. (27)
As previously, setting f " σ in equation [START_REF] Graepel | Solving noisy linear operator equations by gaussian processes: Application to ordinary and partial differential equations[END_REF] shows that the integral above is indeed finite. This does not finish the proof as we need to find a set A with PpAq " 1, independently from φ, such that @ω P A, xU ω , L ˚φy " 0. For this we shall use Lemma 3.6. Let

A :" B X `č φPF A φ ˘, (28)
where the set F is introduced in Lemma 3.6. Then PpAq " 1 since PpBq " 1, PpA φ q " 1 and F is countable. Let ω P A. Since U ω P L 1 loc pDq, equation [START_REF] Graepel | Solving noisy linear operator equations by gaussian processes: Application to ordinary and partial differential equations[END_REF] shows that the map T ω : φ Þ ÝÑ xU ω , L ˚φy is a continuous linear form on DpDq. In particular, Theorem 6.6pcq p. 155 from [START_REF] Rudin | Functional analysis[END_REF] states that T ω is in fact sequentially continuous. Let φ P DpDq and pφ n q Ă F be such that φ n Ñ φ in DpDq, from Lemma 3.6. From the sequential continuity of T ω , T ω pφq " lim nÑ8 T ω pφ n q " 0 since @n P N, T ω pφ n q " 0. That is, we have proved that @ω P A, @φ P DpDq, xU ω , L ˚φy " T ω pφq " 0.

Since PpAq " 1, this shows that piiq ùñ piq.

When U is not centered, consider the centered random field V defined by V pxq " U pxq´mpxq for which the above proof can be applied. Since L is linear and m is assumed to verify Lpmq " 0 in the sense of distributions, the probabilistic sets A U " tLpU q " 0 in the sense of distributionsu and A V " tLpV q " 0 in the sense of distributionsu coincide and thus, A Ă A U . Finally, U and V have the same covariance function kpx, x 1 q. Thus, PpLpU q " 0 in the distrib. senseq " 1 ðñ PpLpV q " 0 in the distrib. senseq " 1 ðñ @x P D, Lpk x q " 0 in the distrib. sense, which finishes the proof in the general case.

Remark 3.7. Distributional solutions are the weakest types of solutions for PDEs. In general, additional regularity conditions have to be imposed to obtain physically realistic solutions, such as Sobolev regularity or entropy conditions as for nonlinear hyperbolic PDEs [START_REF] Serre | Systems of conservation laws[END_REF]. However, every step in the above proof remains valid when replacing φ P DpDq with φ P C n c pDq. Although we have not explicited the usual topology of C n c pDq in this article, we state that this is enough to show that the equalities stated in Proposition 3.5 also hold in C n c pDq 1 , the space of finite order generalized functions of order n, rather than just in D 1 pDq. C n c pDq 1 is a smaller space than D 1 pDq, though less used in functional analysis than D 1 pDq.

We partially recover Proposition 3.2 when the sample paths of U are n times differentiable with locally integrable n th derivative and k P C n,n pD ˆDq. Indeed, in that case one can show that if L " ř |α|ďn a α pxqB α , then we simply have L " L in Proposition 3.2. Additionally, LpU ω q and Lpk x q both lie in FpD, Rq X L 1 loc pDq. In that framework, Proposition 3.2 states that

@x P D, Lpk x q " 0 ðñ PpLpU q " 0q " 1, (29)
where the function equalities of the form Lpf q " 0 in equation ( 29) are valid everywhere on D.

In contrast, for any function g that lies in L 1 loc pDq, we have g " 0 in the sense of distributions ðñ g " 0 a.e.

Equation ( 30) is just another way of saying that the linear map f Þ ÝÑ T f given in ( 12) is injective. Following equation ( 30), Proposition 3.5 states a slightly weaker result than (29), namely that @x P D, Lpk x q " 0 a.e. ðñ PpLpU q " 0 a.e.q " 1.

If we actually have that the sample paths of U lie in C n pDq, nullity almost everywhere implies nullity everywhere and we recover equation ( 29) from equation [START_REF] Kuchment | Mathematics of photoacoustic and thermoacoustic tomography[END_REF]. Instead of having the sample paths of U lie in C n pDq though, one may rather encounter the case where U is mean-square differentiable up to a certain order m. Under some continuity assumptions over the covariance function of U and up to suitable modifications, [START_REF] Scheuerer | Regularity of the sample paths of a general second order random field[END_REF] showed that the sample paths of the mean-square differentiated process are actually weak derivatives of the sample paths of U . As observed after Definition 3.3, we thus expect that the sample paths of the mean-square differentiable random fields verifying Point 2, piiq of Proposition 3.5 are solutions of some weak formulation of the PDE, rather than just distributional solutions.

Example (A first order PDE). Consider a continuous, nondifferentiable one dimensional covariance function k 0 : R ˆR Ñ R, for example k 0 px, x 1 q " expp´|x ´x1 |q. It is then readily checked that the function k : R 2 ˆR2 Ñ R defined by kppx, yq, px 1 , y 1 qq " k 0 px ´y, x 1 ´y1 q is positive definite and verifies Point 2, piiq of Proposition 3.5 for the PDE

B x u `By u " 0 in R 2 . ( 32
)
Consider now a centered second order random field pU px, yqq px,yqPR 2 with covariance function k, passing to a measurable version of U if necessary (it exists from Section 2.1.2, as the continuity of k yields the continuity in probability of U ). Then almost surely, its sample paths verify the PDE [START_REF] Lang | Real and functional analysis[END_REF] in the sense of distributions, even though they are not expected to be differentiable. An example of random field whose covariance function is k as defined above, is the GP pU 0 px ýqq px,yqPR 2 where pU 0 pxqq xPR " GP p0, k 0 q. These formulas can be obtained by viewing the PDE (32) as a transport equation under the condition that U px, 0q " U 0 pxq, following the same approach as in the upcoming Section 4.2.

A heredity property for Gaussian process regression

3.3.1 Gaussian process regression in a nutshell. GPs can be used for function interpolation. Let u be a function defined on D for which we know a dataset of values B " tupx 1 q, ..., upx n qu. Conditioning the law of a GP pU pxqq xPD " GP pm, kq on the database B yields a second GP Ũ given by Ũ pxq :" pU pxq|U px i q " upx i q, i " 1, ..., nq. The law of Ũ is known: p Ũ pxqq xPD " GP p m, kq. m and k are given by the so-called Kriging equations ( 33) and [START_REF] Lange-Hegermann | Linearly constrained Gaussian processes with boundary conditions[END_REF]. Let X " px 1 , ..., x n q T , denote mpXq the column vector such that mpXq i " mpx i q, kpX, Xq the square matrix such that kpX, Xq ij " kpx i , x j q and given x P D, kpX, xq the column vector such that kpX, xq i " kpx i , xq. Suppose that KpX, Xq is invertible, then [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] " mpxq = mpxq `kpX, xq T kpX, Xq ´1pupX q ´mpXqq, (33) kpx, x 1 q = kpx, x 1 q ´kpX, xq T kpX, Xq ´1kpX, x 1 q.

The Kriging standard deviation function is then given by σpxq " kpx, xq 1{2 .

The so-called Kriging mean m plays the role of an approximation of u; in particular, it interpolates u at the observation points: mpx i q " upx i q for all i " 1, ..., n. Moreover, the Kriging covariance k can be used to further control the distance between u and m.

Conditioned

Gaussian processes under linear differential constraints. We can now state the following corollary, which draws the consequences of Proposition 3.5 when applied to GPR.

Proposition 3.8 (Heredity of Proposition 3.5 to conditioned GPs). Let D and L be as defined in Proposition 3.5. Let pU pxqq xPD " GP pm, kq be a Gaussian process that verifies the assumptions of Proposition 3.5. Suppose also that

Lpmq " 0 and @x P D, Lpk x q " 0 both in the sense of distributions.

piq Then whatever the integer p, the vector u " pu 1 , ..., u p q T P R p and the vector X " px 1 , ..., x p q T P D p such that kpX, Xq is invertible, the Kriging mean mpxq and the Kriging standard deviation function σ both lie in L 1 loc pDq, and we have Lp mq " 0 and @x P D, Lp kx q " 0 both in the sense of distributions.

where m and k are defined in equations [START_REF] Lange-Hegermann | Algorithmic linearly constrained Gaussian processes[END_REF] and [START_REF] Lange-Hegermann | Linearly constrained Gaussian processes with boundary conditions[END_REF]. piiq As such, the sample paths of the conditioned Gaussian process `Ũ pxq ˘xPD defined by Ũ pxq " pU pxq|U px i q " u i @i " 1, ..., pq are almost surely solutions of the equation Lpf q " 0 in the sense of distributions:

PpLp Ũ q " 0 in the sense of distributionsq " 1.

Proof. Note first that for all x P D, kpx, xq ď kpx, xq ( [START_REF] Fasshauer | Meshfree approximation methods with MATLAB[END_REF], p. 117). Thus the function σ : x Þ ÝÑ kpx, xq 1{2 also lies in L 1 loc pDq. Point piq is then a direct consequence of the definition of m and k in equations ( 33) and [START_REF] Lange-Hegermann | Linearly constrained Gaussian processes with boundary conditions[END_REF], and the linearity of L. Proposition 3.5 can then be applied conjointly with piq, which yields point piiq since the mean and covariance functions of the GP Ũ are m and k (see equations [START_REF] Lange-Hegermann | Algorithmic linearly constrained Gaussian processes[END_REF] and ( 34)). Proposition 3.8 shows that when U is a GP, the results of Proposition 3.5 are inherited on the conditioned posterior process Ũ . One weak consequence of Proposition 3.8 is that if GPR is performed with a covariance function k that verifies point piiq of Proposition 3.5, then all the possible Kriging means provided by GPR remain solutions of the PDE Lp mq " 0.

Gaussian processes and the 3 dimensional wave equation

The formalism we used in the previous section is necessary to tackle hyperbolic PDEs as in some cases, their solutions only verify the PDE in a weaker sense, e.g. the distributional sense ([16], Sections 2.1.1 and 7.2). Hyperbolic PDEs are typically encountered when describing finite speed propagation phenomena and their prototype is the wave equation (see equation [START_REF] Mendes | Bayesian inference in the numerical solution of Laplace's equation[END_REF]); this equation is central in a number of fields such as acoustics, electromagnetics and quantum mechanics. In this section, we derive a GP model for the solutions of the homogeneous 3D wave equation, with explicit covariance formulas in the form of convolutions. We show on one example that the model we obtain below is capable of dealing with an initial speed v 0 that is piecewise continuous and an initial position u 0 that has piecewise continuous derivatives, when the initial discontinuity surfaces are "nice enough". This is an advantage with reference to the previous models, where the sample paths actually had to be sufficiently differentiable to obtain sample path degeneracy with reference to the PDE. lw " 0 @px, tq P R 3 ˆR˚, wpx, 0q " u 0 pxq @x P R 3 , pB t wqpx, 0q " v 0 pxq @x P R 3 .

(

) 37
Throughout this article, we will refer to u 0 as the initial position and v 0 as the initial speed.

where F t and 9 F t are known generalized functions. That is, the function wpx, tq above is a solution of the system [START_REF] Mendes | Bayesian inference in the numerical solution of Laplace's equation[END_REF], in which t is now allowed to lie in R rather than R ˚. The existence of such an extension is possible because of the time reversibility of the wave equation (in the language of semigroup theory, its semigroup can be embedded in a group, [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF], Theorem 4.5 p. 222). In dimension 3, F t and 9 F t are compactly supported generalized functions of order 0 and 1 respectively. They are given by

F t " σ c|t| 4πc 2 t and 9 F t " B t F t @t P R, (39)
Proposition 4.1. Define the two functions k wave v pz, z 1 q " rpF t b F t 1 q ˚kv spx, x 1 q, (43)

k wave u pz, z 1 q " rp 9 F t b 9 F t 1 q ˚ku spx, x 1 q. ( 44
)
(i) Then U " pU pzqq zPR 3 ˆR and V " pV pzqq zPR 3 ˆR as defined in [START_REF] Raissi | Machine learning of linear differential equations using Gaussian processes[END_REF] are two independent centered GPs with covariance functions k wave u and k wave v respectively. Consequently, pW pzqq zPR 3 ˆR is a centered GP whose covariance function is given by

k W pz, z 1 q " k wave v pz, z 1 q `kwave u pz, z 1 q. ( 45
)
(ii) Conversely, any measurable centered second order random field with covariance function k W has its sample paths solution of the wave equation [START_REF] Mendes | Bayesian inference in the numerical solution of Laplace's equation[END_REF], in the sense of distributions, almost surely.

The formulas ( 43) and ( 44) can easily be derived formally, by running computations as if F t and 9

F t were regular generalized functions (Section 2.2.4). This is somewhat justified because any generalized function can be approximated with a sequence of smooth compactly supported functions, by a "cutting and regularizing" argument ( [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], Theorem 28.2, Chapter 28). However, checking that this procedure passes to the limit everywhere is tedious. Here, we rather make use of representations of F t and 9 Proof. piq : first we prove that U and V are GPs. Since U 0 and V 0 are GPs, LpU 0 q and LpV 0 q are only comprised of Gaussian random variables (see Section 2.1.3). We then rely on the Kirschoff formula [START_REF] Owhadi | Bayesian numerical homogenization[END_REF], writing the integrals as limits of Riemann sums. We start with V , that is, we focus on the first term in Kirschoff's formula [START_REF] Owhadi | Bayesian numerical homogenization[END_REF]. To show that V is a Gaussian process, we only need to show that for any z, V pzq P LpV 0 q as this will ensure the Gaussian process property. Since the sample paths of V 0 are continuous almost surely, there exists a sequence of numbers pa n k q Ă R and points py n k q Ă Sp0, 1q such that for almost any ω P Ω,

F
V pzqpωq " pF t ˚V 0 ω qpxq " t ż Sp0,1q V 0 px ´c|t|γqpωq dΩ 4π " t 4π ż 2π 0 ż π 0 V 0 px ´c|t|γpθ, ϕqqpωq sinpθqdθdϕ " lim nÑ8 n ÿ k"1 a n k V 0 px ´c|t|y n k qpωq.
This shows that V pzq is the a.s. limit of the sequence of centered Gaussian random variables pY n q Ă LpV 0 q, where Y n " ř n k"1 a n k V 0 px´c|t|y n k q; Y n is Gaussian because V 0 is a GP. Almost sure convergence implies convergence in law. From [START_REF] Gall | Mouvement brownien, martingales et calcul stochastique[END_REF], Proposition 1.1, V pzq is normally distributed and the convergence also takes place in L 2 pPq. Therefore, V pzq P LpV 0 q and V is a Gaussian process. From the same proposition, V pzq is centered because the variables Y n are centered. Note that since F t is supported on the compact set Sp0, c|t|q, we only required the sample paths of V 0 to be continuous rather than continuous and compactly supported.

We apply the same reasoning to U , by applying the above steps to the second part of Kirschoff's formula [START_REF] Owhadi | Bayesian numerical homogenization[END_REF]. One's ability to write out the integrals as a limit of Riemann sums is ensured when the sample paths of U 0 lie in C 1 pR 3 q.

Finally, since U 0 and V 0 are independent, LpU 0 q and LpV 0 q are orthogonal in L 2 pPq. Since LpU q Ă LpU 0 q and likewise for V , U and V are independent Gaussian processes as for Gaussian random variables, independence is equivalent to null covariance. Finally, the sum of independent Gaussian random variables is a Gaussian random variable. Therefore LpW q Ă LpU q `LpV q is only comprised of Gaussian random variables and W is a Gaussian process. Now, we prove that ErU pzqU pz 1 qs " rp 9 F t b 9 F t 1 q ˚ku spx, x 1 q.

The main argument is Fubini's theorem for Radon measures. For this we use the fact that 9 F t is a distribution of order 1 and can be identified to a sum of derivatives of measures (see equation ( 14)): for all t P R, there exists tµ t α u αPN 3 ,|α|ď1 a family of Radon measures such that

9 F t " ÿ |α|ď1 B α µ t α in the sense of distributions. (47)
Moreover, 9 F t is compactly supported, therefore all the measures µ t α are also compactly supported. First, we write U ω pzq in integral form:

U ω pzq " `9 F t ˚U 0 ω ˘pxq " x 9 F t , τ ´x Ǔ 0 ω y " A ÿ |α|ď1 B α µ t α , τ ´x Ǔ 0 ω E (48)
"

ÿ |α|ď1 xµ t α , p´1q |α| B α τ ´x Ǔ 0 ω y " ÿ |α|ď1 ż R 3 p´1q |α| B α U 0 ω px ´yqµ t α pdyq. (49)
Before applying Fubini's theorem, we need to check an integrability condition. Let α P N 3 be such that |α| ď 1. Recall that |µ α t | is defined in Section 2.2.5; denote also σ B α U 0 pxq " a VarpB α U 0 pxqq. Since the sample paths of U 0 lie in C 1 pDq a.s, those of B α U 0 lie in C 0 pDq and thus the function x Þ ÝÑ VarpB α U 0 pxqq also lies in C 0 pDq ([5], chapter 1, Section 4.3). Therefore the function x Þ ÝÑ σ B α U 0 pxq also lies in C 0 pDq. We now check that the integral I below is finite. We use Tonelli's theorem and the Cauchy-Schwarz inequality:

I :" ż Ω ÿ |α|ď1 ż R 3 ˇˇB α U 0 ω px ´yq ˇˇ|µ t α |pdyq ÿ |α 1 |ď1 ż R 3 ˇˇB α 1 U 0 ω px 1 ´y1 q ˇˇ|µ t 1 α 1 |pdy 1 qPpdωq " ÿ |α|,|α 1 |ď1 ż R 3 ż R 3 ż Ω ˇˇB α U 0 ω px ´yqB α 1 U 0 ω px 1 ´y1 q ˇˇPpdωq|µ t α |pdyq|µ t 1 α 1 |pdy 1 q " ÿ |α|,|α 1 |ď1 ż R 3 ż R 3 E " |B α U 0 px ´yqB α 1 U 0 px 1 ´y1 q| ‰ |µ t α |pdyq|µ t 1 α 1 |pdy 1 q ď ÿ |α|,|α 1 |ď1 ż R 3 ż R 3 ˆE" B α U 0 px ´yq 2 ‰ E " B α 1 U 0 px 1 ´y1 q 2 ‰ ˙1{2 |µ t α |pdyq|µ t 1 α 1 |pdy 1 q ď ˜ÿ |α|ď1 ż R 3 ˆE" B α U 0 px ´yq 2 ‰ ˙1{2 |µ t α |pdyq ¸ˆ˜ÿ |α|ď1 ż R 3 ˆE" B α U 0 px ´yq 2 ‰ ˙1{2 |µ t 1 α |pdyq ḑ´ÿ |α|ď1 p|µ t α | ˚σB α U 0 qpxq ¯ˆ´ÿ |α|ď1 p|µ t 1 α | ˚σB α U 0 qpx 1 q ¯ă `8.
For all multi-index α, the scalar p|µ t α | ˚σB α U 0 qpxq is finite because x Þ ÝÑ σ B α U 0 pxq is continuous and |µ t α | is compactly supported. Note also that from Assumption pA 2 q, the GP U 0 is mean square differentiable up to order 1, which implies ( [START_REF] Ritter | Average-case analysis of numerical problems[END_REF], Section III.1.4) that we have, for all multi-indexes α, α 1 such that |α|, |α 1 | ď 1, x and x 1 :

E " B α U 0 pxqB α 1 U 0 px 1 q ‰ " B α 1 B α 1 2 k u px, x 1 q. ( 50
)
where B 1 (respectively B 2 ) denotes derivatives with reference to the first (respectively second) argument of k u . We may thus permute integrals and differential operators in E " U pzqU pz 1 q ‰ :

E " U pzqU pz 1 q ‰ " E « ÿ |α|ď1 ż R 3 p´1q |α| B α U 0 px ´yqµ t α pdyqq ÿ |α 1 |ď1 ż R 3 p´1q |α 1 | B α 1 U 0 px ´yqµ t 1 α 1 pdy 1 q ff " ÿ |α|,|α 1 |ď1 ż R 3 ż R 3 p´1q |α| p´1q |α 1 | B α 1 B α 1 2 E " U 0 px ´yqU 0 px 1 ´y1 q ‰ µ t α pdyqµ t 1 α 1 pdy 1 q " ÿ |α|,|α 1 |ď1 ż R 3 ż R 3 p´1q |α| p´1q |α 1 | B α 1 B α 1 2 k u px ´y, x 1 ´y1 qµ t α pdyqµ t 1 α 1 pdy 1 q " " ´ÿ |α|ď1 B α µ t α b ÿ |α 1 |ď1 B α 1 µ t 1 α 1 ¯˚k u ȷ px, x 1 q " rp 9 F t b 9 F t 1 q ˚ku spx, x 1 q,
which proves [START_REF] Rudin | Functional analysis[END_REF]. One proves that E " V pzqV pz 1 q ‰ " rpF t b F t 1 q ˚kv spx, x 1 q the exact same way, which is actually simpler as F t is directly a measure. To conclude,

k W pz, z 1 q " CovpW pzq, W pz 1 qq " ErpW pzqW pz 1 qs " E "`U pzq `V pzq ˘`U pz 1 q `V pz 1 q ˘‰ " E " U pzqU pz 1 q ‰ `E" U pzqV pz 1 q ‰ `E" V pzqU pz 1 q ‰ `E" V pzqV pz 1 q ‰ " rp 9 F t b 9 F t 1 q ˚ku spx, x 1 q `rpF t b F t 1 q ˚kv spx, x 1 q. (51)
The cross terms are null because U pzq and V pz 1 q are independent as well as U pz 1 q and V pzq. piiq : with expression (45), one checks that for any fixed z 1 , the function z Þ ÝÑ k W pz, z 1 q is of the form [START_REF] Narcowich | Generalized Hermite interpolation via matrix-valued conditionally positive definite functions[END_REF] and thus verifies lk x 1 " 0 in the sense of distributions. piiq is then a direct consequence of Proposition 3.5.

Remark 4.2. If U and V are not independent, then the two terms rp 9 F t b F t 1 q ˚kuv spx, x 1 q and rpF t b 9 F t 1 q ˚kvu spx, x 1 q must be added to equation [START_REF] Roques | Spatial statistics and stochastic partial differential equations: A mechanistic viewpoint[END_REF], where k uv px, x 1 q denotes the cross covariance between U and V : k uv px, x 1 q " CovpU pxq, V px 1 qq and k vu px, x 1 q " CovpV pxq, U px 1 qq " k uv px 1 , xq.

More explicitly, we have the following Kirschoff-like integral formulas for k wave v and k wave u :

rpF t b F t 1 q ˚kv spx, x 1 q " tt 1 ż Sp0,1qˆSp0,1q k v px ´c|t|γ, x 1 ´c|t 1 |γ 1 q dΩdΩ 1 p4πq 2 , ( 52
) rp 9 F t b 9 F t 1 q ˚ku spx, x 1 q " ż Sp0,1qˆSp0,1q
´ku px ´c|t|γ, x 1 ´c|t 1 |γ 1 q ´c|t|∇ 1 k u px ´c|t|γ, x 1 ´c|t 1 |γ 1 q ¨γ ´c|t 1 |∇ 2 k u px ´c|t|γ, x 1 ´c|t 1 |γ 1 q ¨γ1 `c2 tt 1 γ T ∇ 1 ∇ 2 k u px ´c|t|γ,

x 1 ´c|t 1 |γ 1 qγ 1 ¯dΩdΩ 1 p4πq 2 . ( 53
)
Above, ∇ 1 k u px, x 1 q is the gradient vector of k u with reference to x, ∇ 2 k u px, x 1 q is the gradient vector of k u with reference to x 1 and ∇ 1 ∇ 2 k u px, x 1 q is the matrix whose entry pi, jq is given by

∇ 1 ∇ 2 k u px, x 1 q ij " B x 1 i B x 2 j k u px, x 1 q. ( 54
)
(B x 1 i is the derivative with reference to the i th coordinate of x, B x 2 j is the derivative with reference to the j th coordinate of x 1 ). to initial conditions u 0 and v 0 with piecewise regularity. The formulas [START_REF] Serre | Systems of conservation laws[END_REF] and [START_REF] Steinwart | Convergence types and rates in generic Karhunen-Loeve expansions with applications to sample path properties[END_REF] are valid in a more general context than that of assumptions (A 1 ) and (A 2 ). We provide below examples where these formulas yield valid covariance functions (in particular, functions defined for all values of px, tq and px 1 , t 1 q) corresponding to initial conditions with some forms of piecewise discontinuities. Assume, for example, that the initial speed v 0 is compactly supported on a ball Bpx 0 , Rq centered on some point x 0 with radius R. This is a natural model when v 0 is assumed to be a localized source. For the process V 0 , this translates as V 0 pxq " 0 a.s. if x is outside the ball Bpx 0 , Rq. One can thus truncate the covariance function of V 0 accordingly, e.g. choosing the following function for k v (see Section 2.1.3 for k 1{2 ) k v px, x 1 q " k 1{2 px, x 1 q1 r0,Rs p||x ´x0 ||q1 r0,Rs p||x 1 ´x0 ||q.

(

) 55
Above, ||x|| denotes the Euclidean norm of x. Such a covariance function indeed verifies k v px, xq " VarpV 0 qpxqq " 0 if ||x ´x0 || ą R and the GP corresponding to k v is V 0 pxq " V 1{2 pxq1 r0,Rs p||x x0 ||q, where V 1{2 is a continuous modification of a GP with covariance function k 1{2 . Note that the sample paths of V 0 are piecewise continuous and that V as defined in ( 42) is well-defined and measurable. The integrals in [START_REF] Serre | Systems of conservation laws[END_REF] still make sense and point piiq from Proposition is still valid: the sample paths of the process V whose covariance function is k wave v (or any other measurable centered second order random field with this covariance function) remains a solution of the wave equation in the distributional sense. One can perform the same kind of discussions on k wave u : for example, equation [START_REF] Steinwart | Convergence types and rates in generic Karhunen-Loeve expansions with applications to sample path properties[END_REF] shows that when k u P C 1,1 pR 3 ˆR3 qzC 2,2 pR 3 ˆR3 q, k wave u is only expected to lie in C 1,1 pR 3 ˆR3 q; the sample paths of the GP with covariance function k wave u will be at most of class C 1 and thus cannot be strong solutions of equation [START_REF] Mendes | Bayesian inference in the numerical solution of Laplace's equation[END_REF]. This is the case when k u is the k 3{2 Matérn covariance function from equation [START_REF] Carrizo-Vergara | A general framework for SPDE-based stationary random fields[END_REF].

More generally, one can incorporate a finite number of discontinuities on k v and on the derivatives of k u so that they remain piecewise continuous: the integrals above will remain well defined and the sample paths of the corresponding GPs will remain distributional solutions to the wave equation, even though they will not be sufficiently differentiable to be strong solutions.

Conclusion and perspectives

In Section 3, we have presented a new result that provides a simple characterization of the measurable second order random fields pU pxqq xPD whose sample paths verify homogeneous linear differential constraints within the framework of generalized functions. This characterization is valid for any linear differential operator L, provided that its coefficients fulfil minimal smoothness requirements, and no stationarity assumptions over pU pxqq xPD are required. Motivated by physical applications, we described in Section 4 a Gaussian process model of the wave equation which is central to describe propagation phenomena. This PDE served as an application case for Proposition 3.5, and the GP model was derived by putting a GP prior on the wave equation's initial conditions. In Proposition 4.1, we presented covariance formulas that are tailored to the wave equation and take the form of convolutions; these expressions are interesting in themselves and call for physics-informed GPR applications for this equation. In particular, we showed that these formulas can model piecewise continuously differentiable solutions for the wave equation. Moreover, this setting provides a natural way to incorporate any type of information, both numerical or experimental. In a forthcoming paper, we will show how to use GPR conjointly with the covariance functions from Proposition 4.1 in a numerical setting, in order to construct approximate solutions of the wave equation based on scattered observations. This in turn provides a natural method for solving different inverse problems, by using the likelihood of GPR as well as the reconstructed solution. An other application concerns the design of transparent boundary conditions (TBC): this provides artificial boundary conditions on a computational domain so that the computed solution is exactly an approximation of a solution on the whole space. Those conditions are usually nonlocal and restricted to simple geometries. A GPR strategy is meshless therefore suitable to design TBC on any type of computational domains. Proposition 3.5 constitutes a first step towards understanding PDE constrained random fields in an weakened sense; different functional analysis frameworks can now be considered, obvious extensions being the weak or variational formulations of equation [START_REF] Adler | The Geometry of Random Fields[END_REF]. These formulations are obtained by transferring only a part of the derivatives of the PDE to the test function and are for instance the canonical way of studying elliptic PDEs ( [START_REF] Evans | Partial differential equations[END_REF], Section 6.1.2). The natural spaces arising from these formulations are Sobolev spaces rather than D 1 pDq. An attached question, as studied in [START_REF] Scheuerer | Regularity of the sample paths of a general second order random field[END_REF], is that of the Sobolev regularity of a given second order random field; a current research topic is whether or not one may relax the continuity assumptions required in [START_REF] Scheuerer | Regularity of the sample paths of a general second order random field[END_REF]. Finally, the matter of using random fields for modelling and approximating solutions of nonlinear PDEs is a natural direction for future research.

  ˙ˆÿ |α|ďn ||B α φ|| 8 ă `8.

4. 1

 1 General solution to the 3 dimensional wave equation Denote the 3D Laplace operator ∆ " B 2 xx `B2 yy `B2 zz and the d'Alembert operator l " 1{c 2 B 2 tt ´∆ with constant wave speed c ą 0. We focus on the general initial value problem in the free space R 3

4. 2 . 1

 21 Extending the covariance functions k wave u and k wave v

  The solution of this problem is unique in the distributional sense ([START_REF] Duistermaat | Distributions[END_REF], p. 164). It can be extended to all t P R ([14], p. 295) and is represented as follow ([START_REF] Duistermaat | Distributions[END_REF], p. 295 again)

	wpx, tq " pF t ˚v0 qpxq `p 9 F t ˚u0 qpxq	@px, tq P R 3 ˆR,

  t thanks to Radon measures (Sections 2.2.5 and 2.2.6) and use Fubini's theorem. We refer to Sections 2.2.6 and 2.2.8 for the definition of 9 F t b 9 F t 1 , and Section 2.2.7 for the definition of p 9 F t b 9 F t 1 q ˚ku .

Acknowledgements. Research of all the authors was supported by SHOM (Service Hydrographique et Océanographique de la Marine) project "Machine Learning Methods in Oceanography" no-20CP07. We thank Rémy Baraille in particular for his personal involvement in the project. We are thankful to the reviewers and the editor for their interesting and constructive remarks, leading to an improved and enriched version of the manuscript.

where σ R is the surface measure of the sphere of center 0 and radius R; 9 F t " B t F t means that for all f P C 1 c pR 3 q, x 9 F t , f y " B t xF t , f y. We make these expressions more explicit in equation [START_REF] Owhadi | Bayesian numerical homogenization[END_REF], using spherical coordinates. It is worth noting that pF t q tPR corresponds to the Green's function of the wave equation, in the sense that it verifies the system (37) with u 0 " 0 and v 0 " δ 0 where δ 0 is the Dirac mass ( [START_REF] Duistermaat | Distributions[END_REF], pp. 294-295). As discussed in Remark 3.4, pF t q tPR is a family of singular measures and this PDE system has to be understood in the distributional sense. Note also that equations [START_REF] Nguyen | Gaussian functional regression for linear partial differential equations[END_REF] show that F t and 9

F t are supported on the sphere of radius c|t|: "the support of F t propagates at finite speed c". This property is known as the Huygens principle for the three dimensional wave equation, see [START_REF] Evans | Partial differential equations[END_REF], p. 80.

Suppose that u 0 P C 1 pR 3 q and v 0 P C 0 pR 3 q, then w as defined in equation ( 38) is a pointwise defined function (Section 2.2.7) and in that case an explicit formula for such convolutions is reminded in equation ( 15) (yet one may actually make sense out of (38) when u 0 and v 0 are only required to be any generalized functions, see [START_REF] Trèves | Topological vector spaces, distributions and kernels[END_REF], Chapter 27).

Equation ( 38) can be written using means over spheres. Denote pr, θ, ϕq, r ě 0, θ P r0, πs, ϕ P r0, 2πs the spherical coordinates, Sp0, 1q the unit sphere of R 3 and γ " psin θ cos ϕ, sin θ sin ϕ, cos θq T the corresponding parametrization of Sp0, 1q (||γ|| 2 " 1). We write dΩ " sin θdθdϕ the surface differential element of Sp0, 1q. The formulas [START_REF] Narcowich | Generalized Hermite interpolation via matrix-valued conditionally positive definite functions[END_REF]

Gaussian process modelling of the solution

Suppose now that u 0 and v 0 are unknown, and only pointwise values of w are observed. We thus model u 0 and v 0 as random functions and put Gaussian process priors over u 0 and v 0 . More precisely, we make the following assumptions.

(A 1 ) Suppose that the initial conditions u 0 and v 0 of Problem [START_REF] Mendes | Bayesian inference in the numerical solution of Laplace's equation[END_REF] are sample paths drawn from two independent Gaussian processes U 0 " GP p0, k u q and V 0 " GP p0, k v q: Dω P Ω, @x P R 3 , u 0 pxq " U 0 ω pxq and v 0 pxq " V 0 ω pxq.

(A 2 ) Suppose that all sample paths of U 0 lie in C 1 pR 3 q and that those of V 0 lie in C 0 pR 3 q, almost surely. A sufficient condition for this is given in [START_REF] Adler | Random Fields and Geometry[END_REF], Theorem 1.4.2. This theorem states that under mild technical assumptions, the paths of pU pxqq xPD " GP p0, kq lie in C l a.s. as soon as k P C 2l pD ˆDq, which we assume from now on. This is e.g. fulfilled by the Matérn covariance functions from equation ( 8), with l " 0 for k 1{2 and l " 1 for k 3{2 .

We now analyse the consequence of these two assumptions. First, they imply that by solving [START_REF] Mendes | Bayesian inference in the numerical solution of Laplace's equation[END_REF], one obtains a time-space stochastic process W px, tq defined by

Here again, V 0 ω denotes the sample path of V 0 at ω P Ω and likewise for U 0 ω . In particular, thanks to assumption pA 2 q, equation (41) defines a random variable for all px, tq. Note the space-time variable z " px, tq and note the random variables V pzq : ω Þ ÝÑ pF t ˚V 0 ω qpxq and U pzq : ω Þ ÝÑ p 9

that is, W pzq " U pzq `V pzq. We show in the next proposition that the random fields U, V and W are GPs as well. In particular we describe their covariance functions.