Effective Models and Numerical Homogenization for Wave Propagation in Heterogeneous Media on Arbitrary Timescales

A family of effective equations for wave propagation in periodic media for arbitrary timescales O(ε-α)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(\varepsilon ^{-\alpha })$$\end{document}, where ε≪1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon \ll 1$$\end{document} is the period of the tensor describing the medium, is proposed. The well-posedness of the effective equations of the family is ensured without requiring a regularization process as in previous models (Benoit and Gloria in Long-time homogenization and asymptotic ballistic transport of classical waves, 2017, arXiv:1701.08600; Allaire et al. in Crime pays; homogenized wave equations for long times, 2018, arXiv:1803.09455). The effective solutions in the family are proved to be ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} close to the original wave in a norm equivalent to the L∞(0,ε-αT;L2(Ω))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm {L}^{\infty }}(0,\varepsilon ^{-\alpha }T;{{\mathrm {L}^{2}}(\varOmega )})$$\end{document} norm. In addition, a numerical procedure for the computation of the effective tensors of arbitrary order is provided. In particular, we present a new relation between the correctors of arbitrary order, which allows to substantially reduce the computational cost of the effective tensors of arbitrary order. This relation is not limited to the effective equations presented in this paper and can be used to compute the effective tensors of alternative effective models.


Introduction
The wave equation in heterogeneous media is widely used in many applications such as seismic inversion, medical imaging or the manufacture of composite materials. We consider the following model problem: let Ω ⊂ R d be a hypercube and let u ε : [0, T ] × Ω → R be the solution of where we require x → u ε (t, x) to be Ω-periodic and the initial conditions u ε (0, x) and ∂ t u ε (0, x) are given. As we allow the domain Ω to be arbitrarily large, (1) can be used to model wave propagation in infinite media. We assume here that the tensor a ε varies at the scale ε 1 while the initial conditions and the source f have wavelength of order O (1). In such multiscale situations, standard numerical methods such as the finite element (FE) method or the finite difference (FD) method are accurate only if the size of the grid resolves the microscopic scale O(ε). Hence, as ε → 0 or as the domain Ω grows the computational cost of the method becomes prohibitive and multiscale numerical methods are needed.
Several multiscale methods for the approximation of (1) are available in the literature. They can be divided into two groups (see [3] for a review). First, the methods suited when the medium does not have scale separation: [26,27,33,34], and [4]. These methods rely on sophisticated finite element spaces relying on the solutions of localized problems at the fine scale. Second, the methods suited when the medium has scale separation (i.e., a special structure of the medium is required). These methods are built in the framework of the heterogeneous multiscale method (HMM): the FD-HMM [11,23] and the FE-HMM [1]. In both methods, the effective behavior of the wave is approximated by solving micro-problems in small sampling domains.
The FD-HMM and the FE-HMM rely on homogenization theory [12,14,17,28,32,36]: they are built to approximate the homogenized equation and thus provide approximations of u ε in an L ∞ (0, T ; L 2 (Ω)) sense. The homogenization of the wave equation (1) is provided in [15]. For a given sequence of tensors {a ε } ε>0 , we have the existence of a subsequence of {u ε } ε>0 that converges weakly * in L ∞ (0, T ; W per (Ω)) to u 0 as ε → 0 (definitions of the functional spaces are provided below). The limit u 0 , called the homogenized solution, solves the homogenized equation with the same initial conditions as for u ε . The homogenized tensor a 0 in (2) is obtained as the G-limit of a subsequence of {a ε } ε>0 (see [19,38]). In general, a 0 depends on the choice of the subsequence and thus no formula is available for its computation. In this paper, we consider periodic media, i.e., we assume that the medium is described by where Y is a reference cell (typically Y = (0, 1) d ). Under assumption (3), a 0 is proved to be constant and an explicit formula is obtained (see, e.g., [12,14,17,28]): it can be computed by means of the first-order correctors, which are defined as the solutions of cell problems (i.e., elliptic equations in Y based on a(y) with periodic boundary conditions). Therefore, in the periodic case the homogenized solution u 0 can be accurately approximated independently of ε. However, for wave propagation on large timescales, u ε develops dispersive effects at the macroscopic scale that are not captured by u 0 . Furthermore, if the initial conditions or the source have high spatial frequencies (in between O (1) and O(ε)), the dispersion appears at shorter times. Hence, to develop numerical homogenization methods for long-time propagation, or in high frequency regimes, new effective models are required.
The study of this dispersion phenomenon has recently been the subject of considerable interest. Analyses for periodic media and timescales O(ε −2 ) are provided in [5,6,8,20,21,29,37] and numerical approaches are studied in [2,10]. A result for locally periodic media for timescales O(ε −2 ) was also obtained in [7]. For arbitrary timescales O(ε −α ), α ∈ N, effective equations were proposed in [9] and [13]. The well-posedness of these equations is obtained using regularization techniques: in [13], the regularization relies on the tuning of an unknown parameter, which poses problems in practice; in [9] a filtering process is introduced (yet not tested in practice).
In this paper, we present two main results first reported in [35,Chap. 5]. The first main result is the definition of a family of effective equations that approximate u ε for arbitrary timescales O(ε −α ). The effective equations, derived by generalizing the technique introduced for timescales O(ε −2 ) in [5], have the form 1 where a 0 is the homogenized tensor and a 2r , b 2r are pairs of nonnegative, symmetric tensors of order 2r + 2 and 2r , respectively, which satisfy constraints based on high order correctors, solutions of cell problems. Note that the correction of the right-hand side generalizes the one introduced in the case α = 2 in [8] and discussed in [5]. For all effective solutionsũ in the family, we prove an error estimate that ensuresũ to be ε close to u ε in the L ∞ (0, ε −α T ; W ) norm (see (5)). In contrast to the effective equations proposed in [9] and [13], the well-posedness of (4) does not rely on regularization but is naturally ensured by the non-negativity of the tensors. The unregularized versions of the effective equations from [9] and [13] do not belong to the family (4) but are closely related: they have the form (4) with a 2r = g 2r , b 2r = 0 and these pairs of tensors satisfy the correct constraints. The issue is however that the sign of g 2r happens to be negative for some r (this is proved for g 2 [18]) which results in the ill-posedness of the unregularized version of these effective equations.
The second main result of the paper is an explicit procedure for the computation of the high order effective tensors {a 2r , b 2r } in (4), for which we provide a new relation between the high order correctors. In particular, while the natural formula to compute a 2r , b 2r requires to solve the cell problems of order 1 to 2r + 1, this relation ensures that only the cell problems of order 1 to r +1 are in fact necessary. The consequence is a significant reduction of the computational cost needed to compute the effective tensors of arbitrary order. We emphasize that this result can also directly be used to reduce the computational cost for the tensors of the effective equations from [9] and [13].
The paper is organized as follows. In Sect. 2, we present our first main result: We derive the family of effective equations and state the error estimate. We then compare the obtained effective equations with the ones from [9] and [13]. In Sect. 3, we construct a numerical procedure to compute the tensors of effective equations. We then present our second main result: a relation between the correctors which allows to reduce the computational cost of the effective tensors. In Sect. 4, we illustrate our theoretical findings in various numerical experiments. Finally, in Sect. 5 we provide the proofs of the main results.

