Solutions of Hyperbolic Stochastic PDEs on Bounded and Unbounded Domains

We treat several classes of hyperbolic stochastic partial differential equations in the framework of white noise analysis, combined with Wiener–Itô chaos expansions and Fourier integral operator methods. The input data, boundary conditions and coefficients of the operators are assumed to be generalized stochastic processes that have both temporal and spatial dependence. We prove that the equations under consideration have unique solutions in the appropriate Sobolev–Kondratiev or weighted-Sobolev–Kondratiev spaces. Moreover, an explicit chaos form of the solutions is obtained, useful for calculating expectations, variances and higher order moments of the solution.


Introduction
Hyperbolic stochastic partial differential equations arise as models of various phenomena used in mathematical physics, economy, molecular biology and many other areas of science, where random fluctuations and uncertainties are incorporated into the equation by white noise or other singular generalized stochastic processes such as Poissonian processes or general Lévy processes. In this paper we will consider two types of hyperbolic problems for suitable differential operators, acting by the Wick product instead of classical multiplication. This is due to the fact that we allow random terms to be present both in the initial conditions and right-hand side of the equations, as well as in the coefficients of the involved operators. Having all these highly random terms will lead to singular solutions that do not allow to use ordinary multiplication, but require its renormalization also known as the Wick product. The Wick product is known to represent the highest order stochastic approximation of the ordinary product [32], and has been used in many models together with the Wiener chaos expansion method [18,19,26,27,29,30,[40][41][42]47].
Powerful tools used for deterministic equations with singular input data are pseudodifferential calculus and Fourier integral operators (see, e.g., [11,22,33,34]), that have experienced a rapid development in the recent years (see, for instance, [1,3,4,10,44] and the references quoted therein). The roots of pseudo-differential operators and Fourier integral operators stem from microlocal analysis (see, e.g., [17,20]). General approaches to solving deterministic hyperbolic equations are presented, e.g., in [22,28], and we will rely on these results and their extensions.
Henceforth, in this paper we will present techniques for solving singular hyperbolic stochastic partial differential equations resulting from the synergy of these, nowadays classical, two powerful methods: chaos expansions and Fourier integral operators.
The first model we will consider is an initial-boundary value problem for a second order, wave-type, differential operator on a bounded open set U ⊂ R d having smooth boundary, and on a time interval I = [0, T ], T > 0, namely a jk (t, x; ω)∂ k , having spacetime stochastic processes a jk as coefficients that are continuously differentiable in time, and satisfying suitable ellipticity conditions, which will be specified in detail in Sects. 2 and 3 below. Here, chaos expansions will be used in connection with well-known representations and estimates for an infinite sequence of associated deterministic initial-boundary value problems on the bounded open set U ⊂ R d . The second model on which we will focus is an initial value (that is, Cauchy) problem for a differential hyperbolic operator of order m ∈ N, which we will study globally on R d , namely, L♦u(t, x; ω) = f (t, x; ω) (t, x; ω) ∈ I × R d × (D k t u)(0, x; ω) = u k (x; ω), k = 0, . . . , m − 1, (x; ω) ∈ R d × , and a jα being smooth space-time stochastic processes, j = 1, . . . , m, α ∈ N d 0 such that |α| ≤ j (see Sects. 2 and 4 for the precise hypotheses on L, f , and u k , k = 0, . . . , m −1). To perform our analysis in this case, we will again use chaos expansions, but this time in connection with the properties of a class of Fourier integral operators, defined through objects globally defined on R d .
Hyperbolic SPDEs via Wiener chaos expansion methods have been studied in [21], but our approach is more powerful and allows more singular input data in the model. The main idea we present in this paper relies on the chaos expansion method (also known as the propagator method): first, one uses the chaos expansion of all stochastic data in the equation to convert the SPDE into an infinite system of deterministic PDEs, then the PDEs are recursively solved, and finally one must sum up these solutions to obtain the chaos expansion form of the solution of the initial SPDE. The crucial point is to prove convergence of the series given by the chaos expansion that defines the solution, and this part relies on obtaining good energy estimates of the PDE solutions, proving their regularity and using estimates on the Wick products. This approach has many advantages, most notably it provides an explicit form of the solution of the SPDE from which one can directly compute the expectation, variance and other moments, and it is convenient also for numerical approximations by truncating the series in the chaos expansion to finite sums.
The second main tool we use in this paper is the SG calculus of Fourier integral operators (further abbreviated as SG-FIOs theory). Applications of the SG-FIOs theory to SG-hyperbolic Cauchy problems were initially given in [12,14]. Many authors have, since then, expanded the SG-FIOs theory and its applications to the solution of hyperbolic problems in various directions. To mention a few, see, e.g., [10,44], and the references quoted there and in [4]. In particular, the results in Theorem 2.11 have been applied in [4] to study classes of SG-hyperbolic Cauchy problems, constructing their fundamental solution {E(t, s)} 0≤s≤t≤T . The existence of the fundamental solution provides, via Duhamel's formula, existence and uniqueness of the solution to the system, for any given Cauchy data in the weighted Sobolev spaces H z,ζ (R d ), (z, ζ ) ∈ R 2 . In short, the Cauchy problem for a linear SG-hyperbolic operator L, satisfying suitable additional hypotheses (which will be explicitly stated in Sect. 4), can be turned into an equivalent Cauchy problem for a first order system. The fundamental solution operator to such a system allows then to write the (unique) solution of the original Cauchy problem in terms of SG-FIOs (modulo remainder terms). A remarkable feature, typical for these classes of hyperbolic problems, is the well-posedness with loss of decay/gain of growth at infinity, observed, e.g., in [2,3,12]. We also mention that random-field solutions of hyperbolic SPDEs via Fourier integral operator methods have been recently studied in [5,8], while function-valued solutions for associated semilinear hyperbolic SPDEs have been obtained in [7].
The plan of our exposition is as follows. In the preliminary section (Sect. 2) we provide a basic overview of the notation and recall some results that are required for further reading. We introduce the basic notions of white noise theory including chaos expansions of generalized stochastic processes, Wick products and stochastic differential operators, and we recall the fundamental notions of pseudo-differential calculus, Fourier integral operators, and weighted Sobolev spaces, within the environment of the so-called SG calculus. In Sect. 3, which represents the first main result of the paper, we prove existence and uniqueness of a local solution to the linear equation (1.1). In Sect. 4, which represents the second main result of the paper, we prove existence and uniqueness of a global solution to the Eq. (1.2) and then prove existence and uniqueness of solutions to systems of first order linear hyperbolic SPDEs of the form (1.2). Finally, Sect. 5 contains some examples of applications of the previous results, including the modeling of seismic wave equations and earthquake motions [37,46], of molecular dynamics [15,16,39,48], and of plasma reactions and dynamics of the inner structure of stars [15,43].
Although our exposition follows the classical Hida-Kondratiev space approach within classical white noise theory, the chaos expansion method we present can easily be extended to model fractional Brownian motion with a Hurst parameter H ∈ (0, 1), fractional Poissonian noise or other fractional versions of stochastic processes. In [19] it was shown that there exists a unitary mapping between Gaussian and Poissonian white noise spaces. Hence, solutions of a stochastic differential equation on the Poissonian white noise space can be obtained by applying this mapping to the solution of the corresponding stochastic differential equation taken on the Gaussian white noise space. Fractional versions of spaces of this type were further studied in [24,25] and can be related to the Gaussian white noise space by a simple change of the Hermite basis (2.1) to another basis of orthogonal polynomials (e.g. Charlier polynomials for the Poissonian noise).

