A spectral ansatz for the long-time homogenization of the wave equation

Consider the wave equation with heterogeneous coefficients in the homogenization regime. At large times, the wave interacts in a nontrivial way with the heterogeneities, giving rise to effective dispersive effects. The main achievement of the present work is a new ansatz for the long-time two-scale expansion inspired by spectral analysis. Based on this spectral ansatz, we extend and refine all previous results in the field, proving homogenization up to optimal timescales with optimal error estimates, and covering all the standard assumptions on heterogeneities (both periodic and stationary random settings).

1. Introduction 1.1. General overview. Let d ≥ 1 be the space dimension and let a be a symmetric coefficient field on R d that satisfies the boundedness and ellipticity properties for some λ > 0. We shall consider both the case when a is periodic and the case when a is a stationary ergodic random field (in the latter case, we restrict to a Gaussian model for illustration, cf. Definition 1. 3). Given an impulse f ∈ C ∞ c ((0, ∞); L 2 (R d )), we consider the ancient solution of the associated linear wave equation for t < 0, (1.2) in the homogenization regime 0 < ε ≪ 1, and we are interested in the accurate description of the long-time behavior of the flow. The reason why we focus on ancient solutions (with u ε = f = 0 for t < 0) is to ensure the well-preparedness of the wave and to avoid propagating time oscillations; see the discussion in Section 1.6. On short timescales t = O(1), the standard theory [7] ensures that the flow can be approximated to leading order by the ancient solution of a homogenized wave equation, where the (constant) effective coefficientā is the same as for the homogenization of the corresponding steady-state problem. This means that homogenization and time evolution decouple to leading order on short timescales. As first understood by Santosa and Symes [25], this is however no longer the case on longer timescales: more precisely, a non-trivial interaction between homogenization and time evolution appears as soon as t ≥ O(ε −2 ) in the periodic setting, leading to a dispersive correction to the naïve homogenized wave equation (1.3).
In the periodic setting, the first rigorous analysis of this phenomenon is based on spectral theory, more specifically on Floquet-Bloch theory, and is due to Lamacz [23] in one space dimension, and to Dohmal, Lamacz, and Schweizer [9,10] in higher dimension. They proved the convergence to some suitable dispersive homogenized wave equation up to times t ≪ ε −3 . Due to the use of the Floquet-Bloch theory, it was not clear that this approach could be applied beyond the periodic setting to other standard frameworks for homogenization (such as the case of quasi-periodic or random coefficient fields). To treat such cases, Benoit and the second author developed in [5] an approximate version of the Floquet-Bloch theory, which was inspired by [2] and by the observation that the derivative of the Bloch wave with respect to the wave number at 0 is a multiple of the standard corrector in homogenization. Extending this to all orders, they introduced a notion of 'Taylor-Bloch waves', which approximately diagonalize the elliptic operator −∇ · a∇ at low wavenumber. In contrast to the standard Floquet-Bloch analysis, this approximate spectral approach is easily transposed to the random setting as it does not rely on the existence of exact Bloch waves (see [14] for an extension of these ideas to other regimes). In the periodic case, this allowed them to derive a whole hierarchy of higher-order homogenized equations that are valid to leading order up to times t ≤ O(ε −ℓ ) for any ℓ ≥ 0. These higher-order homogenized equations are essentially well-posed (up to truncating high-frequencies). In the random case, they also managed to cover the case of random coefficient fields, for which a homogenized description can only be found up to some maximal timescale t = O(ε −ℓ * ). Although this analysis allows to reach long timescales, it does not provide approximations with optimal accuracy. There is indeed a strong limitation in the analysis: since the impulse f in equation (1.2) is not adapted to O(ε) oscillations of the coefficients, Bloch waves at low wave number only describe the solution to leading order, and Bloch waves at higher wave number should further be taken into account for a finer description. The main difficulty is that Bloch waves at higher wave number are not easily related to "correctors" in homogenization, so that it was unclear how to improve the accuracy in [5].
Shortly after [5], following a variant of classical two-scale expansion methods [6], Allaire, Lamacz, and Rauch [3] and Abdulle and Pouchon [24,1] obtained similar results in the periodic setting and did improve the accuracy in the description of the wave flow on long timescales in terms of some two-scale expansion. Interestingly, and as opposed to the equations obtained by approximate spectral theory, the homogenized equations obtained in [3,24,1] have to be significantly reformulated several times before they give rise to a well-posed problem (this is called the "criminal approximation" in [3]). Because the number of such reformulations increases with the order of accuracy, and although arbitrarily long times and optimal accuracy can be reached, the "infinite-order" homogenized operator they implicitly define cannot be inverted even for very smooth functions.
On the one hand, approaches inspired by spectral theory are physically-motivated (for waves equations, the spectrum is of the essence), they yield well-posed equations valid for long times, but were so far limited in terms of accuracy. On the other hand, approaches based on systematic two-scale expansions allow to reach both long times and optimal accuracy, but essentially require as many reformulations as the order of accuracy to obtain a well-posed equation, which prevents one from inverting the associated "infinite-order" homogenized operator (in other words, the bound on the homogenization error is not sharp). To sum up, a physically-motivated two-scale expansion to reach long times and optimal accuracy was still missing.
The main aim of this paper is to introduce a full spectral two-scale expansion for (1.2), that extends [5] to any order and allows to invert the "infinite-order" homogenized operator for smooth enough impulse f . This spectral two-scale expansion is defined in Theorem 1, whereas the infinite-order result is given in Corollary 1. This infinite-order result shows that the analysis we do here is indeed paying off (it cannot be obtained using the systematic two-scale expansions of [3,24,1]). The main insight of this work is encapsulated in Proposition 1.5, which reformulates the explicit formula for the solution of (1.2) based on Floquet-Bloch theory in a way that leverages an intrinsic two-scale expansion (which we call spectral two-scale expansion). A fundamental physical feature of this spectral two-scale expansion is the following: corrections due to the fact that the impulse f is not adapted to O(ε) oscillations of the coefficients are local with respect to f , see Remark 1.1. This original feature should be of interest to the engineering community, as it means in particular that the expansion actually reduces (essentially) to the much simpler one used in [5] outside the support of the impulse. The main merit of this work is to work out this spectral two-scale expansion and its combinatorial structure. The adaptation from the periodic to the random setting is essentially routine to the expert in stochastic homogenization. It is however important and shows the limitation of homogenization techniques for waves in random media -which is why a detailed statement is included in Theorem 3.
Last, we also thoroughly discuss the two-scale approach of [3,24,1] (which we call geometric two-scale expansion,cf. Section 1.5), and we relate it to the new spectral twoscale expansion. This essentially amounts to comparing redundant hierarchies of corrector equations, which we do in an algorithmic way. As an output, we improve the error analysis of [3,24,1], cf. (1.12).
The rest of this introduction is organized as follows: In Section 1.2, we state our main results on long-time homogenization, both in the periodic and in the random settings, comparing the results obtained with the spectral and the geometric approaches. Next, in Section 1.3, we comment on the important question of the well-posedness of the formal homogenized equations (cf. (1.5) and (1.10) below). In Sections 1.4 and 1.5, we motivate the special form of the spectral and the hyperbolic two-scale expansions. We conclude in Section 1.6 with a discussion of the well-preparedness assumption for (1.2), going beyond the framework of ancient solutions.
• N stands for the set of nonnegative integers. For a multi-index n = (n 1 , . . . , n k ) ∈ N k , we let |n| = n 1 + . . . + n k . • E [·] stands for expectation in the random setting, and is also used in the periodic setting as a short-hand notation for averaging on the unit cell We denote by C ≥ 1 any constant that only depends on the dimension d and on the ellipticity constant λ in (1.1). We use the notation (resp. ) for ≤ C× (resp. ≥ 1 C ×) up to such a multiplicative constant C. We write ≪ (resp. ≫) for ≤ C× (resp. ≥ C×) up to a sufficiently large multiplicative constant C. We add subscripts to indicate dependence on other parameters. • The ball centered at x of radius r in R d is denoted by B r (x), and we set B(x) = B 1 (x), B r = B r (0), and B = B 1 (0). • When defining hierarchies of correctors and of homogenized coefficients, we take the convention that all quantities that are not defined are implicitly set to zero: e.g. ψ n = 0 for n < 0 andb n = 0 for n ≤ 0 in Definition 2.1, etc.