Definitions and Notation
Let us start by introducing some definitions and notations used in the paper. Let H 1 per (Ω) be the closure of the space C ∞ per (Ω) for the H 1 norm. We denote the quotient spaces L 2 (Ω) = L 2 (Ω)/R and W per (Ω) = H 1 per (Ω)/R. The space W per (Ω) (resp. L 2 0 (Ω)) is composed of the zero mean representatives of the equivalence classes in W per (Ω) (resp. L 2 (Ω)). The dual space of W per (Ω) (resp. W per (Ω)) is |Ω| Ω v and (·, ·) Ω denotes the standard inner product in L 2 (Ω). We define the following norm on W per (Ω) Using the Poincaré-Wirtinger inequality, we verify that · W is equivalent to the L 2 norm: We denote Ten n (R d ) the vector space of tensors of order n. In the whole text, we drop the notation of the sum symbol for the dot product between two tensors and use the convention that repeated indices are summed. The subspace of Ten n (R d ) of symmetric tensors is denoted Sym n (R d ), i.e., q ∈ Sym n (R d ) iff q i 1 ···i n = q i σ (1) ···i σ (n) for any permutation of order n σ ∈ S n . We define the symmetrization operator S n : The coordinate S n (q) i 1 ···i n is denoted S n i 1 ···i n {q i 1 ···i n }. We denote = S an equality holding up to symmetries, i.e., for p, q ∈ Ten n (R d ) we have A colon is used to denote the inner product of two tensors in Ten n (R d ), p : q = p i 1 ···i n q i 1 ···i n . We say that a tensor q ∈ Ten 2n (R d ) is major symmetric if it satisfies We say that a tensor q ∈ Ten 2n (R d ) is positive semidefinite if and it is positive definite if the equality in (9) holds only for ξ = 0. The tensor product of p ∈ Ten m (R d ) and q ∈ Ten n (R d ) is the tensor of Ten m+n (R d ) defined as ( p ⊗ q) i 1 ···i m+n = p i 1 ···i m q i m+1 ···i m+n . Note that up to symmetries the tensor product is commutative, i.e., p ⊗ q = S q ⊗ p. We use the shorthand notation The derivative with respect to the i-th space variable x i is denoted ∂ i and the derivation with respect to any other variable is specified. For q ∈ Ten n (R d ), we denote the differential operator q∇ n x :=q i 1 ···i n ∂ n i 1 ···i n .