Preliminaries
In this section we recall various notions, which we will use in the sequel.

Chaos Expansions and the Wick Product
Denote by ( , F, P) the Gaussian white noise probability space (S (R), B, μ), where S (R) denotes the space of tempered distributions, B the Borel sigma-algebra generated by the weak topology on S (R) and μ the Gaussian white noise measure corresponding to the characteristic function given by the Bochner-Minlos theorem. We recall the notions related to L 2 ( , μ) (see [19]), where = S (R) and μ is Gaussian white noise measure. We adopt the notation N 0 = {0, 1, 2, . . . }, Define the set of multi-indices I to be (N N 0 ) c , i.e. the set of sequences of non-negative integers which have only finitely many nonzero components. Especially, we denote by 0 = (0, 0, 0, . . .) the multi-index with all entries equal to zero. The length of a multi-index is |α| = ∞ i=1 α i for α = (α 1 , α 2 , . . .) ∈ I, and it is always finite. Similarly, α! = ∞ i=1 α i !, and all other operations are also carried out componentwise. We will use the convention that α − β is defined if α n − β n ≥ 0 for all n ∈ N, i.e., if α − β ≥ 0, and leave α − β undefined if α n < β n for some n ∈ N. We here denote by h n , n ∈ N 0 , the Hermite orthogonal polynomials and by ξ n , n ∈ N, the Hermite functions The Wiener-Itô theorem states that one can define an orthogonal basis {H α } α∈I of L 2 ( , μ), where H α are constructed by means of Hermite orthogonal polynomials h n and Hermite functions ξ n , (2.1) Then, every F ∈ L 2 ( , μ) can be represented via the so called chaos expansion Denote by ε k = (0, 0, . . . , 1, 0, 0, . . .), k ∈ N the multi-index with the entry 1 at the kth place. Denote by H 1 the subspace of L 2 ( , μ), spanned by the polynomials H ε k (·), k ∈ N. All elements of H 1 are Gaussian stochastic processes, e.g. the most prominent one is Brownian motion given by the chaos expansion Denote by H m the mth order chaos space, i.e. the closure of the linear subspace spanned by the orthogonal polynomials H α (·) with |α| = m, m ∈ N 0 . Then the Wiener-Itô chaos expansion states that L 2 ( , μ) = ∞ m=0 H m , where H 0 is the set of constants in L 2 ( , μ). The expectation of a random variable is its orthogonal projection onto H 0 , hence it is given by E(F(ω)) = f (0,0,...) .
It is well-known that the time-derivative of Brownian motion (white noise process) does not exist in the classical sense. However, changing the topology on L 2 ( , μ) to a weaker one, T. Hida [18] defined spaces of generalized random variables containing the white noise as a weak derivative of the Brownian motion. We refer to [18,19,23] for white noise analysis (as an infinite dimensional analogue of the Schwartz theory of deterministic generalized functions).
. . , α n , . . .) ∈ I. We will often use the fact that the series α∈I (2N) − pα converges for p > 1 [19,Proposition 2.3.3]. Define the Banach spaces Their topological dual spaces are given by The Kondratiev space of generalized random variables is (S) −1 = p∈N 0 (S) −1,− p endowed with the inductive topology. It is the strong dual of (S) 1 = p∈N 0 (S) 1, p , called the Kondratiev space of test random variables which is endowed with the projective topology. Thus, The time-derivative of the Brownian motion exists in the generalized sense and belongs to the Kondratiev space (S) −1,− p for p > 5 12 [23, p. 21]. We refer to it as to white noise and its formal expansion is given by . We extended in [40] the definition of stochastic processes also to processes of the chaos expansion form U (t, ω) = α∈I u α (t)H α (ω), where the coefficients u α are elements of some Banach space X . We say that U is an X -valued generalized stochastic process, i.e.
The notation ⊗ is used for the completion of a tensor product with respect to the π -topology (see [50]). We note that if one of the spaces involved in the tensor product is nuclear, then the completions with respect to the π -and the ε-topology coincide. It is known that (S) 1 and (S) −1 are nuclear spaces [19,Lemma 2.8.2], thus in all forthcoming identities ⊗ can be equivalently interpreted as the ⊗ π -or ⊗ ε -completed tensor product. Thus, when dealing with the tensor products with (S) 1, p and (S) −1,− p , we work with the π -topology.
The Wick product of two stochastic processes F = α∈I f α H α and G = β∈I g β H β ∈ X ⊗ (S) −1 is given by and the nth Wick power is defined by F ♦n = F ♦(n−1) ♦F, F ♦0 = 1. Note that H nε k = H ♦n ε k for n ∈ N 0 , k ∈ N. The Wick product always exists and results in a new element of X ⊗ (S) −1 , moreover it exhibits the property of E(F♦G) = E(F)E(G) holding true. The ordinary product of two generalized stochastic processes does not always exist and E(F · G) = E(F)E(G) would hold only if F and G were uncorrelated.
One particularly important choice for the Banach space X is X = C k [0, T ], k ∈ N. We proved in [41] that differentiation of a stochastic process can be carried out componentwise in the chaos expansion, i.e. due to the fact that (S) −1 is a nuclear space it holds that C k ([0, T ], (S) −1 ) = C k [0, T ]⊗(S) −1 . This means that a stochastic process U (t, ω) is k times continuously differentiable if and only if all of its coefficients The same holds for Banach space valued stochastic processes i.e. elements of C k ([0, T ], X ) ⊗ (S) −1 , where X is an arbitrary Banach space. By the nuclearity of (S) −1 , these processes can be regarded as elements of the tensor product space In order to solve (1.1) and (1.2) we will choose some special Banach spaces, for example if U is an open subset of R d , then some appropriate choices may be etc. depending on the input data in the SPDEs.
In general, the function spaces that we will adopt as those where to look for the solutions to (1.1) and (1.2) will be of the form where I ⊂ R is an interval of the form [0, T ] or [0, ∞), and G k , k = 0, 1, 2, . . . , l, or k ∈ Z + , are suitable Hilbert spaces (or Banach spaces) such that where → denotes dense continuous embeddings. We can also consider the topological duals of G j , j ∈ Z + , denoted by G − j , respectively, and write For example, if G 0 = L 2 (U ), G 1 = H 1 0 (U ), then G −1 = H −1 (U ) and they form a Gelfand triple. Hence, one has In particular, for the spaces in (2.2) and in (2.3) we have, respectively,