1.2.
Main results: long-time homogenization. Our main results yield long-time error estimates for the two-scale expansion of the heterogeneous wave equation (1.2) with optimal accuracy up to the optimal maximal timescale. We separately consider the periodic and the random settings; the case of quasi-periodic coefficient fields could be treated as well but is skipped for shortness.
1.2.1. Periodic setting. We start with the case when the coefficient field a is periodic on the unit cell Q = (− 1 2 , 1 2 ) d . The main result of this contribution provides a two-scale expansion with optimal error estimate up to times t ≤ O(ε −ℓ ) for any ℓ ≥ 0. This is obtained by extending the spectral approach of [5] to higher-order accuracy in form of a suitable two-scale expansion: while the error estimate in [5] saturated at O(ε), we now reach accuracy O(ε ℓ t), cf. (1.6). The formal homogenized equation (1.5) takes the form of a dispersive correction of (1.3) and the discussion of its well-posedness is postponed to Section 1.3; note that the homogenized differential operator is necessarily symmetric in the sense thatb k = 0 for all k even. 1 The proof of this main result is displayed in Section 2, together with the definition of spectral correctors. Theorem 1. Let a be Q-periodic. There exist sequences of spectral correctors {ψ n } n and {ζ n,m } n,m obtained as solutions of elliptic problems on the periodic cell Q, there exist a sequence of homogenized tensors {b n } n , and a sequence of Fourier multiplier {γ n } n with |γ n (ξ)| ≤ 1, cf. Definition 2.1, such that the following holds. For all ℓ ≥ 1 and for any impulse f ∈ C ∞ (R; H ∞ (R d )) with f = 0 for t < 0, the ancient solution u ε of the heterogeneous wave equation (1.2) is accurately described by the 'spectral two-scale expansion' whereū ℓ ε is an ancient solution of the following formal homogenized equation on R × R d , in one of the meanings provided in Lemma 2.10, (Here, ψ n is an nth-order tensor field, ζ n,m is an (n + 1)th-order tensor field, andb k is a matrix-valued (k −1)th-order tensor.) More precisely, we have the following error estimate: for all ℓ ≥ 1 and t ≥ 0, . ♦ Remark 1.1. The spectral two-scale expansion (1.4) has an important property of physical interest: the second sum contains a series of terms that are all local with respect to the impulse f , and this local contribution vanishes outside the support of f . This specific form for the expansion owes to the Bloch wave analysis of Section 1.4. This also illustrates the superiority of the present analysis over [5]: although outside the support of the impulse the expansion in [5] has essentially the same form as (1.4) (up to the pseudo-differential operator γ ℓ (ε∇)), it does not reach accuracy beyond O(ε). ♦ Remark 1.2. The above two-scale expansion (1.4) involves the pseudo-differential operator γ ℓ (ε∇). Although convenient in the analysis, it might not be so in practice (e.g. for numerical purposes). As shown in (2.6), we can expand γ ℓ (ε∇) = 1 + ∞ k=2 ε k γ k ℓ ⊙ ∇ k for some explicit tensors {γ k ℓ } k , so that for the purpose of (1.6) it can be approximated to the required accuracy O(ε ℓ ) by a finite-order differential operator. Note that the pseudodifferential operator has some intrinsic spectral interpretation, playing the role of a normalization of Bloch waves, cf. Section 1.4. ♦ The scaling of the above two-scale expansion error estimate (1.6) with respect to the order ℓ is expected to be optimal: it indeed implies the summability of the two-scale expansion (or the invertibility of the associated "infinite-order homogenized operator" as pointed out in the introduction), and it allows to deduce for instance the following corollary.
consider the unique solution of the associated heterogeneous wave equation Then, in terms of the two-scale expansion (1.4), withū ℓ ε now denoting the corresponding solution of the formal homogenized equation (1.5) in the sense of Lemma 2.10, we have for all t ∈ R, In particular, if for instance the impulse takes the form f t (x) = f 1 (t)f 2 (x), where f 1 has a smooth and compactly supported Fourier transform on R and where f 2 has a compactly supported Fourier transform on R d , then the two-scale expansion is summable in the following sense: for 0 < ε ≪ f 1 small enough (only depending on d, λ, f ), for all T ∈ R, For comparison, we display the corresponding result that can be obtained instead of Theorem 1 when using a more standard "geometric" approach to devise a two-scale expansion as in [3,24,1] (see Section 1.5 for an explanation of the naming "geometric"). We slightly improve the error estimates of [3] thanks to the use of suitable flux correctors (see Remark 3.7). Yet, the scaling with respect to ℓ in the error estimate (1.11) is much worse than the one in Theorem 1 by a factor ℓ ℓ , thus showing the advantage of the spectral approach. (In particular, Corollary 1 does not follow from Theorem 2.) The proof is displayed in Section 3, together with the definition of hyperbolic correctors.
Theorem 2. Let a be Q-periodic. There exists a sequence of hyperbolic correctors {φ n,m } n,m obtained as solutions of elliptic problems on the periodic cell Q, and there exists a sequence of homogenized coefficents {ā n,m } n,m , cf. Definition 3.1, such that the following holds. For all ℓ ≥ 1 and for any impulse f ∈ C ∞ (R; H ∞ (R d )) with f = 0 for t < 0, the ancient solution u ε of the heterogeneous wave equation (1.2) is accurately described by the 'hyperbolic two-scale expansion' wherev ℓ ε is an ancient solution of the following formal homogenized equation on R × R d , up to a suitable revamping, cf. (3.7), in one of the meanings provided in Lemma 2.10, (1.10) (Here, φ n,m is an nth-order tensor andā n,m is a matrix-valued (n − 1)th-order tensor.) More precisely, we have the following error estimate: for all ℓ ≥ 1 and t ≥ 0, It is instructive to compare the spectral and geometric two-scale expansions (1.4) and (1.9). Outside the support of f , (1.4) provides an approximation of u ε to order O(ε ℓ ) by a sum of ℓ + 1 terms, whereas (1.9) reaches a similar order of approximation with a sum of 1 4 (ℓ + 1)(ℓ + 4) terms (note that φ n,m vanishes for m odd, cf. Definition 3.1). This difference between O(ℓ) and O(ℓ 2 ) terms in the expansions illustrates the more intrinsic character of the spectral two-scale expansion and its superiority in terms of estimates. There is obviously a link between spectral and hyperbolic correctors, and spectral correctors {ψ n , ζ n,m } n,m can indeed be recovered as linear combinations of hyperbolic correctors {φ n,m } n,m with coefficients that are nonlinear functions of hyperbolic homogenized coefficients {ā n,m } n,m . Working out the precise algorithmic relation between the spectral and geometric approaches is quite involved and necessarily algorithmic. This is the object of Section 4. In particular, in combination with Theorem 1, one can improve (1.11) a posteriori to , (1.12) thus removing the spurious factor ℓ ℓ . This constitutes a significant strengthening of the error analysis of [3] and could not be obtained from the geometric two-scale expansion approach only.
1.2.2. Random setting. We turn to the case of a stationary ergodic random coefficient field a. For simplicity and illustration, we shall focus on the following Gaussian model. Definition 1.3. The coefficient field a is said to be Gaussian with parameter β > 0 if it has the form a = h(G) for some h ∈ C 1 (R k ) d×d with k ≥ 1 and for some R k -valued centered stationary Gaussian random field G such that the covariance function has β-algebraic decay at infinity, |c(x)| ≤ (1 + |x|) −β . ♦ The analysis of Theorem 1 can be repeated in this Gaussian setting, and leads to the following two-scale expansion result with optimal error estimate up to the optimal maximal timescale. Note that the maximal timescale depends on the decay rate β for correlations and saturates in case of integrable decay β > d (which corresponds to the strongest mixing possible in the Gaussian setting and yields the central limit theorem scaling for large-scale averages of the coefficients). This result extends [5] to higher-order accuracy. The proof is displayed in Section 2.8. Exactly as in the periodic case above, a corresponding result could also be obtained in terms of the hyperbolic two-scale expansion; we skip the detail for shortness.
Theorem 3. Let a be Gaussian with parameter β > 0 in the above sense, and define ℓ * := ⌈ β∧d 2 ⌉. We can construct spectral correctors {ψ n } n≤ℓ * and {ζ n,m } n+2m≤ℓ * −3 as well-behaved solutions of some hierarchy of elliptic problems on the probability space, we can construct homogenized tensors {b n } n≤ℓ * , and a Fourier multiplier γ ℓ with |γ ℓ (ξ)| ≤ 1, cf. Appendix A, such that the following holds. For any impulse f ∈ C ∞ (R; H ∞ (R d )) with f = 0 for t < 0, the ancient solution u ε of the heterogeneous wave equation (1.2) is accurately described by the corresponding spectral two-scale expansion S ℓ ε [ū ℓ ε , f ] given in (1.4), whereū ℓ ε is an ancient solution of the corresponding formal homogenized equation (1.5) in one of the meanings provided in Lemma 2.10. More precisely, we have the following error estimate: for all t ≥ 0 and q < ∞, : d even, (1.13) Up to such timescales, the above result can be used in particular to derive ballistic transport properties of the wave, see [5,14], as well as to derive spectral information in a suitable low-energy regime, see [12].