Settings of the Problem
Recall assumption (3): a ε (x) = a x ε , where y → a(y) is a d × d symmetric, Yperiodic tensor. In addition, we assume that a(y) is uniformly elliptic and bounded, i.e., there exists λ, Λ > 0 such that Without loss of generality, let the reference cell be Y = (0, 1 ) × · · · × (0, d ). We assume that the hypercube Ω = (ω l In particular, (12) ensures that for a Y -periodic function γ , the map x → γ x ε is Ω-periodic (γ is extended to R d by periodicity). Note that the integers n i in (12) can be arbitrarily large. In particular n i can be of order O(ε −α ).

First Main Result: Family of Effective Equations and a Priori Error Estimate
In this section, we present the family of effective equations and provide the corresponding a priori error estimate. In Sect. 2.1, we derive the family in three steps: (i) We discuss the ansatz on the form of the effective equations; (ii) using asymptotic expansion we derive the high order cell problems; (iii) we obtain the constraints on the effective tensors by investigating the well-posedness of the cell problems. In Sect. 2.2, we define rigorously the family of effective equations and state the a priori error estimate. Finally, in Sect. 2.3 we compare the obtained equations with the other effective equations available in the literature. For the sake of readability, we postpone the technical proofs to Sect. 5.

Derivation of the Family of Effective Equations
In the whole derivation, we assume that the data are as regular as necessary. The specific requirements are stated in Theorem 2. Note that we consider here timescales ε −α T with α ≥ 2. For timescales ε −α T with α < 2 it can be shown following similar techniques that the standard homogenized equation is a valid effective model (see [35,Sect. 5.1.1]).

Ansatz on the Form of the Effective Equations
We first discuss the ansatz on the form of the effective equations, which has a major importance in the derivation. We assume that the effective equations have the form where a 0 is the homogenized tensor (26), a 2r ∈ Ten 2r +2 (R d ), b 2r ∈ Ten 2r (R d ) are tensors to be defined and Q is a differential operator to be defined (the construction of Q f is discussed in Remark 9). As discussed in [5], if the set of considered equations is too small, we end up with ill-posed equations. In particular, without the operators −b 2r ∇ 2r x ∂ 2 t in (14), our derivation would lead to the same ill-posed effective equations obtained in [9] and [13] (the unregularized versions).
Following the classical Faedo-Galerkin method (see [24,Chap. 7]), we prove the following well-posedness result for (14). We define the bilinear forms and the associated Banach spaces We call a functionũ ∈ L ∞ (0, Theorem 1 Assume that the tensors a 2r ∈ Ten 2r +2 (R d ), b 2r ∈ Ten 2r (R d ) are positive semidefinite (9) and satisfy the major symmetries (8). Furthermore, assume that the data satisfy the regularity Then, there exists a unique weak solution of (14).
Let us provide a short sketch of the proof. We look for successive approximations of a weak solution in the form u m (t) = m k=0 u m k (t)ϕ k , where {ϕ k } k∈N is a smooth basis of W per (Ω). For each m, u m (t) is obtained as the solution of a well-posed ordinary differential equation. We then prove that the sequence {u m } m≥0 is bounded in L ∞ (0, T ε ; V). In particular, note that the sign assumptions on the tensors ensure We thus obtain the existence of a subsequence that weakly * converges in L ∞ (0, T ε ; V). We can then prove that the weak * limit is the unique weak solution.

Asymptotic Expansion, Inductive Boussinesq Tricks
We make the ansatz that u ε can be approximated by an adaptation ofũ of the form where u k are to be defined and the map y → u k (t, x, y) is Y -periodic. We split the error as and follow the argument presented in [5] based on the error estimate of Lemma 9 (see also [35,Sect. 4.2.2]): for u ε − B εũ L ∞ (0,ε −α T ;W ) to be of order O(ε), we need the terms involvingũ in the remainder to be of order O(ε α+1 ) in the L ∞ (0, ε −α T ; W * per (Ω)) norm (it is sufficient that this holds in the stronger L ∞ (0, ε −α T ; L 2 (Ω)) norm). We now expand r ε : using the equation for u ε (13) and the form of the adaptation (16), we obtain where the operators A yy , A xy , A x x are defined as and the remainder is Classical two-scale asymptotic expansion [14] advises to look for u k of the form where the components of the tensor χ k are Y -periodic functions to be defined. The next step is the main difficulty of the derivation: We must use the effective equation (14) to substitute all the time derivatives in the terms of order O(ε 0 ) to O(ε α ) in (18). Various versions of this process have been used in related works [5,6,8,21,25,30] to obtain well-posed effective equations from ill-posed ones. To refer to this kind of manipulation, the term Boussinesq trick was coined in [8] (see [16] where such tricks are used for several Boussinesq's type of equations).
As these inductive Boussinesq tricks represent a technical challenge, let us explain here the case α = 2 and f = 0 and postpone the general case to Sect. 5.1.1 (Theorem 5 and Lemma 8). For α = 2 and f = 0, the effective equation (14) can be written as We now use the equation to substitute ∂ 2 tũ in the last term and obtain Using (20) and the two last equalities, we find that the time derivatives in (18) can be written as where the remainder R εũ has order O(ε 3 ) in the L ∞ (0, ε −α T ; L 2 (Ω)) norm (forũ sufficiently regular). Using this expression in (18), we obtain the desired development in the case α = 2 and f = 0. This process is generalized in Theorem 5 and Lemma 8. With these results, we are able to rewrite r ε in (18) without time derivative in the terms of order O(ε 0 ) to O(ε α ): defining the tensors c r ∈ Ten 2r +2 (R d ) inductively and p k ∈ Ten k+2 (R d ) as we obtain Provided sufficient regularity ofũ, the remainder R εũ has order O(ε α+1 ) in the L ∞ (0, ε −α T ; L 2 (Ω)) norm. Furthermore, provided sufficient regularity of f , the remainder S ε f has order O(ε) in the L ∞ (0, ε −α T ; L 2 (Ω)) norm (this is sufficient forũ and u ε to be ε close, see Theorem 2). Using the definition of u k in (20) and of A yy , A xy , A x x , we verify that (for k = 0 use u 0 =ũ and χ 0 = 1) where e i denotes the i-th vector of the canonical basis of R d . Hence, canceling successively the terms of order O(ε −1 ) to O(ε α ) in (22), we obtain the cell problems: the correctors {χ k i 1 ··i k } α+2 k=1 are the functions in W per (Y ) such that where the tensors p k are defined in (21) and the symmetrization operator is defined in (6). Note that the cell problems for χ 1 i and χ 2 i j are the same as the cell problems obtained with two-scale asymptotic expansion [14].

Remark 1
In (23), we have chosen symmetric right-hand sides. This is possible thanks to the symmetry of ∇ n xũ in (22), as a term of the form q∇ n xũ can be rewritten as S n (q)∇ n xũ . This choice ensures that the correctors are symmetric tensors functions. In particular, χ k has only k+d−1 k distinct entries instead of d k if it was not symmetric. When it comes to numerical approximation, each distinct entry of the corrector corresponds to a cell problem to solve and this symmetrization saves computational time. In addition, this choice ensures that the odd order cell problems are well-posed unconditionally (see below).

Constraints on the High Order Effective Tensors
The last step in the derivation of the family of effective equations is to obtain the constraints on the effective tensors by imposing the well-posedness of the cell problems. In order to investigate the solvability of (23), let us state the following classical result (obtained with the Fredholm alternative or the Lax-Milgram theorem combined with the characterization of W * per (Y )).

Lemma 1 For an elliptic and bounded tensor a(y), consider the following variational
We now proceed to the two following tasks: First, we verify that the odd order cell problems satisfy unconditionally the solvability condition (25); second, we impose the solvability condition (25) on the right-hand sides of the even order cell problems to obtain the constraints on the effective tensors.
First note that (23a) is well-posed as its right-hand side unconditionally satisfies (25). Next, we verify that the cell problem for χ 2 is well-posed as the homogenized tensor satisfies which ensures the solvability condition (25) to hold. We then continue this process to derive the constraints on the higher order effective tensors imposed by the wellposedness of the higher order cell problems.
We assume that the cell problems are well-posed up to order 2r , r ≥ 1. Consider the cell problem for χ 2r +1 . Recalling that p 2r −1 = 0 and as the correctors χ 1 , . . . , χ 2r have zero mean, we verify that the solvability condition (25) is equivalent to the following relation between the correctors. Lemma 2 For any 1 ≤ r ≤ α/2 , the correctors χ 2r and χ 2r −1 satisfy the equality Remark 2 A proof of Lemma 2 can be found in [35,Lemma 5.2.5]. A similar result is also known in the context of Bloch wave theory (see, e.g., [13,20] and the references therein). As discussed, Lemma 2 guarantees that the odd order cell problems are wellposed unconditionally. Pursuing the reasoning, this result ensures that no operator of odd order is needed in the effective equations. It is also the reason why no additional correction is required in the effective equation when increasing the timescale from an even integer α to α + 1.
Finally, consider the cell problem for χ 2r +2 . Imposing the solvability condition (25) on the right-hand side, we obtain the following constraint on the tensor p 2r : This constraint can be rewritten in terms of the effective tensors a 2r , b 2r using the definition of c r and p 2r in (21): The constraints (28) characterize the family of effective equations. Indeed, if the effective tensors a 2r , b 2r satisfy (28), the cell problems for χ 1 to χ α+2 are well-posed because their right-hand sides satisfy the solvability condition (25). Therefore, the adaptation B εũ given by (16) and (20) is well defined and can be used to prove an error estimate between u ε −ũ. This result is presented in the next section and rigorously proved in Sect. 5.2.

A Priori Error Estimate for the Family of Effective Equations
In the previous section, we derived the constraints on the effective tensors for the effective equations to approximate u ε . We present here an error estimate that ensures that the solutions of the derived effective equations are ε close to u ε in the L ∞ (0, ε −α T ; W ) norm.
Let us first rigorously define the family of effective equations derived in Sect. 2.1.

Definition 1
The family E of effective equations is the set of equations (14), where Q f is defined as and for 1 ≤ r ≤ α/2 the tensors a 2r , b 2r satisfy the following requirements: (i) a 2r ∈ Ten 2r +2 (R d ), b 2r ∈ Ten 2r (R d ); (ii) a 2r and b 2r are positive semidefinite (see (9)); (iii) a 2r and b 2r satisfy the major symmetries (8); (vi) a 2r and b 2r satisfy the constraints (28).

Remark 3
In the case α = 2, the correction Q f of the right-hand side in the effective equations was introduced in [8] and discussed in [5]. We verify that the definition of Q f in (29) leads to a lower constant multiplying f L 1 (0,ε −α T ;H r (α) (Ω)) in estimate (30). More details on the origin of (29) are given in Remarks 9 and 11.
For the effective solutions in family E, we prove in Sect. 5.2 the following error estimate.
Theorem 2 Let d ≤ 3, α ≥ 2 and letũ belong to the family of effective equations E (Definition 1). Furthermore, assume that a(y) ∈ L ∞ (Y ) and that the data andũ satisfy the following regularitỹ where r (α) = α + 2 α/2 + 2. Then, the following error estimate holds where C depends only on T , r =1 , and Y , and the norm · W is defined in (5).

Remark 4
For d ≥ 4, the result of Theorem 2 holds provided a higher regularity on the effective solution and f are assumed. Specifically, assuming that m is a sufficiently large integer for the embedding H m per (Ω) → C 0 (Ω) to hold, the statement of Theorem 2 is true for r (α) = α + 2 α/2 + m.

Remark 5
From the derivation in Sect. 2.1, we may hope that B εũ is a better approximation of u ε thanũ. For example, in the elliptic case the error between u ε and an adaptation can be estimated in the energy norm (see, e.g., [28], note that under (12) there is no boundary layer). However, in the case of the wave equation, the answer is not as simple. In fact, an error estimate in the energy norm can be proved only in some specific settings: for example if f = 0 and for well prepared initial position u 0 of the form B εū 0 for someū 0 . Note that this issue is related to the lack of convergence of the energy of the fine scale wave toward the energy of the homogenized wave (see [15]).

Comparison with Other Effective Equations in the Literature
In [13], effective equations for arbitrary timescales are derived. The settings are slightly different as the wave equation (13) is considered in the whole space R d and with a tensor that can be of a more general nature: periodic, almost periodic, quasiperiodic and random (we refer to [13] for the specific definitions). In these circumstances, [13] proposes effective equations that have the form The effective tensors in this equation are indeed verified to match the tensors g 2r defined in (27) (see [35,Sect. 5.2.6]). Under some weak regularity assumption on the initial condition, an error estimate between u ε andū is proved in the norm. This estimate is a strong theoretical result as it holds in the norm we expect in the context of homogenization of the wave equation. However, the use of the effective equation (31) in practice is problematic as no procedure for the computation of the regularization parameter γ is available. In fact, numerical tests indicate that the range of acceptable values for γ is narrow: if γ is too small, the equation is ill-posed and if γ is too large, the solutionū of (31) does not describe u ε accurately (see [35,Sect. 5.4.3]).
In [9], another effective equation for arbitrary timescales is proposed. The settings of this result are the following: The wave equation (13) is considered in the whole space R d with a periodic tensor, it includes an oscillating density and admits nonzero source and trivial initial conditions. In the particular case of a constant density (as considered in the present paper), the effective equation proposed in [9] reads where S ε 1 , S ε 2 are filtering differential operators that ensure the equation to be wellposed. The main result is an error estimate in energy norm between u ε and an adaptation ofû. In particular, the results enable the adaptation to approximate u ε as accurately as one wants by increasing α accordingly. Nevertheless, the filtering process used in (32) to obtain a well-posed equation has yet to be tested in practice.
Despite the inherent difference between (31) and (32) and effective equations (14) in family E, their respective effective tensors are tied through the relation (see (28)): Hence, if we let b 2r = 0 for all r in the effective equations in E, relation (33) reads g 2r = (−1) r a 2r and we end up with the unregularized versions of (31) and (32) (i.e., γ = 0 and S ε 1 = S ε 2 = Id). Similarly, Eqs. (31) and (32) without their regularization devices satisfy requirement (vi) in Definition 1. However, while the well-posednesses of (31) and (32) rely on their respective regularization process, the well-posedness of the effective equations in E is guaranteed by requirements (ii) and (iii). To fulfill these requirements, we have an explicit and constructive algorithm described in Sect. 3.1.
Finally, before approximating any of the effective equations (14), (31) or (32), the effective tensors g 2r must be computed. For this calculation, a substantial gain of computational time is obtained by using the formula provided by our second main result in Theorem 3, Sect. 3.2.

Second Main Result: Computation of Effective Tensors and Reduction of the Computational Cost
In this section, we provide a numerical procedure for the computation of the tensors of some effective equations in the family E. In particular, in Sect. 3.2 we present a new relation between the correctors that allows to reduce significantly the computational cost for the effective tensors. We emphasize that this result is not limited to the family E and can be used to compute the effective tensors of the alternative effective models from [9] and [13] discussed in Sect. 2.3. The final algorithm is provided in Sect. 3.4.
Our goal is now to construct tensors {a 2r , b 2r } α/2 r =1 satisfying the requirements (i) to (vi). Note that ifq r was positive semidefinite, the pair a 2r =q r and b 2r = 0 would trivially satisfy (i) to (vi). However, although the sign ofq r is unknown for r ≥ 2, it is known thatq 1 is negative definite (see [5,8,18]). Hence, the main challenge of the construction lies in the sign of the tensors and to build valid tensors, we need two basics from the tensor world. First, we use the following result, proved in Sect. 5.3.

Lemma 3
For any n ≥ 1, the tensor S 2n (⊗ n a 0 ) is positive definite.
Second, we use a "matricization" operator, which linearly maps a symmetric tensor q ∈ Sym n (R d ) to a symmetric matrix M(q) whose sign is the same as q, i.e., where ν is a simple tensor-to-vector transformation. One construction for such an operator is provided in Sect. 3.3. With Lemma 3 and the operator M, we are able to build tensors that satisfy (i) to (vi). Lemma 4 For r ≥ 1, assume that a 0 and {a 2s , b 2s } r −1 s=1 have already been computed and letq r be the tensor defined in (35). Define the tensor R ∈ Ten 2r (R d ) as where λ min (·) denotes the minimal eigenvalue of a symmetric matrix and {·} + = max{0, ·}. Then, the tensors satisfy the requirements (i) to (vi) of Definition 1.
Proof First, note that the orders of the tensors in (i) are correct. Second, as a 2r , b 2r are fully symmetric tensors, they trivially satisfy the major symmetries and (iii) is verified. Next, we verify (vi), i.e., a 2r , b 2r satisfy (34) (recall the meaning of = S in (7)): We are left with (ii): verifying that a 2r , b 2r are positive semidefinite. The positive semidefiniteness of b 2r = R follows directly from the non-negativity of δ * and Lemma 3. Finally, let us verify that a 2r is positive semidefinite. Using (36) and the definitions in (37), we have where we denoted v ξ = ν(ξ ) and | · | the Euclidean norm. The definition of δ * in (37) ensures that the right-hand side is nonnegative, proving that a 2r is positive semidefinite.
That concludes the verification of (i) to (vi) and the proof of the lemma is complete.
Lemma 4 allows to compute the tensors of one effective equations in the family E. We emphasize that this is one possible construction among others. To obtain other effective equations, we have many options. For example, note that replacing δ * in (37) with any δ ≥ δ * provides other valid pairs of tensors. Alternatively, we can replace S 2r (⊗ r a 0 ) in the definition of R with another positive definite tensor. Finally, let us note the following alternative (used in [5] in the case α = 2).

Remark 6
Assume that the tensorq r in (35) can be decomposed asq r = Sq r ,1 − q r ,2 ⊗a 0 , whereq r ,2 ∈ Ten 2r (R d ) is positive semidefinite. Then we verify that the pairs (38), where R is defined as define effective equations in the family.

A New Remarkable Relation Between the Correctors to Reduce the Cost of Computation of the Effective Tensors
Let us discuss the cost of the procedure described in the previous section. For an integer α, assume that we want to construct the tensors of an effective equation for a timescale O(ε −α ): a 0 , {a 2r , b 2r } s r =1 , where s = α/2 . The main computational cost of this construction lies in the calculation of the tensors {q r } s r =1 in (35), which involves the tensors {S 2r +2 (g 2r )} s r =1 , where Following this natural-but naive-formula, we thus need approximations of the correctors χ 1 to χ 2s+1 . In this section, we present a new relation between the correctors ensuring that in fact the correctors χ 1 to χ s+1 are sufficient (Theorem 3). Let us quantify the computational gain achieved thanks to this result. As each χ r is a symmetric tensor valued function, it has r +d−1 r distinct components (see Remark 1). The number of cell problems to solve to have all the distinct components of χ 1 to χ k is thus Hence, to compute {q r } s r =1 Theorem 3 allows to avoid the approximation of χ s+2 , . . . , χ 2s+1 , i.e., we spare the approximation of

Remark 7
It is important to note that this gain is not specific to the effective equations presented here. Indeed, formula (41) can directly be used to compute the (sufficient) symmetric part of the high order tensors in the effective models from [9] and [13] discussed in Sect. 2.3.
The new relation between the correctors is presented in the following result, proved in Sect. 5.4. Note that this result was obtained independently in [22,Remark 2.5] in the context of stochastic homogenization.
k=1 be the zero mean correctors defined in (23). Then, for 1 ≤ r ≤ α/2 , the tensor g 2r defined by (40) satisfies the decomposition where the tensors k r , h r ∈ Ten 2r +2 (R d ) are defined as where the second double sum in the definition of h r vanishes for r = 1.
Observe that the tensor h r only depends on the correctors χ 1 to χ r so that decomposition (41) indeed guarantees that S 2r +2 (g 2r ) can be computed from χ 1 to χ r +1 .

Matrix Associated to a Symmetric Tensor of Even Order
In this section, we construct the "matricization" operator used in the construction of Sect. 3.1. This operator maps a given symmetric tensor of even order q to a matrix whose sign is the same as q.
We consider the bilinear map We verify that the cardinality of n) be a bijection. We define then the bijective mapping With these notations, we rewrite the map defined in (43) as Finally, we define the matrix associated to a tensor as (45) We verify that for any ξ, η ∈ Sym n (R d ), qξ : η = M(q)ν(ξ ) · ν(η). Hence, q is positive definite (resp. semidefinite) if and only if M(q) is positive definite (resp. semidefinite). In particular, property (36) holds.

Algorithm for the Computation of the High Order Effective Tensors
We present here the algorithm for the computation of the effective tensors of one equation in the family E. The procedure relies on the construction explained in Lemma 4, the formula provided by Theorem 3, and the "matricization" operator M defined in (45). The set of index I (d, n) is defined in (44).

Algorithm 4
Compute the tensors of an effective equation (14) in the family E.

Numerical Experiments
In this section, we present numerical experiments to illustrate the result of Theorem 2: the effective equations in the family capture the long-time behavior of u ε . Note that reporting numerical error in the approximation of u ε in a pseudo-infinite medium for very large timescales is not conceivable as computing a reference solution is out of reach even for one-dimensional problems. We can however consider a small periodic domain as this setting is covered by our theory. In addition, we will illustrate that high order effective models are also useful when we deal with high frequency regimes.
We consider the one-dimensional model problem (13) given by the data u 0 (x) = e −4x 2 , u 1 = f = 0, a(y) = √ 2 − cos(2π y) (we verify that a 0 = 1) with ε = 1/10, and the periodic domain Ω = (−L, L), L = 84. Let us give some insight on the macroscopic evolution of u ε : the central pulse u 0 separates into left-and right-going waves packets with speed a 0 = ±1; these packets meet at x = L (equivalently x = −L) for t = L + 2k L, k ∈ N, and at x = 0 for t = 2k L, k ∈ N. As time increases, dispersion appears in each packet. As the domain has a periodic boundary, at some point in time the packets start to superpose. In our settings, t = ε −4 is the larger timescale until the bulk of the dispersion spans approximately half of the domain Ω. For α = 0, 2, 4, we compare u ε with the corresponding effective solutions of order s = α/2 , denotedũ α/2 (ũ 0 is the homogenized solution). We approximate u ε using a spectral method grid of size h = ε/16 and a leap frog scheme for the time integration with Δt = h/16. The effective solutions {ũ s } 2 s=0 are approximated with a Fourier method (the effective coefficients are computed using Algorithm 4) on a grid of size h = ε/16 (no time integration). Details on the numerical methods can be found in [35,Sect. 5.3].
The results are displayed in Fig. 1. In the top-left plot, we compare u ε ,ũ 0 ,ũ 1 at t = ε −2 = 10 2 . We observe thatũ 0 does not capture the macroscopic dispersion developed by u ε , whileũ 1 accurately describes it. In the top-right plot, we compare u ε , {ũ s } 2 s=0 , at t = ε −4 = 10 4 (observe that right-and left-going wave packets are in fact superposed). At this timescale, the first-order effective solutionũ 1 does not capture accurately the dispersion anymore, as expected, but the second orderũ 2 does. In the bottom plot, we compare the normalized errors err(ũ s )(t) = (u ε −ũ s )(t) L 2 (Ω) / u ε (t) L 2 (Ω) for the effective solutions {ũ s } 2 s=0 on the time interval [0, ε −4 ] (the x-axis is in log scale). We observe thatũ 0 is accurate up to ε −1 and then deteriorates. The first order u 1 is accurate up to ε −3 , then deteriorates. Finally, the second orderũ 2 has a satisfying accuracy over the whole time interval (0, ε −4 ). The observations of this experiment corroborate the result of Theorem 2: the effective solutionũ α/2 accurately describe u ε up to O(ε −α ) timescales.
In two dimensions, computing a reference solution for the previous experiment has a huge computational cost. However, we can illustrate an interesting fact: high order effective models are useful in high frequency regimes. Let us first give some insights on this fact. We consider the following specific settings: let a x ε be a given tensor with ε > 0 fixed, g(x) = e β 2 |x| 2 be a Gaussian with β = O(1) and ν > 0 a scaling parameter. In (13), we let f = 0 and the initial conditions u 0 (x) = g(νx) and u 1 = 0 (we assume that Ω is arbitrarily large). To specify the dependence of u ε on the parameter ν, we denote it as u ε ν . In the equation for u ε ν , making the changes of variablesx = νx andt = νt and introducingû ε ν (t,x) = u ε ν (t/ν,x/ν), we find and we conclude that u ε ν (t/ν,x/ν) =û ε ν (t,x) = u νε 1 (t,x). In other words, for ν > 1 (i.e., an increase in the frequencies of the initial wave) the long-time effects of u νε 1 can be observed at a shorter time in u ε ν (modulo a contraction of space). However, for high values of ν we meet situations where Theorem 2 does not provide a satisfactory error estimate: on the one hand, in the estimate for u ε ν , the increase of ν deteriorates the error constant in Theorem 2; one the other hand, in the estimate for u νε 1 , the constant is good, but the period of the tensor ε = νε is to far from the asymptotic regime for homogenization to be meaningful. In practice, we observe that the increase of the frequencies of the initial wave (increase of ν) leads to additional dispersive effects in u ε ν (or u νε 1 ). Furthermore, the use of higher order effective models allows to capture these additional effects. To illustrate this, we consider the two dimensional model problem given by the data f = u 1 = 0 and u 0 (x) = g(νx), g(x) = e −20|x| 2 , ν = 5 1/3 , a(y) = 1 − 0.5 cos(2π y 2 ) 0 0 1− 0.5 cos(2π y 2 ) , ε = 1/10.
We compute the effective tensors using Algorithm 4 and approximate u ε ν and {ũ s ν } 3 s=1 at time t = 20. As above, u ε is approximated with a spectral method on a grid of size h = ε/16 and leap frog scheme with Δt = h/100 and {ũ s } 2 s=0 are computed with a Fourier method on a grid of size h = ε/16 (see [35,Sect. 5.3]). In Fig. 2, we compare the solution in the periodic medium u ε ν (top-left) with the effective solutions of order 1,ũ 1 ν (top-right), and order 2,ũ 2 ν (bottom-left). We observe thatũ 2 ν captures more accurately the dispersion developed by u ε ν thanũ 1 ν . This is even better seen in the 1d cut at {x 1 = 0} in the bottom-right plot of Fig. 2. Furthermore, even though distinguishing the higher order dispersion from the ε-scale oscillation is not easy in this regime, in the zoom we can guess that the model of order 3 is better than the order 2.

Proofs of the Main Results
In this section, we provide the proofs of the main results of the paper. In Sect. 5.1, we prove the inductive Boussinesq tricks: Theorem 5 and Lemma 8. In Sect. 5.2, the a priori error estimate for the family of effective equations of Theorem 2 is proved. In Sect. 5.3, we prove the result on positive definite tensors from Lemma 3 and in Sect. 5.4, we prove the new relation between the correctors from Theorem 3.

Inductive Boussinesq Tricks for the Derivation of the Family of Effective Equations
In this section, we present the technical task that was postponed in the derivation of the family of effective equations in Sect. 2.1. Specifically, we provide the result allowing to substitute the time derivatives in the terms of order O(ε 0 ) to O(ε α ) in r ε (18). To that end, the main challenge is to proceed to inductive Boussinesq tricks (Theorem 5). This result is then used in Lemma 8, which provides the specific relation used in Sect. 2.1.

Inductive Boussinesq Tricks
The following theorem is the key result of this section. (14) and let N be an even integer such that 0 ≤ N ≤ 2 α/2 . If the right-hand side of (14) is defined as

Theorem 5 Letũ be the solution of
thenũ satisfies where c r is the tensor defined inductively as and the remainders R ε Nũ and S ε N f are defined in (67) and satisfy the following estimates: for any integer n such that N + n ≤ α (49) where r (α, N , n) = 2 α/2 + N + n and the constants depend only on the tensors a 0 , {a 2r , b 2r } α/2 r =1 .

Remark 9
We verify that the definition (46) of Q f ensures the remainder S ε N f to have maximal order in terms of ε in the second estimate in (49). More specifically, among the right-hand side of the form Q f = α/2 for all even N such that 0 ≤ N ≤ 2 α/2 . This affirmation is proved in Remark 10.
The result of Theorem 5 relies on inductive Boussinesq tricks. Let us summarize this process. We start from the expression of ∂ 2 tũ given by the effective equation (14), which we rewrite as where Q f is defined in (46) (see Remark 9). The result for N = 0 trivially follows (50) with S ε 0 f = G 0 (0) and R ε 0ũ = T 0 (0) (defined in (55) below). Let us then assume that N is an even number such that 2 ≤ N ≤ α. In the right-hand side of (50), we inductively substitute ∂ 2 tũ (using (50) itself) in the terms of order O(ε 2 ) to O(ε N ). To understand this technical process at best, let us define inductively some quantities and functions. Note that these definitions arise in the result of Lemma 5, which exhibits the decomposition process of one Boussinesq trick.
We define the tensors and Associated with the tensors A r ( j) and B r ( j), we define for the given N the functions With these notations, we verify that (50) reads (recall the definition of Q f in (46)) where the remainders are T N (0) = G N (0) = 0 if N = α and otherwise.
The proof is divided into three steps. In the first step (Lemma 5), we apply one Boussinesq trick: in S N ( j), we use (50) to substitute ∂ 2 tũ and decompose the result into terms of interest plus remainders. In the second step (Lemma 6), the decomposition is used inductively to obtain a first version of Theorem 5, which involves the tensor c r instead of c r . The final step is to prove thatc r and c r are equal (Lemma 7).

Lemma 6
The effective solution of (14)ũ satisfies the equality wherec r is the tensor defined as and the remainders R ε Nũ and S ε N f are defined in (67) and satisfy the estimates (49). Lemma 7 For 1 ≤ r ≤ α/2 , the tensorc r , defined in (58), equals the tensor c r , defined in (48).

Remark 10
Let us explain how the definition of Q f is obtained. Assuming that Q f is unknown, (54) reads ∂ tũ = Q f + R N (0) + S N (0) + T N (0). Using inductively (56a) and (56c), which hold independently of Q f , we obtain (57) with S ε N f defined as We want to build Q f so that S ε N f has order O(ε N +2 ) for all even 0 ≤ N ≤ 2 α/2 . Inserting the ansatz Q f = α/2 s=0 (−1) s ε 2s q s ∇ 2s x f , where q s ∈ Ten 2s (R d ) are unknown, we obtain (after some work) Canceling the terms of order O(1) to O(ε N ), we obtain q 0 = 1 and an inductive definition for q m , m ≥ 1. We then verify by induction that q m = b 2m for 1 ≤ m ≤ α/2 , i.e., Q f is defined as (46).

Proof of Lemma 5
To slightly simplify the notation, Let us denoteᾱ = 2 α/2 . We first prove (56a). Using (50), we substitute ∂ 2 tũ in S N ( j) and obtain where the remainder is T N 1 ( j + 1) = 0 if N =ᾱ and otherwise. In order to decompose the terms in (59), we study the two index sets (denoting = N /2) where we recall that 0 ≤ j ≤ − 2. We verify that these sets can be written as Hence, we rewrite (59) as where the remainder is (61) Using the definitions of A r ( j) in (52) and B r ( j) in (51), we verify that the two sums in (60) match R N ( j + 1) and S N ( j + 1), respectively. That proves (56a).
Next, let us prove (56b). Using the definition of Q f in (46) and recalling (53c) and (53d), we havẽ As done above, we decompose the sum and use the expression of J m ( j) to obtaiñ where the remainder is Using the definitions of B m ( j + 1) in (51) and F N ( j + 1) in (53c), we obtain (56b).
Next, we prove (56c). From the definition of S N ( j) in (53b), we use (50) to get where the remainder is The first term of the right-hand side of (63) matches the definition of R N (ᾱ/2) in (53a). That proves (56c). Finally, let us prove (56d). From the definition ofF N ( j), we verify that (56d) holds with G N (N /2) defined as That concludes the proof of Lemma 5.

Proof of Lemma 6
Starting from (54), we use inductively the decompositions in (56) and obtain We define the remainders as  T N ( j) and G N ( j), we verify that R ε Nũ and S ε N f satisfy the estimates in (49). Next, we have to develop N /2 j=0 R N ( j) in (66). Using the definition of R N ( j) in (53a) and exchanging the sums, we find that Using the definition of A r ( j) in (52), we verify that for r = 0, r j=0 A r ( j) = a 0 =c 0 , and for 1 ≤ r ≤ α/2, Combining (66), (67), (68), and (69), we obtain (57) and the proof of Lemma 6 is complete.

Proof of Lemma 7
We prove by induction on r that the tensorc r , defined in (58), equals the tensor c r , defined in (48). The base case is trivially verified asc 0 = a 0 = c 0 . Then, for r ≥ 2, we assume thatc s = c s for 1 ≤ s ≤ r − 1 and we have to verify that the tensor equalsc r , defined in (58). Using the induction assumption and (58), we have Let us denote the triple sum T and its summand x r , j,s = a 2s ⊗B −s ( j − 1)⊗b 2(r − ) . Applying the change of indices m = r − and exchange the sums twice to get Using the definition of B r ( j) in (51), we then have We change the index k = j + 1 in this equality to obtain from (70) (recall that This expression matches the definition ofc r in (58). Hence, we have proved that c r = c r for all 0 ≤ r ≤ α/2 and the proof of Lemma 7 is complete.

Use of the Inductive Boussinesq Tricks
In the asymptotic expansion in Sect. 2.1, we need to substitute the terms involving ∂ 2 tũ in the development of r ε in (18). This task is performed in the following result, obtained thanks to Theorem 5.

Proof of Lemma 8 Let us define
Using Theorem 5 and the definition of p k , we obtain for any choice of even integers 0 ≤ N (k) ≤ α: where z j k and the remainders are defined as with S ε N (k) f and R ε N (k)ũ the remainders defined in (67). The first estimate in (49) implies that Hence, as our goal is for this quantity to have order O(ε α+1 ) (see the discussion in Sect. 2.1), we need to set N (k) such that Dealing with even and odd indices separately, we verify that the smallest even integer satisfying the above inequality is We now rewrite the double sum in the right-hand side of (73) separately for even and odd index k. For the sum on even index k = 2 : 0 ≤ ≤ α/2 , we have where in the first equality we changed the index m = r + and in the second we exchanged the order of summation and used that z 2s+1 k = 0 for any s, k. Similarly, for the sum on odd index k = 2 − 1 : 1 ≤ ≤ α/2 , we have Using the two last equalities in (73), we gather the sums to find Replacing z k− j j by its definition in (74), we obtain (72) and that concludes the proof of Lemma 8.

Proof of the a Priori Error Estimate for the Family of Effective Equations (Theorem 2)
In this section, we prove Theorem 2. Note that the main ingredient of the proof is an adaptation operator based on the adaptation constructed in Sect. 2.1. Let us first introduce some notations. We use a bracket [v] to denote the equivalence class of v ∈ L 2 (Ω) in L 2 (Ω), and a bold face letter v to denote elements of W per (Ω).

Note that the Hilbert space
. We define the following norm on W per (Ω) In particular, we verify that a function w ∈ W per (Ω) satisfies w W = [w] W , where · W is defined in (5).
The proof of Theorem 2 is structured as follows. First, based on the adaptation B εũ defined by (16) and (20), we define the adaptation operator B εũ . Then, using the triangle inequality, we split the error as and estimate the two terms of the right-hand side separately. Let us first discuss the consequences of the assumptions made in the theorem. The fact that the effective equation belongs to the family E (Definition 1) has two major implications: first, the equation is well-posed; second, the tensors {a 2r , b 2r } α/2 r =1 satisfy the constraints (28) and thus the cell problems (23) are well-posed. Hence, we have the existence and uniqueness of the effective solutionũ and of the correctors χ 1 , . . . , χ α+2 . Inductively, we can show that the correctors satisfy the following bound where the constant C depends only on the ellipticity constant λ and the reference cell Y . Let us then investigate the regularity ofũ and f . As we assume d ≤ 3, the embedding H 2 per (Ω) → C 0 per (Ω) is continuous and the regularity assumption ensures that where we recall that r (α) = α + 2 α/2 + 2 (as α ≥ 2 we have r (α) ≥ α + 4).
We are now able to define the adaptation operator. Thanks to the regularity of f and the correctors, S ε f defined in (74) belongs to L 2 (0, ε −α T ; L 2 (Ω)). We can thus define ϕ as the unique solution in L ∞ (0, ε −α T ; W per (Ω)) of where A ε is defined in (17). The adaptation operator is then defined as where B εũ is the adaptation constructed in Sect. 2.1 (see (16) and (20)): Assumption (12) ensures that x) is Ω-periodic. Using the regularity of the correctors andũ, we verify that the adaptation (notice that the regularity of ϕ and its derivatives are weaker than that of B εũ ). (22) together with the cell problems (23) ensure that where R ε iniũ (t) is defined in (19) and R εũ in (74). The following error estimate can then be used to quantify η ε (see [5,Corollary 2.2] or [35,Corollary 4.2.2]).

Lemma 10
Let γ ∈ L 2 per (Y ) and v ∈ H 2 per (Ω). Then, the following estimate holds where the constant C depends only on Y and d.
As Ω satisfies (12), the numbers N i = ω l i −ω r i i ε are integers and the cells constituting Ω belongs to the set {ε(n · + Y ) : Hence, almost every x ∈ Ω can be written as x = ε(ξ + y) for some ξ ∈ Ξ, y ∈ Y . For such triplet (x, ξ, y), the Y -periodic function γ satisfies γ x ε = γ (ξ + y) = γ (y). Let Z ⊂ R d be an open set with a C 1 boundary that contains Y and is contained in its neighborhood, i.e., As Z has a C 1 boundary and d ≤ 3, Sobolev embedding theorem ensures that the embedding H 2 (Z ) → C 0 (Z ) is continuous. Hence, there exists a constant C Y , depending only on Y , such that We now prove (78). Using (79), we have where we made the change of variables x = ε(ξ + y).
where v ξ,ε is the function of C 0 (Ȳ ) defined by v ξ,ε (y) = v ε(ξ + y) . Using (80) gives where we used that every cell ε(ξ + Y ) belongs to the neighborhoods of (6d − 3) cells (including itself). This proves (78) and the proof of the lemma is complete.
With Lemma 10, we can estimate the remainders in (77). Using (74) and estimate (49), we obtain From (75), we verify that α − 1 ≤ N (k) + k ≤ α and thus Similarly, using (49) to estimate R ε ini (19) and S ε f (74), we verify that Remark 11 From (74), we verify that Referring to Remark 9, the definition of Q f ensures that the second term has order O(ε α+1 ). Hence, Q f ensures a lower constant in the second estimate in (83).
With these estimates, we are able to prove Theorem 2.

A Symmetrized Tensor Product of Symmetric Positive Definite Matrices is Positive Definite (Proof of Lemma 3)
Lemma 3 states that the tensor S 2n (⊗ n a 0 ) is positive definite, where a 0 is the homogenized tensor. We prove here that this property is true for any symmetric, positive definite matrix. A first important result is the following.

Lemma 11
Let R ∈ Ten 2n (R d ) be a positive definite tensor and let A ∈ Sym 2 (R d ) be a symmetric, positive definite matrix. Then the tensor of Ten 2n+2 (R d ) defined by A i 1 i 2n+2 R i 2 ··i 2n+1 is positive definite.
Proof As A is symmetric positive definite, the Cholesky factorization provides an invertible matrix H such that A = H T H . For ξ ∈ Sym n+1 (R d ), we thus have (86) As R is positive definite, the equality holds if and only if H r j ξ ji 2 ··i n+1 = 0 for all r , i 2 , . . . , i n+1 ∈ {1, . . . , d}. Let i 2 , · · · , i n+1 be arbitrarily fixed and denote v j = ξ ji 2 ··i n+1 . Hence, we have H r j v j = 0 for all r , which is equivalent to H T v = 0. As H T is regular, we obtain that v = 0. We have proved that the equality in (86) implies ξ = 0. Hence, the tensor is positive definite and the proof of the lemma is complete.
With Lemma 11 at hand, we are able to prove the following result, which implies Lemma 3.

Lemma 12
If A ∈ Sym 2 (R d ) is a symmetric, positive definite matrix, then the tensor S 2n (⊗ n A) ∈ Sym 2n (R d ) is positive definite.
Proof We proceed by induction on n. The case n = 1 is ensured by [5,Lemma 4.1] (it can also be deduced from Lemma 11). We assume that the result holds for 1, . . . , n − 1 and prove it for n. Let ξ ∈ Sym n (R d )\{0}. First, assume that n is odd. Then, the product S 2n (⊗ n A)ξ : ξ is composed of terms of the form i.e., one of the factor A i r i s share on index with both entities of ξ . The induction hypothesis combined with Lemma 11 ensure that all these terms are strictly positive and thus S 2n (⊗ n A) is positive definite. Second, we assume that n is even. Then, the product S 2n (⊗ n A)ξ : ξ is composed of terms of two forms. First, there are terms of the form (87). By the same induction argument as before, they are strictly positive. Second, terms of the form Altogether, we verify that S 2n (⊗ n A)ξ : ξ > 0 and the proof of the lemma is complete.

Proof of the New Relation Between the Correctors (Theorem 3)
We prove the result for r ≥ 1. Note that in the case r = 1, we adopt the convention that empty sums vanish, i.e., 0 k=1 x k = 0. For the sake of clarity, we assume that |Y | = 1 so that (·, 1) Y = · Y .
In Sect. 2.1, we explained that the cell problems (23) are well-posed if and only if p j = S g j . Let us then replace p j by g j in the expression of the cell problems (23). For r ≥ 1, let χ 1 , . . . , χ 2r +1 be the 2r + 1 first zero mean correctors. We define the tensors A k , B k , C k ∈ Ten 2r +2 (R d ) as Note that the symmetry of a ensures the following symmetry relations for A k and C k : Furthermore, using the test function w = χ 2r +2−k i k+1 ··i 2r +2 in the cell problem for χ k i 1 ··i k in (23) and the symmetry of a, we obtain the following relations Define then the tensor T ∈ Ten 2r +2 (R d ) as Using the definition of σ k and (88), we verify that where in the second equality we changed the index k = 2r − m. Using (89), we decompose the tensor T as where the tensors U 1 , U 2 and U 3 are The rest of the proof relies on the following result.

Lemma 13
The tensors U i defined above satisfy the following relations: where h r is the tensor defined in (42).
While the proof of (i) and (ii) is direct, the proof of (iii) requires preliminary work, done in the two following lemmas. We verify that the expression of h r in (42) corresponds precisely to the sum on these two sets of index, i.e., H = S h r .
We now prove Lemma 13.
Finally, we prove (iii). With the definition of σ k , we find where we made the changes of index k = n − 1 and k = 2r + 1 − m, respectively. As (−1) 2r +2−k = (−1) k , using Lemma 14 we obtain As χ 0 = 1 and the correctors have zero mean, we verify that D 0 = 0. Using Lemma 15, we obtain U 3 = −h r and (iii) is proved.
With Lemma 13 at hand we can prove Theorem 3.
This equality matches the decomposition (41) and the proof of Theorem 3 is complete.

Conclusion
In this paper, we presented a family of effective equations for wave propagation in periodic media for arbitrary timescales. In particular, for any given α ≥ 0, our main result (Theorem 2) ensures the effective solutions to be close to u ε in the L ∞ (0, ε −α T ; W ) norm. As emphasized, the effective equations are well-posed without requiring any regularization process. In addition, we described a numerical procedure to compute the effective tensors of equations in the family. We showed that the computational cost of this procedure can be significantly reduced by using a new relation between the correctors. This relation should also be used to compute the tensors of the alternative effective model available in the literature. One question that is raised is how to find the best effective equation in the family? More precisely, can we find a criterion to find an optimal equation in the family and can we build such an equation explicitly? These interrogations also concern the effective models from [9] and [13]. Answers to these questions are exciting topics left for future research.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.