Stochastic Operators and Differential Operators with Stochastic Coefficients
Let X be a Banach space endowed with the norm · X . Consider X ⊗ (S) −1 with elements u = α∈I u α H α so that α∈I u α 2 X (2N) − pα < ∞ for some p ≥ 0. Let D ⊂ X be a dense subset of X endowed with the norm · D and A α : D → X , α ∈ I, be a family of linear operators on this common domain D. Assume that each A α is bounded i.e., In case when D = X , we will write L(X ) instead of L(D, X ). The family of operators A α , α ∈ I, gives rise to a stochastic operator A♦ : D ⊗(S) −1 → X ⊗(S) −1 , that acts in the following manner In the next two lemmas we provide two sufficient conditions that ensure the stochastic operator A♦ to be well-defined. Both conditions rely on the l 2 or l 1 bounds with suitable weights. They are actually equivalent to the fact that A α , α ∈ I, are polynomially bounded, but they provide finer estimates on the stochastic order (Kondratiev weight) of the domain and codomain of A♦.
Proof For u ∈ X ⊗ (S) −1,−r , we have by the generalized Minkowski inequality that A similar example may be constructed with D = L 2 (R) and X = H −1 (R). Note that in these examples, we could have written the operator also in the form Let us now consider the differential operator that governs Eq.
This operator acts in the following way: for an element Hence, we may identify the operator A♦ with the family of operators A α : and thus the operator A♦ given by is well-defined by Lemma 2.1.