Remarks 1.4 (Extensions).
• Nontrivial initial data: In Theorems 1 and 3, we focus on well-prepared data, or equivalently, on ancient solutions of the heterogeneous wave equation, cf. (1.2). If we rather consider the initial-value problem , then, because of ill-preparedness, an O(ε) contribution with almost-periodic time oscillations with O(ε −1 ) frequency is expected to appear and to maintain forever: this is formally described in Section 1.6 below and shows that a twoscale description cannot hold beyond accuracy O(ε). This issue is naturally by-passed by rather considering the time-averaged solution Indeed, setting t 0 := inf(supp θ), an integration by parts ensures that z ε,θ is the ancient solution of ( 14) in terms of the time-averaged impulse f t θ (x) :=´∞ 0 θ(t − s) f s (x) ds. Hence, an effective approximation for the time-averaged solution z ε,θ is obtained as a direct consequence of Theorems 1 and 3 for ancient solutions.
• Heterogeneous mass density: We may replace ∂ 2 t by ρ( x ε )∂ 2 t in the wave equation (1.2). Provided that the weight function ρ satisfies the uniform non-degeneracy condition 1 C 0 ≤ ρ(x) ≤ C 0 for some constant C 0 > 0, it does not change much in the analysis once the definitions of correctors are suitably adapted. The necessary changes are transparent and we skip the detail for shortness. because of dispersive corrections. Indeed, the next-order homogenized coefficientb 3 can be proved to be non-negative, cf. [25], so that equation (1.5) is ill-posed in general. This is not new to homogenization, and a similar difficulty occurs in the elliptic setting when studying higher-order two-scale expansions, see e.g. [6,15]. In this case, one typically uses an inductive method, which, for the wave equation, would read as follows: for ℓ ≥ 1, set w ℓ ε := ℓ k=1 ε k−1wk , wherew 1 is the solution of ( and where for 2 ≤ k ≤ ℓ we inductively definew k as the unique solution of It is easily checked that thisw ℓ ε indeed satisfies (1.5). However, as originally observed in [3] (in the similar setting of (1.10)), this notion of solution displays an immoderate time growth, which destroys any hope of using it for an accurate description of (1.2) on long timescales. More precisely, the energy norm ∇w ℓ;t ε L 2 (R d ) is expected to behave like O( εt ℓ−1 ), which would make the approximation u ε ∼ S ℓ ε [w ℓ ε , f ] trivially false on long timescales t ≫ ε −1 . This time growth appears as a snowball effect as corrector terms in the above hierarchy of equations for {w k } k≥1 have the preceding profiles as sources.
Instead of this naïve inductive method, one must look for another way to rearrange equation (1.5). In fact, asā ≥ λ Id, we note that the Fourier symbol of the operator (1.15) remains positive in a fixed Fourier support for ε small enough. Therefore, if the spatial Fourier transform of the impulse f is compactly supported (uniformly in time), then for ε small enough the Duhamel formula allows us to define a unique solution that keeps the same compact support in Fourier space at all times. Without this additional assumption on f , the operator (1.15) needs to be modified at high frequencies O( 1 ε ) to ensure ellipticity, and there are different ways to proceed. We briefly describe high-frequency filtering, highorder regularization, and the so-called Boussinesq trick: these three approaches happen to be essentially equivalent up to higher-order O(ε ℓ ) errors, and we refer to Section 2.3 for the details.
(I) High-frequency filtering: In [3], the authors proposed to use a low-pass truncation, which amounts to filtering out high frequencies of the impulse: equation (1.5) is then replaced by for some α ∈ (0, 1) and some χ ∈ C ∞ c (R d ) with χ| 1 2 B = 1 and χ| R d \B = 0. Provided that ε is small enough to ensure that the Fourier symbol of the operator (1.15) is positive on ε −α B, the above equation admits a unique ancient solution that has spatial Fourier transform supported in R × ε −α B. Note that in this approach the choice of χ and α is arbitrary.
(II) Higher-order regularization: The alternative method used in [5] amounts to regularizing the operator (1.15) by adding a high-order positive operator κ ℓ (ε|∇|) ℓ (−△), where the factor κ ℓ > 0 is chosen for instance as the smallest value that ensures the following uniform positivity, With this addition, the obtained operator is uniformly elliptic, so there is a unique ancient solution to the modified equation In this approach, the power ≥ ℓ of the added high-order operator and the constant 1 2 λ ∈ (0, λ) in the uniform positivity condition are again arbitrary. (III) Boussinesq trick: This last method proceeds by rearranging the ill-posed equation (1.5) and is inspired by the standard perturbative procedure to rearrange the ill-posed Boussinesq equation in the theory of water waves, see e.g. [8]. This so-called Boussinesq trick was first adapted to the present setting by Lamacz [23] in one space dimension for ℓ = 3.
It was extended in [9] to higher dimension for ℓ = 3, and further extended to all orders ℓ ≥ 3 by Abdulle and Pouchon [1]. It is somehow of the same spirit as the higher-order regularization above, but with the additional twist that it further uses the wave equation itself. This approach is slightly more intrinsic than the previous two ones, but it also has the disadvantage of involving derivatives of the impulse. Let us illustrate the main idea for ℓ = 3. We first choose κ 3 > 0 as the smallest value such that We then decompose the operator (1.15) as and we use that at leading order the equation (1.5) yields to the effect that one may reformulate (1.5) as up to an error of order O(ε 4 ). By our choice of κ 3 , this equation is well-posed. This method extends to arbitrary order and we refer to Section 2.3 for the details.
Next, we turn to the well-posedness of the corresponding formal homogenized equation (1.10) obtained for the hyperbolic two-scale expansion. While for equation (1.5) the difficulty only came from the lack of ellipticity of the spatial differential operator (1.15), the existence theory for equation (1.10) is much more delicate as this equation involves higherorder mixed space-time derivatives. Just as for equation (1.5), the inductive method (1.16) used in the elliptic setting leads to secular growth of the approximate solution and is of no use in the present situation. Before being able to use high-frequency filtering, high-order regularization, or the Boussinesq trick, we need to get rid of higher-order mixed space-time derivatives in (1.10), which can be done by iteratively using the equation itself (quite in the spirit of the above presentation of the Boussinesq trick for ℓ = 3). This is called the criminal approximation in [3]. The thorough revisiting of this idea is the object of Section 3.3: more precisely, ifv ℓ ε solves (1.10), then it is shown to satisfy an equation of the form where {b n } n coincides with the spectral homogenized coefficients and where {c n } n is some family of nonlinear combinations of the hyperbolic coefficients {ā n,m } n,m ; see Lemma 3.8 for a precise statement. Now the differential operator in the left-hand side of this equation is the same as in (1.5): it displays a lack of ellipticity, but the same approach to wellposedness can be repeated, using either high-frequency filtering, high-order regularization, or the Boussinesq trick. Note that the right-hand side in the above reformulation of the homogenized equation (1.10) differs from the homogenized equation (1.5) obtained with the spectral approach, but there is no contradiction as the spectral and hyperbolic twoscale expansions also differ: this demonstrates the actual complexity of the link between spectral and hyperbolic correctors; see Section 4.

1.4.
Spectral approach and two-scale expansion. This section constitutes the main insight of this contribution: the form of the spectral two-scale ansatz (1.4), For that purpose, we focus on the periodic setting and first proceed to a fine Floquet-Bloch analysis of the solution u ε of the heterogeneous wave equation (1.2). Starting point is an application of the Floquet transform, which is known to transform the heterogeneous wave operator −∇ · a( · ε )∇ on L 2 (R d ) into a family of fibered wave operators on the periodic Bloch space L 2 (Q). As in [5,14,16], it is useful to enrich the structure by considering shifts of the periodic coefficient field and by augmenting the physical space R d to include such shifts: we defineũ ε ∈ L ∞ loc (R + ; L 2 (R d × Q)) such that, for all q ∈ Q,ũ ε (·, q) is the ancient solution of the following shifted wave equation, In this augmented setting, the ε-Floquet transform [14,16] dy, which is Q-periodic with respect to q. The Fourier inversion formula takes on the following guise, cf. [14, Lemma 2.2],w This leads to a direct integral decomposition 2 via Fourier modes e ξ (x) := e iξ·x . Under this decomposition, the above wave equation (1.18) is equivalent to the following family of wave equations on the unit cell Q, for all ξ ∈ R d , where we recall thatf stands for the spatial Fourier transform of f . Solving this equation by means of Duhamel's formula and using (1.19) to invert the Floquet transform, we get with the short-hand notation L ξ := −(∇ + iξ) · a(∇ + iξ). In other words, the solution can be decomposed as a superposition of fibered evolutions on L 2 (Q), and, since the force f is not oscillating, only the fibered spectral measures associated with the constant function 1 do matter. It remains to evaluate the space-time oscillating factor in this formula (1.20) and to extract an effective behavior as ε ≪ 1.
for t < 0, and assume for simplicity that the spatial Fourier transformf is compactly supported uniformly in time. For all ξ, the self-adjoint operator L ξ on L 2 (Q) has discrete spectrum and we denote its smallest eigenvalue by λ ξ ≥ 0. For |ξ| ≪ 1, this eigenvalue is simple and we denote by w ξ a corresponding normalized eigenfunction. We then set is analytic for |ξ| ≪ 1, and that there holds for |ξ| 1 and ε ≪ 1, With this notation, for ε ≪ f 1 (depending on the Fourier support of f ), the above formula (1.20) forũ ε can be expanded as follows: for all n ≥ 1, which for |ξ| 1 is analytic with respect to ε ≪ 1 and satisfies Ψ m We emphasize the structure of the above expansion (1.23). Since in (1.20) only the fibered spectral measures associated with the constant function w 0 ≡ 1 do matter, which is the ground state of L 0 , the main contribution in (1.23) is naturally given by the perturbed ground state w εξ ∝ π εξ 1. However, as the impulse f is not oscillating, hence is not adapted to oscillations of the ground state w εξ , higher modes also create another non-vanishing contribution. In other words: -The first term in (1.23) corresponds to the contribution of the ground state of the fibered operators {L ξ } ξ and the time dependence is expressed by some effective evolution determined by the fibered ground eigenvalues {λ ξ } ξ .
-The second term in (1.23) is only of order O(ε 3 ) and is induced by higher modes. More precisely, the ill-preparedness of the impulse f creates a non-trivial oscillatory contribution in Duhamel's formula due to higher modes, which amounts after time integration to a contribution that is local with respect to f . In particular, note that this term vanishes outside the support of f . For comparison, the Bloch wave approach in [23,9,10,5] rather focused on the first term in (1.23), thus neglecting the O(ε 3 ) contribution of higher modes, and further replaced the oscillating factor π εξ 1 by a simpler (not normalized) proxy that is easier to characterize but that yields an additional O(ε) error.
Proof of Proposition 1.5. In order to evaluate the space-time oscillating factor in (1.20), we must investigate the spectrum of L ξ in the perturbative regime |ξ| ≪ 1. As this selfadjoint operator has compact resolvent by Rellich's theorem, it has discrete spectrum. We denote by λ ξ its smallest eigenvalue, which is nonnegative as L ξ is. Note that for ξ = 0 the smallest eigenvalue of L 0 is λ 0 = 0 and is simple (with constant eigenfunction). Since the perturbation L ξ − L 0 is L 0 -bounded with relative norm 1 + |ξ| 2 , standard perturbation theory [22] together with the discreteness of the spectrum of L 0 ensures that the smallest eigenvalue λ ξ remains simple for |ξ| ≪ 1 small enough. Moreover, the branch of eigenvalues ξ → λ ξ is analytic for |ξ| ≪ 1, and there is a corresponding analytic branch of eigenfunctions. Recall the definition of the corresponding projectors π ξ , π ⊥ ξ in the statement, and note that π 0 1 = 1, so that (1.22) follows. Now expanding the constant function 1 with respect to those projectors, identity (1.20) turns intõ The first right-hand side term is already of the desired form, cf. (1.23). We turn to the second term, which captures the contribution of higher modes. Due to the discreteness of the spectrum, for |ξ| ≪ 1, the operator L ξ has a spectral gap L ξ | ℑπ ⊥ ξ 1. Combined with (1.22), this yields the following bound for (1.24), For n ≥ 1, noting that iterated integrations by parts in the time integral yield for all λ > 0, where we used the vanishing assumption for the impulse at initial time s = 0, it then follows from functional calculus that . Inserting this into (1.25) yields the conclusion (1.23).
We turn to the applicability of this spectral computation beyond the periodic setting. In case of a stationary ergodic random coefficient field a, a similar Floquet decomposition (1.20) can be justified, cf. [14,16], but the corresponding fibered operators {L ξ } ξ are then defined on L 2 (Ω), where Ω is the underlying probability space, and typically have non-discrete spectrum. For instance, the spectrum of L 0 = −∇ · a∇ on L 2 (Ω) is expected to be made of a simple eigenvalue at 0 embedded at the bottom of an absolutely continuous spectrum, cf. [16]. In this setting, the Floquet-Bloch theory fails and the above perturbative spectral computation cannot be adapted. As shown in [5,14], however, an 'approximate spectral theory' can be developed: formal Rayleigh-Schrödinger series can be approximately constructed up to a certain accuracy, leading to approximate Bloch waves that can be used to approximately diagonalize the heterogeneous wave operator and describe the flow on long timescales. Equivalently, we may start from a two-scale ansatz given by the formal ε-expansion of the spectral formula (1.23), and then use PDE techniques to show that it indeed provides a good approximation of the flow to a certain accuracy. This approach is the one we use for the proof of Theorems 1 and 3: it allows us to forget about the underlying spectral interpretation, which is only used to devise an educated guess for a two-scale ansatz that can be used as well in the random setting.
Finally, let us describe the ε-expansion of the spectral formula (1.23), showing that it takes the form of the spectral two-scale ansatz (1.17), and let us derive the relevant hierarchy of PDEs for its coefficients. For that purpose, following [5], we first consider the (not normalized) ground state ψ εξ of L εξ satisfying E [ψ εξ ] = 1. Expanding formally and separating powers of ε in the eigenvalue relation L εξ ψ εξ = λ εξ ψ εξ , we then find that the mapsψ n ξ : Q → C and coefficientsλ n ξ ∈ C are defined inductively byψ 0 ξ = 1,λ 0 ξ = 0, and for all n ≥ 1, ) .
(Recall that by convention we implicitly setψ n ξ ≡ 0 for n < 0.) This hierarchy of equations uniquely defines {ψ n ξ ,λ n ξ } n by the Fredholm alternative since by induction the periodic average of the right-hand side of (1.27) precisely vanishes; cf. Section 2.1. Note that we findλ 1 ξ = 0 in agreement with (1.22). Next, we normalize ψ εξ to get a normalized ground state w εξ and we define the corresponding projections π εξ , cf. (1.21), and, in terms of (1.26), (1.29) We turn to the ε-expansion of {Ψ m ξ,ε } m , cf. (1.24): we write their expansions as and it remains to write PDEs to characterize the coefficients {ζ n,m ξ } n,m . For m = 0, inserting (1.29) and separating powers of ε in the defining equation (1.32) (Recall again that by convention we implicitly setζ n,0 ξ ≡ 0 for n < 0.) Again, these objects are well-defined by induction by the Fredholm alternative since our choice (1.32) precisely ensures that the right-hand side of (1.31) has vanishing periodic average; cf. Section 2.1. Next, we argue iteratively for m ≥ 1: starting from the defining equation with nontrivial integration constant fixed as Putting all this together, using expansions (1.26) and (1.30), and extracting the polynomial ξ-dependence of the coefficients in the notation, the spectral formula (1.23) appears to be precisely equivalent to the spectral two-scale ansatz (1.17), with Fourier multiplier γ given by and withū ε satisfying the associated formal homogenized equation, cf. (1.5), where dispersive corrections correspond to derivatives of the fibered ground states {λ ξ } ξ at ξ = 0. We have also derived the PDE hierarchies defining correctors {ψ n , ζ n,m } n,m , cf. (1.27), (1.31), and (1.33). For the proof of Theorems 1 and 3, we replace infinite series by finite sums and show by PDE techniques that these are still good approximations for the wave flow. In the random setting, just as in the elliptic case [18,4,21,15], we can only solve a finite number of the above corrector equations, which is why a homogenized description is only obtained up to some maximal timescale and accuracy, depending both on space dimension and on mixing properties of the coefficient field.
1.5. Geometric approach and hyperbolic two-scale expansion. Instead of starting from the above spectral analysis, another way to describe oscillations of the solution u ε of the heterogeneous wave equation (1.2) is to appeal more directly to two-scale expansion techniques [6] and rather postulate the following natural hyperbolic two-scale ansatz, . These correctors can be viewed as a refinement of usual elliptic correctors: for all n ≥ 0, the nth-order hyperbolic corrector φ n,0 and homogenized tensorā n,0 defined in Section 3 indeed coincide with their elliptic counterparts [6,21,15]. As in the elliptic setting, hyperbolic correctors have a natural geometric interpretation. We focus on the periodic case to simplify the presentation. The first corrector φ 1,0 , which is the same as in the elliptic setting, is defined to correct Euclidean coordinates ∇ iv is then viewed as an intrinsic Taylor expansion of the limiting profilev in terms of a-harmonic coordinates. In order to describe oscillations with finer accuracy, higher-order correctors are iteratively defined to correct higher-order polynomials and make them adapted to the heterogeneous wave operator. More precisely, higher-order correctors {φ n,m } n,m are defined in such a way that, for any polynomialq in space-time variables x, t, the corrected polynomial captures oscillations of the heterogeneous wave operator in the sense that has no periodic oscillations any longer. In that case, this quantity can automatically be written as for some suitable family {ā n,m } n,m of constant tensors; see Proposition 3.6. This reflects the fact that on large scales the heterogeneous constitutive relation ∇u → a∇u is replaced by the effective relation ∇ū →ā∇ū at leading order, while additional corrective terms must be included when looking for finer accuracy, The difference with the elliptic setting is that in the present hyperbolic setting both space and time variables need to be corrected alike. Comparing to the spectral approach, note that the presence of mixed space-time derivatives in the resulting homogenized equation (1.10) leads to additional well-posedness issues: naïve notions of solution display a secular growth, which was first avoided in [3] as explained at the end of Section 1.3.
1.6. Ill-prepared data. Up to now, we have focussed on well-prepared data, or equivalently, on ancient solutions of the heterogeneous wave equation (1.2). We now briefly discuss the effect of ill-prepared data by means of Floquet-Bloch theory, which provides further insight on the claims of Remark 1.4. For that purpose, as in Proposition 1.5, we assume that the coefficient field a is periodic and that the impulse f has compactly supported spatial Fourier transform: in this setting, for ε small enough, the operator L εξ has discrete spectrum and its lowest eigenvalue λ εξ is simple for all ξ in the Fourier support of f . We then show that, if initial data do not fit spatial oscillations of the ground state, their projection on higher modes propagates and maintains forever, leading to an O(ε) contribution that consists of a superposition of typically incommensurate time oscillations with O(ε −1 ) frequency. This almost-periodic structure prohibits any approximate description by a two-scale expansion beyond accuracy O(ε). More precisely, we consider the initial-value problem By Floquet-Bloch theory, arguing as for (1.25), we now get where the remainder r ε contains all the effects of initial ill-preparedness, Denoting by {ν n εξ } n≥0 the non-decreasing sequence of eigenvalues of L εξ on L 2 (Q) (repeated according to multiplicity, with ν 1 εξ > ν 0 εξ = λ εξ ), and denoting by {γ n εξ } n≥0 a corresponding sequence of normalized eigenfunctions, the above remainder can be written as in terms of κ n εξ := E γ n εξ = O(εξ). As claimed, this shows that the remainder r ε is of order O(ε) and oscillates both in space and time with O(ε −1 ) frequency. Moreover, at a fixed Fourier mode ξ, the time dependence is (typically) almost-periodic, which prohibits any approximate description by means of a two-scale expansion beyond accuracy O(ε).
As explained in Remark 1.4, these complicated oscillations are naturally removed by taking time averages. In terms of the above remainder r ε , this amounts to noting that, so that the above flow decomposition is precisely turned into an expansion of the form (1.23) for the time-averaged flow (1.14).

Spectral approach and two-scale expansion
This section is devoted to the definition of spectral correctors and to the proof of Theorems 1 and 3 and of Corollary 1, including the well-posedness of the formal homogenized equation (1.5).  recalling that we implicitly setψ n ξ ≡ 0 for n < 0 for notational convenience. Factoring out powers of iξ in the above, we may then define the real-valued symmetric tensor fields {ψ n } n and symmetric tensors {b n } n via We shall also use the notationā :=b 1 as the latter coincides with the homogenized coefficient for the associated elliptic equation. Next, for all ℓ ≥ 1, we can define the following Fourier multiplier, • As the choice (2.2) precisely ensures that the right-hand side of (2.1) has vanishing average, an iterative use of the Poincaré inequality on the unit cell Q easily yields the following estimates, for all n, |b n | + ψ n H 1 (Q) ≤ C n . • Since by definition we have E ℓ n=0 ψ n ⊙ (iξ) ⊗n = 1, Jensen's inequality yields E | ℓ n=0 ψ n ⊙ (iξ) ⊗n | 2 ≥ 1, which ensures that the definition (2.4) of the Fourier multiplier γ ℓ (ξ) indeed makes sense and satisfies γ ℓ (ξ) ≤ 1 for all ξ. In addition, in view of (2.5), we can expand, for |ξ| ≪ 1 small enough, for some coefficients {γ k ℓ } k . In addition, this can be truncated to any order n ≥ 0, • The uniform ellipticity of a, cf. (1.1), ensures thatā =b 1 is elliptic, and therefore we haveλ 2 ξ = ξ ·āξ ≥ λ|ξ| 2 for all ξ ∈ R d . As a consequence of the self-adjointness of fibered operators {L ξ } ξ in Section 1.4, their first eigenvalues {λ ξ } ξ are real, hence, in view of (1.26) and (2.3), we deduce thatb n must vanish for all n even. Equivalently,λ n ξ vanishes for n odd. This can alternatively be proven by a direct computation starting from definition (2.2), cf. [5, Proposition 1]; a similar (more involved) argument will be provided in Proposition 3.5, so we skip the detail for now. Proposition 2.3. We haveb n = 0 for all n even. ♦ As is common in the elliptic setting, it is useful to further introduce suitable flux correctors, which will allow to refine error estimates by directly exploiting cancellations due to fluxes having vanishing average.
Next, we turn to the construction of corresponding correctors {ζ n,m } n for m ≥ 1, as motivated in (1.33). The proof of this lemma is analogous to the one above and is skipped for shortness. Again, it is useful to further introduce associated flux correctors, which allow us to refine error estimates by directly exploiting cancellations due to fluxes having vanishing average.