Lemma 2.3 In particular, if the operator has deterministic coefficients, i.e. ifÃ is of the formÃ
Proof Clearly, we may identifyÃ with the constant net of operators for all u(t, x; ω) = α∈I u α (t, x)H α (ω) and hence we obtainÃ : C l (I , In Sect. 3 we will assume other types of conditions on the operator A♦; in particular if we want better regularity of the solutions, e.g. u ∈ C l (I , L 2 (U )) ⊗ (S) −1 instead of u ∈ C l (I , H −1 (U )) ⊗ (S) −1 , then we must assume that some of its components A α , α ∈ I are differential operators only of order one. Precise conditions will be provided in Sect. 3.
Considering the differential operator L that governs Eq. (2.2), we will make special choices for the domain D and range X involving the so called weighted Sobolev spaces H z,ζ (R d ) and many other types of spaces that stem from pseudodifferential calculus. These will be introduced in the next section.

The Global SG Calculus of Pseudodifferential and Fourier Integral Operators
We here recall some basic definitions and facts about the SG-calculus of pseudodifferential and Fourier integral operators, through standard material appeared, e.g., in [4] and elsewhere (sometimes with slightly different notational choices). In the sequel we will often use the so-called Japanese bracket of y ∈ R d , given by y = 1 + |y| 2 .
The class S m,μ = S m,μ (R d ) of SG symbols of order (m, μ) ∈ R 2 is given by all the functions a(x, ξ) ∈ C ∞ (R d × R d ) with the property that, for any multiindices hold (see [11,31,38]). For m, μ ∈ R, ∈ N 0 , is a family of seminorms, defining the Fréchet topology of S m,μ . The corresponding classes of pseudodifferential operators Op(S m,μ ) = Op(S m,μ (R d )) are given by The operators in (2.5) form a graded algebra with respect to composition, i.e., which implies that the symbol c equals a · b modulo S m 1 +m 2 −1,μ 1 +μ 2 −1 . Note that For any a ∈ S m,μ , (m, μ) ∈ R 2 , Op(a) is a linear continuous operator from S(R d ) to itself, extending to a linear continuous operator from S (R d ) to itself, and from (here D ζ is understood as a pseudodifferential operator) with the naturally induced Hilbert norm. When z ≥ z and ζ ≥ ζ , the continuous embedding H z,ζ → H z ,ζ holds true. It is compact when z > z and ζ > ζ . Since as well as, for the space of rapidly decreasing distributions, see [6] and [45, Chap. VII, § 5], The continuity property of the elements of Op(S m,μ ) on the scale of spaces H z,ζ (R d ), (m, μ), (z, ζ ) ∈ R 2 , is expressed more precisely in the next theorem.
where [t] denotes the integer part of t ∈ R.

Remark 2.6 1. Trivially, any
The following characterization of the class O(−∞, −∞) is often useful. An operator A = Op(a) and its symbol a ∈ S m,μ are called elliptic (or S m,μ - for some constant C > 0. If R = 0, a −1 is everywhere well-defined and smooth, and where I denotes the identity operator. In such a case, A turns out to be a Fredholm operator on the scale of functional spaces H z,ζ , We now recall the class of SG-phase functions. A real valued function ϕ ∈ C ∞ (R 2d ) belongs to the class P of SG-phase functions if it satisfies the following conditions: For any a ∈ S m,μ , (m, μ) ∈ R 2 , ϕ ∈ P, the SG FIOs are defined, for u ∈ S(R n ), as and Here the operators Op ϕ (a) and Op * ϕ (a) are sometimes called SG FIOs of type I and type II, respectively, with symbol a and (SG-)phase function ϕ. Note that a type II operator satisfies Op * ϕ (a) = Op ϕ (a) * , that is, it is the formal L 2 -adjoint of the type I operator Op ϕ (a).
The following theorem summarizes composition results between SG pseudodifferential operators and SG FIOs of type I that we are going to use in the present paper, see [13] for proofs and composition results with SG FIOs of type II.
To consider the composition of SG FIOs of type I and type II some more hypotheses are needed, leading to the definition of the classes P δ and P δ (λ) of regular SG-phase functions (see [13]). We recall now their definition. Let λ ∈ [0, 1) and δ > 0. A function ϕ ∈ P belongs to the class P δ (λ) if it satisfies the following conditions: If only condition (1) holds, we write ϕ ∈ P δ . Let ∈ N 0 . In [22] the following seminorms are defined: We define the following subclass of the class of regular SG phase functions: Let λ ∈ [0, 1), δ > 0, ≥ 0. Recall [22]: A function ϕ belongs to the class P δ (λ, ) if ϕ ∈ P δ (λ) and J ≤ λ for the corresponding J .

Remark 2.10
Similarly to Theorem 2.4, chosen a ∈ S m,μ (R d ), (m, μ) ∈ R 2 , and a regular SG phase function ϕ, for any (z, ζ ) ∈ R 2 there exists a constant C > 0, depending only on d, m, μ, z, ζ , such that where p dmμ depends continuously on a finite collection of seminorms of a and ϕ, as symbols in S m,μ (R d ) and S 1,1 (R d ), respectively.
Let M ≥ 2. The study of the composition of SG FIOs of type I Op ϕ j (a j ) with regular SG-phase functions ϕ j ∈ P δ (λ j ) and symbols a j ∈ S m j ,μ j , j = 1, . . . , M, has been done in [4]. The result of such composition is still an SG-FIO with a regular SG-phase function ϕ given by the so-called multi-product ϕ 1 · · · ϕ M of the phase functions ϕ j , j = 1, . . . , M, and symbol a as in Theorem 2.11 here below, a corollary of the main Theorem in [4].

Solutions on Bounded Domains
First we consider the local problem (1.1), i.e., the problem on a bounded spatial domain.
By decomposing u, f , a jk into their chaos expansions, and using the properties of the Wick product, we find that (1.1) is equivalent to an infinite sequence of systems: We also decompose the initial values: Thus (1.1) reduces to the so called propagator system together with the initial and boundary value conditions The initial boundary value problems given by (3.2) together with (3.3) can be solved recursively by induction over the length of the multiindex γ . To this aim, we will use the following well-known result, see [ Assume that a(t; u, v), t ∈ I , is a family of continuous sesquilinear forms on V , such that Consider the Cauchy problem If f ∈ L 2 (I , H ), u 0 ∈ V , u 1 ∈ H , then, there exists a unique u solution of (3.7) such that u ∈ L 2 (I , V ), ∂ t u ∈ L 2 (I , H ), ∂ 2 t u ∈ L 2 (I , V ). Moreover, the energy inequality holds true for any t ∈ I , for a suitable C > 0 depending only on A(t) and T .
Notice that, in addition to Theorem 3.1, after possibly a modification on a set of measure zero, the next result also holds true, see again, e.g., [  Moreover, the mapping ( f , u 0 , u 1 ) → (u, ∂ t u), associating the right-hand side and initial data of (3.7) to its solution and first t-derivative, is a well-defined, linear and continuous application The next one is our first main theorem. The main assumption is that the component of A associated with the multiindex (0, 0, . . .) is the principal part of A actually governing the equation (assumption (A0)), and that all the subsequent terms of the expansion of A are, at most, operators of order one (assumption (A1)).
Proof Clearly, the hypotheses imply that a(t, u, v) satisfies also (3.4) and (3.5), so that Theorems 3.1 and 3.2, as well as Corollary 3.3, can be applied to (3.2). The solution to (3.2) can be written in the Duhamel form where E 0 (t), E 1 (t), E(t, s) depend only on A (0,0,...) and max{ E 0 (t) , E 1 (t) , E(t, s) } ≤ C by (3.8). Also, by the regularity of the solutions and the fact that all operators A γ with γ > 0 are, at most, of order one, for each δ ∈ I there exists a constant K δ > 0 (actually K δ = A δ L(V ,H ) ) such that A δ u λ (t) H ≤ K δ u λ (t) V , t ∈ I , for all λ ∈ I. Now, by (3.8), for some other constant C > 0, depending only on A (0,0,...) , H and V , that is, Clearly, H ) .
Also, we observe that Thus, we find By the assumptions, Also, by the generalized Hölder inequality, we obtain where, by assumption (3.10), Let us chooseT small enough so that 1 − 2C M 2 AT > 0 and letĨ 0 = [0,T ]. Then, We . Thus, by similar computation as onĨ 0 , we find Note that I = [0, T ] can be covered by finitely many intervals of the formĨ k = [kT , (k +1)T ], say with r intervals. We construct piecewise solutions on each of these intervals, each one continuously extending the previous one, resulting in a solution on the entire [0, T ] interval. Clearly, This establishes that the solution u belongs to L 2 (I , V )⊗(S) −1,− p . After a possible modification on a set of measure zero, similarly as in Corollary 3.3, we obtain that u is in The following corollary is a straightforward consequence of the construction of the solution via the system (3.1). It establishes the stochastic unbiasedness of the solution: The expectation of the solution is the deterministic function which is the unique solution of the (deterministic) PDE obtained by taking expectations of all the stochastic coefficients and stochastic initial data.  In this section we treat linear hyperbolic equations of arbitrary order m ∈ N, and linear hyperbolic systems of first order, whose coefficients are globally defined on the whole space R d .

Linear Hyperbolic Equations with Polynomially Bounded Coefficients
We first consider linear operators of the form The operator L acts onto u(t, where Similarly as in Sect. 3 where we assumed in (A1) of Theorem 3.4 that only the principal part of the operator is of second order and the remaining ones are of first order, we will assume here that only L (0,0,...) is a differential operator of order m, while all other operators L β are of order m − 1 at most. Hence, L β take on the form We also assume f (t, The hyperbolicity of L means that the symbol L m (t, x, τ, ξ) of the SG-principal part L m of L, defined here below, satisfies . . , m. The latter means that, for any α, β ∈ N d 0 , k ∈ N 0 , there exists a constant C jkαβ > 0 such that x, τ, ξ) = 0 with respect to τ are usually called characteristic roots of the operator L.
We will deal with the following three classes of equations of the form (1.2), and corresponding operators L: 1. strictly hyperbolic equations, that is, L m satisfies (4.5) with real-valued, distinct and separated roots τ j , j = 1, . . . , m, in the sense that there exists a constant C > 0 such that (4.6) 2. hyperbolic equations with (roots of) constant multiplicities, that is, L m satisfies (4.5) and the real-valued, characteristic roots can be divided into n groups (1 ≤ n ≤ m) of distinct and separated roots, in the sense that, possibly after a reordering of the τ j , j = 1, . . . , m, there exist l 1 , . . . l n ∈ N with l 1 + . . . + l n = m and n sets satisfying, for a constant C > 0, (4.7) notice that, in the case n = 1, we have only one group of m coinciding roots, that is, L m admits a single real root of multiplicity m, while for n = m we recover the strictly hyperbolic case; the number l = max j=1,...,n l j is the maximum multiplicity of the roots of L m ; 3. hyperbolic equations with involutive roots, that is, L m satisfies (4.5) with realvalued characteristic roots such that

Remark 4.1
Recall that roots of constant multiplicities are always involutive, while the converse statement is not true in general, as shown, e.g., in [8].
Definition 4. 2 We will say that the (linear) operator L in (4.1) and the associated Cauchy problem (1.2) are strictly (SG-)hyperbolic, weakly (SG-)hyperbolic with constant multiplicities, or weakly (SG-)hyperbolic with involutive roots, respectively, if such properties are satisfied by the roots of L m , as explained above.
The next one is a key result in the analysis of SG-hyperbolic Cauchy problems by means of the corresponding class of Fourier operators. Given a symbol and consider the eikonal equation with 0 < T 0 ≤ T . By the theory developed in [1,12], it is possible to prove that the following proposition holds true.

Remark 4.4
Of course, if additional regularity with respect to t ∈ [0, T ] is fulfilled by the symbol κ in the right-hand side of (4.9), this reflects in a corresponding increased regularity of the resulting solution ϕ with respect to (t, s) ∈ T 0 . Since here we are not dealing with problems concerning the t-regularity of the solution, we assume smooth t-dependence of the coefficients of L.
In the approach we follow here, which is the same used in [14] and elsewhere, a further key result is the next proposition, an adapted version of the so-called Mizohata Lemma of Perfect Factorization 2 . formulated for the SG-hyperbolic operator L (0,0,...) in the case of roots with constant multiplicities. Of course, it holds true also in the more restrictive case of strict hyperbolicity, which coincides with the situation where l = max j=1,...,n l j = 1 ⇔ n = m. with  Under the hypotheses of weak SG-hyperbolicity with constant multiplicities or with involutive roots, plus a suitable Levi condition 3 , or of strict SG-hyperbolicity, it is possible to show that the Cauchy problem (4.14) can be solved recursively by induction on the length of the multiindex γ . This follows from the next Theorem 4.6, which summarizes some of the main results proved in [1,4,8,12,14], applied to L (0,0,...) .   associating the right-hand side and the initial data of (4.15) to its solutions and its first m derivatives with respect to t, is a well defined, linear and continuous application uniformly with respect to γ ∈ I.

The entries of the matrix-valued operator T glob can be expressed (modulo smoothing reminders) by means of finite linear combinations 4 of SG-FIOs of the form
Op ϕ j (t) (a j (t)), with smooth regular phase functions families ϕ j (t), obtained as solutions of the eikonal equations (4.9), with the real, distinct characteristic roots θ j of L m in place of κ, and suitable smooth amplitude families a j (t), j = 0, . . . , n, see [1,4,12,14]. The continuity of T glob then follows by Theorem 2.9 and Remark 2.10. The uniform continuity with respect to γ ∈ I is an immediate consequence of the fact that such SG-FIOs, as well as T > 0, only depend on L (0,0,...) . 3. By the continuity of T glob , we obtain an energy-type estimate for the solution u γ and its first m derivatives with respect to t, namely for a suitable constant C > 0, independent of γ ∈ I.
We can now prove the second main result of the paper, which is Theorem 4.8. where l denotes the maximum multiplicity of the distinct characteristic roots of the principal symbol L m in (4.5) (in particular, l = 1 for strictly hyperbolic operators, and l = m for weakly hyperbolic operators with involutive characteristics). In the weakly hyperbolic cases, assume also the corresponding Levi condition, as in Theorem 4.6. We also assume that there exists r ≥ 0 such that (4.18) 4 More precisely, in the case of hyperbolic operators with involutive roots, also compositions of a finite number of operators of the type Op ϕ j (t) (a j (t)), as described in Theorem 2.11, are involved.

Finally, assume, for (s, σ )
Then, there exists a a time-horizon T ∈ (0, T ] such that the Cauchy problem (1.2) Proof By (4.16) in Theorem 4.6, with an argument analogous to the one followed in the proof of Theorem 3.4, we obtain an infinite dimensional system equivalent to (1.2), whose solutions are given by Thus, we find (for some new constantC > 0) By the assumptions, Also, by straightforward estimates, we obtain where, by assumption (4.18), Thus, By possibly further reducing T > 0, so that 1 −C M 2 L T 2 > 0, we find which gives the claim.

Linear Hyperbolic Systems of First Order with Polynomially Bounded Coefficients
We now turn our attention to hyperbolic first order linear systems with coefficients having at most polynomial growth at spatial infinity. Namely, let L now denote the operator . . . , N ) is a parameter-dependent, (N × N )-dimensional, diagonal operator matrix, whose entries λ j (t, x, D x ; ω), j = 1, . . . , N , are pseudodifferential operators with parameter-dependent symbols and R = (R jk ) j,k=1,...,N is a parameter-dependent, (N × N )-dimensional operator matrix of pseudo-differential operators with symbols We consider the Cauchy problem  We proceed as in the previous sections, obtaining a family of Cauchy problems for SG-hyperbolic systems indexed by γ ∈ I, namely,  O(k, k)), s ∈ [0, T ), exists (see [11]), and can, in general, be expressed as a limit of matrices of Fourier integral operators (see [22,49]; see Section 5 of [4] for the SG case). In view of the properties of the family λ j(0,0,...) , j = 1, . . . , N , in the three cases considered here E(t, s) can actually be reduced, modulo smoothing operators, to a finite linear combination of (compositions of) SG Fourier integral operators, see [1,[12][13][14]. The next Theorem 4.10 is the analogue for systems of Theorem 4.6.
Theorem 4.10 Let L be SG-hyperbolic of the form in (4.19), either strictly hyperbolic, weakly hyperbolic with constant multiplicites, or weakly hyperbolic with involutive characteristics, according to Definition 4.9. Assume also that such that, for any k, l ∈ N 0 , ∂ k t ∂ l s E(t, s) belongs to C( T , O(k + l, k + l)).

It follows that, for any
, and, for any M ∈ N, the mapping (W γ , G γ ) → (U γ , D t U γ , . . . , D M t U γ ), associating the right-hand side and the initial data of (4.24) to its solutions and its first M derivatives with respect to t, is a well defined, linear and continuous application uniformly with respect to γ ∈ I. 3. By the continuity of T M glob , we obtain an energy-type estimate for the solution U γ and its first M ∈ N derivatives with respect to t, namely for a suitable constant C > 0, independent of γ ∈ I. 4. Theorem 4.6 is a consequence of Theorem 4.10. Indeed, in the three cases considered here, the Cauchy problem (4.15) is turned, by suitable techniques, into an equivalent Cauchy problem of the form (4.24). The Levi conditions on the lower order terms of the involved operator play a crucial role in this aspect (see, e.g., [1,4,8,12,14,22,[33][34][35][36]) We can now state the third main result of the paper, the next Theorem 4.12. Up to a few minor details, the proof is obtained by the same argument employed to prove Theorem 4.8, and is left for the reader.
Similarly as in Corollary 3.5 we may observe that the solution exhibits the unbiasedness property, i.e. its expectation coincides with the solution of the associated PDE obtained by taking expectations of all stochastic elements in (4.20).

Examples
In this concluding section we provide some examples of hyperbolic SPDEs with singularities and possible applications of the previously obtained results in linear problems that arise in physics, geology, cosmology, engineering and ample other areas of science. Apart from Examples 5.1 and 5.2, which serve as a toy model in order to facilitate our algorithm for solving, the other examples serve to provide a motivation for the study of hyperbolic SPDEs and to illustrate how randomness may occur, but their full solution will be part of a forthcoming paper.
Example 5. 1 We start with the most prominent example, the wave equation as the prototype of hyperbolic PDEs. For technical simplicity, we consider the one-dimensional case where a nice analytic formula is known for the solutions, to serve as a simple illustration for our method. Consider the wave equation with a random wave speed c(t, x, ; ω) which has expectation E(c) = c (0,0,...) > 0 with no time/space dependence (for instance, stationary processes have constant expectations): x), i = 1, 2, involves only first order partial x − derivatives. The solution is hence given by

Remark 5.3
However, it is important to note that assumption (A1) in Theorem 3.4, as well as assumption (4.4) in Theorems 4.8 and 4.12, are sufficient conditions but they are not necessary conditions. If one may obtain other good estimates on the regularity of the solutions of the PDEs defining the propagator of the system, then the highest order x-derivatives may be included as well into the SPDE and the desired convergence of the chaos expansion will follow by similar methods as presented. These alternative estimates are out of the scope of this paper and will be presented elsewhere.

Example 5.4
Other similar examples would include the stochastic Helmholtz equation with a random wave number u tt (t, x; ω)−k(ω)♦u(t, x; ω) = f (t, x; ω), with suitable assumptions on the expectation of k and its coefficients in line with Theorem 3.4; or its multidimensional counterparts with a time-space dependent wave number k(t, x; ω). Note that in case of the purely random (not dependent on time and space variables) wave speed c(ω) and wave number k(ω), the solutions obtained by our chaos expansion method coincide with the solutions obtained in [47], where general equations of the form P(ω, D)♦u(t, x; ω) = f (t, x; ω) were considered.

Example 5.5
In this example we provide some instances of operators with variable (time-space depending) coefficients that provide an insight into the fine difference between strict and weak hyperbolicity. The examples are based on [8].
The following equations serve as a motivation for the further study of hyperbolic SPDEs. Under additional conditions they may be adopted to our setting with (A1) in a similar manner as in the previous examples.

Example 5.6
The elastic wave equation, that accounts both for longitudinal and transverse motion in three dimensions, has the general form ρu tt = (λ + 2μ)∇(∇ · u) − μ∇ × (∇ × u) + f . It describes the propagation of waves in solid elastic materials, e.g. seismic waves in the Earth and ultrasonic waves used to detect flaws in materials [46]. In the previous equation u denotes the displacement vector, f the driving force, ρ the density of the material and λ, μ denote the Lamé parameters related to the elastic properties of the medium describing the strain-stress relations. In most physical models these parameters are constant, but since they are subject to some measuring errors (either due to instrument errors or reading errors), it is more convenient to treat them as random variables. The driving force may also incorporate some randomness, both pure randomness and measuring uncertainty, hence it is modeled as a stochastic process. Hence we arrive at a stochastic hyperbolic equation ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ρ(ω)♦u tt (t, x; ω) = (λ(ω) + 2μ(ω)) ♦∇(∇ · u(t, x; ω)) − μ(ω) ♦∇ × (∇ × u(t, x; ω)) + f (t, x; ω) (t, x; ω) ∈ R × R × u(0, x; ω) = 0, u t (0, x; ω) = 0, (x; ω) ∈ R × , where we assumed zero initial displacement and zero initial velocity. If the material consists of different layers (e.g. the Earth's crust does have this property) then it would be more appropriate to consider the coefficients ρ, λ, μ to be stochastic processes ρ(t, x; ω), λ(t, x; ω), μ(t, x; ω) depending on the layer and on the time evolution of the system. Theorem 4.8 provides the necessary tools to solve the equation in its most general form.

Example 5.7
Other examples where classical hyperbolic PDEs may be replaced by SPDEs with random coefficients involve the telegrapher's equations, where voltage (V ) and current (I ) along a transmission line can be modeled by the wave equation u x x (t, x; ω) = k 1 (ω)♦u tt (t, x; ω) + k 2 (ω)♦u t (t, x; ω) + k 3 (ω). Here u stands either for voltage or current (both follow the same pattern), and k 1 , k 2 , k 3 are random variables depending on conductance, resistance, inductance and capacitance (characteristics of the wire material) that incorporate some randomness due to measuring errors and due to unpredictable inhomogeneity of the material, or they might even be regarded as stochastic processes with time-space dependence.
The telegrapher's equation is formally equivalent to the so called hyperbolic heat conduction equation (relativistic heat conduction equation) sometimes used instead of classical parabolic heat conduction to account for the fact that speed of propagation cannot be infinite and must be bounded by the speed of light in vacuum.
Some other interesting models related to the telegrapher's equation concerning random planar movements of a particle driven by Poissonian forces of the fluid are given in [37].

Example 5.8
The next example is provided in [15] and concerns the study of the internal structure of the sun. The Solar & Heliospheric Observatory (SOHO) project was run by NASA and ESA by launching a spacecraft in 1995 with the mission to measure the motion of the sun's surface. From the pulsating waves around the sun's surface scientists would deduce the location of the origin of the shock waves and gain a certain insight into the inner structure of the sun. Assuming that the shock sources are randomly located on a sphere of radius R inside the sun, the dilatation is governed by the Navier equation given by where B(0, R) is the ball centered at the origin with radius R, c(x; ω) is the speed of wave propagation at position x, ρ(x; ω) is the density at position x (we account for measuring errors by letting c and ρ be random), and f (t, x; ω) models the shock that originates at time t at position x. According to the SOHO website the spacecraft was meant to operate only until 1998, but it was so successful that ESA and NASA decided to prolong its mission and it is sending data obtained from the sun up to this date.
Finally we note the very important fact that in the theory of general relativity Einstein's equations can be converted into a symmetric hyperbolic system of equations [43]. Some other papers also consider as an advantage to apply a stochastic approach and treat diffusions in a Lorentzian manifold using stochastic differential equations in the orthonormal frame bundle of the manifold [9]. Hence the results of Theorem 4.12 may be applied into a newly developing field of general relativity where the space geometry incorporates randomness.