2.2.
Spectral two-scale expansion. Given a smooth functionw, we consider its twoscale expansion of order ℓ ≥ 0 associated with the above-defined spectral correctors, as motivated by spectral considerations in Section 1.4, We show by PDE techniques that this expansion is indeed well-adapted to describe the local behavior of the solution to the hyperbolic equation in the following sense: the heterogeneous hyperbolic operator applied to S ℓ ε [w, f ] is equivalent to a higher-order effective operator applied tow, up to error terms of formal order O(ε ℓ ). The proof is postponed to Section 2.4. Proposition 2.9 (Spectral two-scale expansion). Let ℓ ≥ 1, ε > 0, and letw, f be smooth functions satisfying Then, the associated spectral two-scale expansion of order ℓ, defined in (2.14), satisfies the following relation in the distributional sense in R × R d , However, as explained in the introduction, the symbol of the operator may vanish, so equation (2.16) is ill-posed in general. As described in Section 1.3, several high-order modifications of this equation can then be used to ensure well-posedness: highfrequency filtering as in [3], high-order regularization as in [5], or the Boussinesq trick as in [1]. Let us precisely define each of them: (I) High-frequency filtering. Let α ∈ (0, 1), and let χ ∈ C ∞ c (R d ) be a cut-off function with χ| 1 2 B = 1, χ| R d \B = 0.
(III) Boussinesq trick. Set κ 1 = 1, κ 2j = 0 for all j, and for j > 0 we define inductively κ 2j+1 ≥ 0 as the smallest value such that for all ξ ∈ R d , Note that (2.5) entails |κ l | ≤ C l for all l. Then, for all ε > 0, we can defineū (III),ℓ ε as the unique solution in R × R d of We analyze these three modifications of the formal homogenized equation (2.16) and show that they are well-posed and all equivalent up to higher-order errors. The proof is postponed to Section 2.5. (i) If the spatial Fourier transform of f is supported in R × B R for some R ≥ 1, and provided that εR ≪ 1 is small enough (independently of ℓ), then the formal effective equation (2.16) admits a unique ancient solutionŪ ℓ ε ∈ L ∞ loc (R; L 2 (R d )) with spatial Fourier transform supported in R × B R . Moreover, it satisfies for all r, t ≥ 0, .
and for (⋆) = (III), (iii) If the spatial Fourier transform of f is supported in R × B R for some R ≥ 1, and provided that εR ≪ 1 is small enough (independently of ℓ), we have for all r, t ≥ 0, Proof of Proposition 2.9. By scaling, it suffices to consider the case ε = 1 and we omit the subscript ε = 1 for notational convenience. The two-scale expansion (2.14) can be decomposed as We split the proof into three steps, separately deriving equations for each part. Step A direct calculation yields for all n ≥ 0, Combined with the defining equation (2.1) for ψ n , this entails and thus, after summation over 0 ≤ n ≤ ℓ, Recalling the definition of flux correctors, cf. Definition 2.4, this means ..i ℓ+2w , and the claim (2.23) follows.
Step 2. Equation for S ℓ 2 [f ]: we show that A direct calculation yields for all m, n, For m = 0, inserting the defining equation for ζ n,0 , cf. (2.8), this entails and thus, after summation over 0 ≤ n ≤ ℓ − 3, Proceeding similarly for m > 0, we find Combining the above two identities, we are led to Recalling the definition of flux correctors, cf. Definition 2.7, as well as (2.10) and (2.12), and reorganizing the terms, the claim (2.24) follows. Note that the flux correctors τ n,m 's are nontrivial even for n = −1, but those appear only in the case when ℓ − 3 is odd.
Combining the results of the last two steps, reorganizing the terms, and recalling the relation (2.15) betweenw, f , we are led to (2.25) and it remains to reformulate the second and fourth right-hand side terms. For the second right-hand side term, we recall the definition (2.4) of γ ℓ , which yields For the fourth right-hand side term in (2.25), we use the auxiliary correctors of Definition 2.4 to write for ℓ ≥ 2, Combining these identities yields the conclusion.

2.5.
Proof of Lemma 2.10. We split the proof into five steps.
We turn to the proof of a priori estimates. As the equation is linear and has constant coefficients, derivatives of the solution satisfy the same equation up to replacing f by its corresponding derivatives. It is therefore enough to prove the stated estimates with r = 0, Moreover, we note that the L 2 -estimate (2.29) directly follows from (2.28): recalling that D = (∂ t , ∇), we indeed get from (2.28), by integration, which implies Integrating in time and appealing to (2.26), this yields the claim (2.28). Note that uniqueness follows from these a priori estimates by linearity. For Step 5, we shall also need an a priori estimate forŪ ℓ ε in terms of theḢ −1 -norm of the impulse. Multiplying (2.27) with the complex conjugate of (| · | 2 + δ) −1 ∂ t F[Ū ℓ ε ] and repeating the above argument, we infer that Integrating in time, using the monotone convergence theorem to pass to the limit δ ↓ 0, and appealing again to (2.26), we deduce and similarly, due to the constant coefficients and linearity of the equation, we get for all r ≥ 0, Step 2. Proof of (ii) for high-frequency filtering. Let α, χ be fixed. We appeal to (2.26) with R = ε −α : provided that ε 1−α ≪ 1 is small enough, we deduce for |ξ| ≤ ε −α , Hence, as in Step 1, replacing f by χ(ε α ∇)f , there is a unique solution u (I),ℓ ε of (2.18) in L ∞ loc (R; L 2 (R d )) with spatial Fourier transform supported in R × ε −α B, and the claimed a priori estimates similarly follow.
Step 4. Proof of (ii) for Boussinesq trick. In terms of the symbol Note that (2.32) makes sense as κ l ≥ 0 for all l. In addition, the choice (2.21) of {κ l } l precisely ensures that all the terms of the sum over n in the numerator of (2.32) are nonnegative (and actually vanish for n even due to Proposition 2.3). Only keeping the term for n = 1, and recalling κ 1 = 1 by definition, we deduce the lower bound which is pointwise positive. We can then define a solution of (2.33) via Duhamel's formula and thus, upon inverse Fourier transformation, this provides a weak solution of (2.22) in L ∞ loc (R; L 2 (R d )). We turn to the proof of a priori estimates. As in Step 1, it suffices to establish the estimate in energy norm. For that purpose, we start from the following equivalent formulation of (2.33), in terms of the symbols . Arguing as in Step 1, and using that β ℓ ε (ξ) ≥ 1 and γ ℓ ε (ξ) ≥ λ|ξ| 2 , we find that any ancient solution of (2.34) satisfies Inserting the upper bound we are led to the claimed a priori estimate on the energy norm.
We turn to the case (⋆) = (II). By definition, the differencev (II),ℓ ε satisfies the following equation, Hence, combining (2.30) and the a priori estimate of item (ii), together with the bound κ ℓ ≤ C ℓ , we get It remains to treat the case (⋆) = (III). Starting from (2.16) and (2.22), and recalling κ 1 = 1, we get the following equation for the corresponding difference, and thus, further using the equation forŪ ℓ ε in the right-hand side, we get after reorganizing the terms, Combining (2.30) and the a priori estimate of item (ii), together with the bounds |b k | ≤ C k and κ ℓ ≤ C ℓ , we get ≤ (εC) ℓ t D 2ℓ+r−2 f L 1 ((0,t);L 2 (R d )) , and the conclusion follows.
2.6. Proof of Theorem 1. Let a be Q-periodic. We split the proof into three steps. We first establish (1.6) for the energy norm, before turning to the L 2 -estimate, which requires some additional care. We start by assuming that suppf ⊂ R × B R for some R ≥ 1, and then conclude with the general case in the last step.
Step 1. Proof of (1.6) for the energy norm in case suppf ⊂ R × B R with εR ≪ 1.
For simplicity, we start by assuming momentarily that the corrector estimates in Lemma 2.8 hold uniformly in the sense of and (ζ n,m , τ n,m ) W 1,∞ (Q) ≤ C n+m+1 , for all n, m ≥ 0. (2.37) As suppf ⊂ R × B R with εR ≪ 1, we can consider the solutionŪ ℓ ε of the formal effective equation (2.16) as given by Lemma 2.10(i). From Proposition 2.9 and Lemma B.1, using the assumed uniform corrector estimate (2.37), we then obtain , and thus, combining this with the a priori bounds of Lemma 2.10(i), It remains to replaceŪ ℓ ε byū (⋆),ℓ ε in the left-hand side for (⋆) = (I), (II), or (III). For that purpose, recalling the definition of the spectral two-scale expansion, cf. (2.14), and using the assumed uniform corrector estimates (2.37), we note that , hence, by Lemma 2.10(iii), Combined with (2.38), this proves the claim (1.6) for the energy norm in case suppf ⊂ R × B R with εR ≪ 1, provided that (2.37) holds.
It remains to treat the case when the uniform boundedness assumption (2.37) for correctors is not satisfied. In that case, we rather appeal to the Sobolev embedding to estimate products with correctors: for any periodic corrector or corrector gradient ϕ ∈ {ψ n , ∇ψ n , σ n , ∇σ n , ζ n,m , ∇ζ n,m , τ n,m , ∇τ n,m } n,m , we can estimate, for any function g, Up to a fixed loss of derivatives in the estimates, we may then appeal to the L 2 corrector estimates in Lemma 2.8, and the above proof of (1.6) for the energy norm is adapted directly.
Step 2. Proof of (1.6) for the L 2 -norm in case suppf ⊂ B R with εR ≪ 1. As in Step 1, we aim to apply Proposition 2.9 and Lemma B.1, together with corrector estimates. However, the following terms are a priori problematic in the right-hand side of the equation for the spectral two-scale expansion given by Proposition 2.9, Indeed, these terms are not total derivatives and involveŪ ℓ ε itself: when applying Lemma B.1 to estimate the L 2 -norm of the two-scale expansion error, these terms would therefore con- with a prefactor t 2 instead of t . In order to improve on this, we shall reformulate T ℓ ε as a total time-derivative up to terms that depend only locally on f . Using the short-hand notationL As in the proof of Lemma 2.10(i), cf. (2.26), the assumption εR ≪ 1 precisely ensures that the operatorL ℓ ε : L 2 (R d ) →Ḣ −2 (R d ) can be inverted when restricted to functions with spatial Fourier transform supported in B R . As by definition bothŪ ℓ ε and f have spatial Fourier transform supported in B R , we may then write By uniform ellipticity (2.26), we have for any function g with suppĝ ⊂ B R , (2.42) Now using (2.41) to reformulate T ℓ ε , cf. (2.40), we get Using this to replace the corresponding terms in the equation for the two-scale expansion error in Proposition 2.9, appealing to Lemma B.1 to estimate its L 2 -norm, using the corrector estimates of Lemma 2.8, using the Sobolev embedding to estimate products with correctors as in (2.39), and using (2.42), we get for a > d 2 , and thus, by the a priori estimate of Lemma 2.10(i), It remains to argue as in Step 1 to replaceŪ ℓ ε byū (⋆),ℓ ε in the left-hand side. For that purpose, recalling the definition of the spectral two-scale expansion, cf. (2.14), and using again the corrector estimates of Lemma 2.8 together with the Sobolev embedding as in (2.39), we note that for a > d 2 , hence, by Lemma 2.10(iii), Combined with (2.43), this proves the claim (1.6) for the L 2 -norm in case suppf ⊂ B R with εR ≪ 1.
Step 3. Approximation argument: proof for general f . For R ≥ 1, consider the truncated impulse

and letū
(⋆),ℓ ε,R be the solution of the modified effective equations as given by Lemma 2.10(ii) with impulse f replaced by f R . Recalling the definition of the spectral two-scale expansion, cf. (2.14), and using again corrector estimates, we then note that for a > d 2 , and thus, by linearity and by the a priori estimates of Lemma 2.10(ii), .
By definition of f R , as in (2.35), the right-hand side can now be estimated as follows, for any k ≥ 0, Combining this with the results of Steps 1 and 2, the conclusion follows for instance with the choice R = ε −1/2 and k = 2ℓ.
in favor of which we presently argue. Recall that we assume here where f 1 has a smooth and compactly supported Fourier transform on R and where f 2 has a compactly supported Fourier transform on R d .
The assumption on f 2 entails while the assumption on f 1 yields for s ≤ 0, wheref 1 stands for the temporal Fourier transform of f 1 . The claim (2.44) follows.
2.8. Proof of Theorem 3. We briefly explain how the above proof of Theorem 1 is adapted to the random setting. As for Theorem 1, we may assume suppf ⊂ R × B R with εR ≪ 1, and the general conclusion follows by approximation. As explained in Appendix A, the only difference with respect to the periodic setting is that only a finite number ℓ * = ⌈ β∧d 2 ⌉ of correctors can be defined with stationary gradient, and the highest-order corrector has a nontrivial sublinear growth. Using Proposition 2.9 in combination with Lemma B.1 to estimate the energy norm of the two-scale expansion error of order ℓ ≤ ℓ * , and using the corrector estimates of Appendix A and the Sobolev embedding, we find for a > d 2 , where the weight µ * ℓ originates in the growth of correctors and is defined in (A.1). A novelty with respect to the proof of Theorem 1 is that we now need weighted energy estimates forŪ ℓ ε with sublinear weight µ * ℓ . For that purpose, as the homogenized equation (2.16) has constant coefficients, we note thatŪ ℓ ε displays ballistic transport, and therefore . This is easily obtained by interpolation, using that the weight x corresponds to a derivative in Fourier space and using the Duhamel formula forŪ ℓ ε ; see e.g. [5, Proof of Proposition 3] for the details. The above then becomes . Note that the error estimate for ℓ = ℓ * − 1 is occasionally better than the one for ℓ = ℓ * . Optimizing between the results for ℓ = ℓ * and for ℓ = ℓ * − 1, we easily conclude As in the proof of Theorem 1, we can replaceŪ ℓ * ε in the left-hand side by the solution of any well-posed modification of the formal homogenized equation, and we can derive a similar estimate for the L 2 -norm.

Geometric approach and hyperbolic two-scale expansion
This section is devoted to the definition of hyperbolic correctors and to the proof of Theorem 2, including the well-posedness of the formal homogenized equation (1.10). This essentially constitutes a rewriting of [3,24,1] and is needed to rigorously relate those works to the spectral two-scale expansion, cf. Section 4.  In particular, there holds φ n,m = 0,ā n,m = 0, and q n,m = 0 whenever m is an odd integer. 3 ♦ For all n, we note that φ n,0 coincides with the elliptic corrector of order n, cf. [6]. As is common in the elliptic setting, it is useful to further introduce suitable flux correctors, which indeed allow to refine error estimates and slightly improve on the result of [3]. More precisely, flux correctors are designed to allow a direct optimal exploitation of cancellations due to fluxes having vanishing average E q n,m = 0. We start with the definition of flux correctors for m ≥ 1. In case m = 0, as correctors φ n,0 coincide with elliptic correctors, they display more structure than general hyperbolic correctors. We recall how this structure can be exploited to construct a suitable flux corrector σ n,0 that is skew-symmetric: the following lemma is essentially borrowed from [15]. • We setφ 0,0 := 1 and for n ≥ 1 we defineφ n,0 := (φ n,0 j 1 ...jn ) 1≤j 1 ,...jn,≤d withφ n,0 j 1 ...jn ∈ H 1 per (Q) the periodic scalar field that has vanishing average E φ n,0 j 1 ...jn = 0 and satisfies −∇ · a∇φ n,0 j 1 ...jn = ∇ · aφ n−1,0 j 1 ...j n−1 e jn + e jn ·q n−1,0 j 1 ...j n−1 .

(3.2)
When symmetrizing indices j 1 , . . . , j n , the skew-symmetry of σ n−2,0 allows to drop the corresponding right-hand side term, and we may conclude by induction that φ n,0 coincides with its modified versionφ n,0 up to symmetrization as stated.
An iterative use of the Poincaré inequality on the unit cell Q ensures the well-posedness of the above objects and further provides the following a priori estimates.

3.2.
Hyperbolic two-scale expansion. Given a smooth functionw, we consider its twoscale expansion associated with the above-defined hyperbolic correctors, as in (1.35), and we show that it is well-adapted to describe the local behavior of the solution to the hyperbolic equation in the following sense: as explained in (1.37), the heterogeneous hyperbolic operator applied to H ℓ ε [w] is equivalent to some higher-order effective operator applied tow (up to O(ε ℓ ) error terms). The proof is postponed to Section 3.6. Proposition 3.6 (Hyperbolic two-scale expansion). Let ℓ ≥ 1, ε > 0, and letw, f be smooth function satisfying Then, the associated hyperbolic two-scale expansion of order ℓ, defined in (3.3), satisfies the following relation in the distributional sense in R × R d , Note that the last two right-hand side terms in the above equation for the two-scale expansion H ℓ ε [w] are total derivatives (with respect to time or space): this is not a trivial fact for the last term, as it is based on the possibility of constructing skewsymmetric flux correctors {σ n,0 } n . This happens to be crucial when applying Lemma B.1 in order to deduce an optimal L 2 error estimate. This slightly refines the analysis of [3]. ♦ 3.3. Homogenized equations and secular growth problem. As motivated in Proposition 3.6, cf. (3.4), we consider the following formal homogenized equation, for ℓ ≥ 1, However, this equation mixes higher-order space and time derivatives, and its well-posedness is problematic. To avoid the secular growth problem described in Section 1.3, we follow [3] and first show that the differential operator in (3.5) can be rewritten in such a way that it does no longer mix space and time derivatives. This is achieved by iteratively using the equation to eliminate time derivatives up to higher-order terms, which is referred to as the 'criminal method' in [3]. The proof is postponed to Section 3.7. Note that the new homogenized coefficients {b n } n in this reformulation automatically coincide with the coefficients given by the spectral approach: this can for instance be deduced a posteriori by comparing the corresponding homogenization results, Theorems 1 and 2; this allows us to use already here the same notation {b n } n as for spectral homogenized coefficients. Note that by definitionb 1 =ā 1,0 =ā.
3.5. Proof of Proposition 3.5. We split the proof into two steps.
Inserting this into (3.12) yields the conclusion.
3.7. Proof of Lemma 3.8. We split the proof into three steps.
(ii) Hyperbolic correctors: The correctors {φ n,m , σ n,m } n+m<ℓ * can be uniquely defined by the corrector equations of Section 3.1 as centered stationary random fields, and in addition the correctors φ n,m , σ n,m with n + m = ℓ * can be uniquely defined as (nonstationary) random fields with centered stationary gradient. The homogenized coefficientā n,m is well-defined for n + m ≤ ℓ * . Moreover, for n + m ≤ ℓ * , the following moment bounds hold for all q < ∞ and x ∈ R d , |(φ n,m , σ n,m )| 2 We state the following general a priori estimates for linear wave equations, which are used throughout; a short proof is included for convenience. In contrast with the situation in the elliptic setting, we emphasize that putting the impulse in divergence form essentially only brings an improvement when estimating the L 2 -norm, and not the energy norm.
Lemma B.1 (A priori estimates). Let L be a self-adjoint operator on L 2 (R d ) satisfying the bound −△ ≤ L ≤ −C 0 △ for some constant C 0 < ∞. Given F 1 ∈ L 1 loc (R + ; L 2 (R d )) and F 2 , F 3 ∈ W 1,1 loc (R + ; L 2 (R d )) with F 2 | t=0 = F 3 | t=0 = 0, let z be the solution of the wave equation Then, for all t ≥ 0, we have z t L 2 (R d ) C 0 t F 1 L 1 ((0,t);L 2 (R d )) + (F 2 , F 3 ) L 1 ((0,t);L 2 (R d )) , and Dz t The assumption −△ ≤ L ≤ −C 0 △ entails that √ L defines a bounded linear operator L 2 (R d ) →Ḣ −1 (R d ) with bounded inverse. We may therefore defineF 2 as the solution of √ LF 2 = ∇ · F 2 , which satisfies F 2 L 2 (R d ) C 0 F 2 L 2 (R d ) . The solution z of the wave equation can then be represented in terms of Duhamel's formula, Integrating by parts in the integral for F 3 , with F 3 | t=0 = 0, this can be rewritten as and the stated L 2 -estimate follows from spectral calculus. For the energy estimate, we rather use energy conservation in form of and the conclusion follows.