A backscattering model based on corrector theory of homogenization for the random Helmholtz equation

This work concerns the analysis of wave propagation in random media. Our medium of interest is sea ice, which is a composite of a pure ice background and randomly located inclusions of brine and air. From a pulse emitted by a source above the sea ice layer, the main objective of this work is to derive a model for the backscattered signal measured at the source/detector location. The problem is difficult in that, in the practical configuration we consider, the wave impinges on the layer with a non-normal incidence. Since the sea ice is seen by the pulse as an effective (homogenized) medium, the energy is specularly reflected and the backscattered signal vanishes in a first order approximation. What is measured at the detector consists therefore of corrections to leading order terms, and we focus in this work on the homogenization corrector. We describe the propagation by a random Helmholtz equation, and derive an expression of the corrector in this layered framework. We moreover obtain a transport model for quadratic quantities in the random wavefield in a high frequency limit.


Introduction
This work is motivated by the study of electromagnetic wave propagation in sea ice.The latter is a formidably complex multiscale material, as a composite of pure ice, brine pockets, and air inclusions of different sizes, shapes, and contrasts.Sea ice is often represented as a layered medium, where the volume fraction and the nature of the inclusions vary from layer to layer.Sea ice is dispersive, mostly because of the high salt content of brine, and anisotropic due to the particular elongated shape of the brine pockets.Its physical properties, e.g.complex permittivity, depend on various parameters such as the temperature and the salinity, and there is quite a large literature on experimental estimations of such quantities, see e.g. the monograph [25].From a theoretical viewpoint, sea ice is modeled by a background (the pure ice), with randomly located inclusions of brine and air with shapes following some appropriate statistical distributions.There is usually a distinction between the young ice that is a few meters thick, and the multi-year ice that can be up to ten meters thick [17].
The problem we are interested in here is the estimation, from radar data, of the thickness of the sea ice layer together with its effective permittivity.The particular experimental configuration we consider is that of a plane, equipped with a radar system, flying at a given height and along a given path.The main question is then to determine whether the layer is sufficiently thick for the plane to land.The essential difficulty in doing so is that the beam emitted by the radar system is not perpendicular to the sea ice layer, but impinges on it with a given angle (the reason for this is to remove the so-called left-right ambiguity in some radar systems).
This difficulty is explained as follows.In the frequency range the radar operates, say from 100MHz to 1GHz (higher frequencies do not penetrate well in the sea ice, see e.g.[9], and are therefore not adequate for our purpose), there is a clear separation of scales between the wavelength and the size of the inclusions: typical numbers given in [26] are 5mm for the air bubbles radius, and 5mm for the major axis of brines pockets with 0.025mm for the transversal axis; since the wavelength ranges from about 30cm to a few meters, it is then always significantly larger than the inclusions.This corresponds to a homogenization regime, where wave propagation in a highly heterogeneous medium is accurately approximated by wave propagation in a constant, effective medium.Note that even though the volume fraction of the brine pockets is not too large (10% for brine pockets, 25% for air bubbles in some appropriate conditions [26]), their effect is not negligible due to their very large contrast: while the dielectric constant of pure ice is about 3.14 with essentially no absorption, that of brine is about 70 with an absorption of 20 [26].
If the plane is at an altitude of 50 to 100 meters, we can assume that waves propagating in the air are in a high frequency regime, and therefore that ray theory applies.The consequence of this is, within the homogenization/ray theory picture and assuming the air/sea ice interface is flat, that there should simply be no backscattered signal at the antenna as an application of Snell-Descartes laws.This is of course not the case in practice, and means that what is measured at the detector corresponds to correctors to these approximations.Deriving models for these measurements is the main objective of this work.
There are various sources of corrections to the homogenization/ray theory description: (i) surface roughness at the air/sea ice interface (ii) surface roughness at the sea ice/seawater interface (iii) surface roughness between layers in the sea ice (iv) correctors to the geometrical optics approximation (v) correctors to the homogenization limit.Besides, a complete description of the electromagnetic wave propagation in the sea ice involves the resolution of Maxwell's equations, in particular in order to account for anisotropy and polarization effects.Our ultimate goal is then to treat (i)-(v) in Maxwell's picture to model the backscattered signal.This is an extremely challenging task to account for all these phenomena at once, all the more in the context of Maxwell's equations.We decided therefore to start our program by designing a simplified, yet realistic, model for wave propagation that retains the essential feature of the true physical situation, namely that the measured data are essentially obtained from corrections to the homogenization/ray theory approximation.
We will then describe the propagation by a scalar wave equation, namely the Helmholtz equation, ignoring anisotropy and polarization.We will also suppose that we are in a regime where the effective (homogenized) coefficients of sea ice are simply given by the expectation of the random coefficients.This allows us to obtain a relatively direct characterization of the corrector to homogenization as a Gaussian field.We will provide analytical conditions for this assumption to hold, and also identify experimental situations where the assumption seems reasonable.Note that when the homogenized coefficients are not given by the expectation, the characterization of the corrector is considerably more difficult.There has been a lot of progress in this direction over the last years, mostly in the context of elliptic (diffusion) equations, see [1,14,13] for a few references.One of our future research plans is then to adapt these new results to our present problem.Finally, the characterization of the corrector in the stochastic homogenization of Maxwell's equations is an even more difficult question; to the best of our knowledge a precise characterization of the corrector in that setting remains open to this day.
Two types of asymptotic analysis will be performed in this work: first, a fluctuation theory to characterize the limit of the corrector to the homogenized model; then a high frequency limit via Wigner transforms to obtain a simplified transport model for correlations of the wavefield.We will provide conditions for these two limits to hold at once.In the mathematical analysis, we will focus on the corrector to homogenization aspects by assuming that the sea ice consists of only one layer with flat interfaces.We will then generalize our results without proofs to a multi-layer non-flat interfaces situation assuming there is a separation of scales between the sea ice inhomogeneities and the rough boundaries.Our proofs are based on adaptations of fairly standard methods, and what seems original in this work is their combination and the application to the sea ice problem.
The paper is structured as follows: In section 2, we set up the model, define the scalings and the random field modeling the heterogeneities.In section 3, we state our main results: the characterization of the corrector to homogenization is given in Theorems 3.1 and 3.2, and a transport model for the correlations of the corrector in Theorem 3.3.In section 3.3, we provide generalizations to more complex settings of multi-layer ices and non-flat interfaces.A conclusion is offered in section 5.The other sections of the article are dedicated to the proofs of our main results, in particular to estimates on the Green's function of the homogenized problem.
Acknowledgments.WJ's work is partially supported by the NSF of China under Grant No. 11701314.OP's work is supported by NSF CAREER Grant DMS-1452349.OP is indebted to Professor Margaret Cheney for bringing up the problem addressed in this work and for many stimulating discussions.

The Mathematical Model
Our starting point is the Helmholtz equation in R d , with d = 2, 3: The function c(x) is the speed, and S ω is the source with (angular) frequency ω.A point x ∈ R d will be written as x = (x ′ , z), where z ∈ R and x ′ ∈ R d−1 .We assume that the medium of propagation consists of three layers separated by two horizontal flat interfaces: the upper half-space (the air) {(x ′ , z) ∈ R d : z > 0}, the sea ice {−L 0 < z < 0}, and the seawater {z < −L 0 } below the sea ice.The speed function c(x) hence has the form 1) needs to be supplemented with appropriate conditions at the infinity in order to have a unique solution.In the layered (waveguide) structure we consider here, it is not fully straightforward to adapt the classical Sommerfeld radiation conditions.The matter is discussed in [8] in two dimensions, and the correct conditions are found by studying separately the different kind of modes supported by the problem (e.g.evanescent or radiative).In order to simplify the mathematical analysis, we bypass the question of the boundary conditions and add some (small) artificial absorption to (2.1), which immediately yields a unique solution in L 2 (R d ) provided that S ω ∈ L 2 (R d ).This has essentially no effects on our results since our estimates are made independent of the artificial absorption.The equivalence between imposing Sommerfeld conditions and adding small absorption is not trivial to establish and is investigated in [6] for instance.
When the radar system is carried by a plane flying at height H, with velocity V p along a straight line in the direction of x 1 axis, S ω has the form (for d = 3), Since V p is usually much smaller than c 0 , the source is considered as fixed.Letting t vary yields different sets of sources and measurements, which will be useful when calculating statistical averages provided that the underlying random medium is stationary and ergodic.Above, g models the frequency distribution of the source with central frequency ω 0 and bandwidth B.
The function χ 0 is typically the characteristic function of a rectangle appropriately oriented (this models a beam with a direction orthogonal to the rectangle).
The simple model we just set up can be readily generalized as follows: firstly, while we considered above only one type of heterogeneities at the scale ℓ c (for instance the brine pockets), additional (independent) random fields could be added to account for more inclusions (e.g. the air bubbles); secondly, a multi-layer model in the sea ice could be considered, to separate for instance young ice to multi-year ice.We will perform the mathematical analysis in the simplest case, and state the generalized results without proofs since the augmented models do not require further essential modifications except for more involved algebra.Finally, we remark on the situation when the interfaces are not flat and vary at a scale larger than that of media heterogeneities.From this separation of the scales, we expect that our results hold with just a rewriting of the interfaces.Nevertheless, our method of proof would have to be adapted since it relies on the translational invariance in the transverse variables, which would no longer hold.
The next section is devoted to the comparison between various length scales of the problem and to the discussion of appropriate scalings.

Scalings
The largest scale in the problem is the height of plane H (tenths of meters), and the smallest one is the typical size ℓ c of the heterogeneities; we recall that the latter is around the millimeter for the air bubbles and for the long axis of the brine pockets.The central wavelength is In the regime we consider here, the central frequency is between 100MHz and 1Ghz, which corresponds to wavelengths between roughly 3m and 30cm.As a consequence, Since the wavelength in pure ice is typically between λ 0 /2 and λ 0 , this means that the length scale of heterogeneity is much smaller than the wavelength, and this is the typical homogenization regime: the wavefield in the heterogeneous sea ice can be well approximated by waves in a homogeneous effective medium.Deriving the homogenized equation and the first-order corrector is the first step of our analysis.We assume that the above relation between the length scales holds for all frequencies in the bandwidth of the source.In other words, we impose that, with ω In terms of B, this means B ≪ ω 0 .
In the second step of the analysis, we derive transport equations for the signals scattered by the sea ice.They are the ones associated with the corrector term.We will see that the latter is a Gaussian field, and, as a consequence, the information is contained in averages of quadratic quantities of the field which can be asymptotically described by transport equations.This corresponds to a high frequency regime.With H of order of tenths of meters, and the central wavelength λ 0 between 3m and 30cm, we clearly have H ≫ λ 0 and the high frequency assumption holds in the air.In the sea ice layer, of thickness between 5m and 10m (these are numbers associated with multi-year ice, younger ice is typically less than 2m thick and does not seem able to withstand the weight of a plane and is therefore not considered), the wavelength in the effective medium is around λ e = λ 0 /2 for a source with range 100Mhz-1Ghz.Hence, for frequencies larger than 300Mhz, the ratio λ e /L 0 is less than 0.1 (this is the worst case), and we can assume as well that the high frequency assumption holds.We then define the following non-dimensional parameters and we set β = εη ≪ 1.
We rescale the spatial variable as x = x ′ H (where x ′ is the new variable though, in the presentation, we drop the primes).With k = 2πH/λ 0 = 2π/η ≫ 1 the rescaled wavenumber, L = L 0 /H, and f (x) = H 2 S ω (Hx), (2.1) becomes where the complex refractive index n β (x) reads Above, α ≪ 1 is the artificial absorption, n 0 = 1, n 1 (resp.n 2 ) is the real part of the index of pure ice (resp.seawater), and κ 1 and κ 2 the absorption in pure ice and in seawater.Note that (n j , κ j ), j = {1, 2}, depends on the frequency even though we will not make this explicit.We chose the above representation of n β for simplicity as it allows us to factor out the term k 2 .The absorptions will be rescaled by a factor k later on.The parameter L is of order 0.1 if H ∼ 50m and L 0 ∼ 5m, which we consider to be large compared to η of order 0.01 for λ 0 ∼ 50cm, and to ε of order 0.005 for ℓ c ∼ 1mm.We will therefore keep L fixed in the expansions.
Some notations.We denote by S = {x = (x ′ , z) ∈ R d : z ∈ (−L, 0)} the (rescaled) sea ice layer, and will use the principal square root of a complex number, defined by, for a = u + iv with v = 0, and denote by q β the rescaled field q(•/β), with β defined before, and q r and q i the real and imaginary parts of q.Above, σ is real and σ 2 is the variance of V (and therefore q has variance one), which does not depend on x by the stationarity hypothesis below.We make the following assumptions on V : (A1) V is stationary with mean n 2 V + iκ V , with κ V > 0. This means that there exists a random variable Ṽ and a family of measure-preserving translations {τ x } x∈R d , such that E{ Ṽ } = n 2 V + iκ V and V (x, θ) = Ṽ (τ x θ).(A2) There exists C 0 > 0 deterministic, such that q L ∞ x ≤ C 0 a.e. in Ω. (A3) There exists n m and κ m such that ℜ{n 2 β (x)} ≥ n 2 m > 0 and ℑ{n 2 β (x)} ≥ κ m > 0, for all x ∈ S.
The boundedness assumption (A2) is made to avoid technicalities and can be weakened; assumption (A3) is there to ensure that n 2 β (x) makes sense from a physical standpoint for all realizations.For the model (2.3), the effective (i.e.homogenized) coefficient is simply given by the expectation of the random coefficient, that is, using the stationarity of V , where n 2 e = n 2 1 +n 2 V and κ e = κ 1 +κ V .Later in this section, we will compare certain properties of this homogenized model with experimental data in the physics literature, and provide a more quantitative criterion for the validity of a model in which the homogenized coefficients are just the averages of the random coefficients.
We further assume that the random field q gets decorrelated fast enough.To describe this, we introduce the maximal correlation function ̺ : (0, ∞) → [0, 1], defined to be where C denotes the space of compact sets in R d , dist(U, V ) = min x∈U,y∈V |x−y| is the distance from U to V , F U denotes the sub-σ-algebra of F generated by q r | U and q i | U , and L 2 0 (F U ) is the space of square integrable mean zero real random variables measurable with respect to F U .We assume that: (A4) the maximal correlation function satisfies In addition, we will need the following matrix M : where Note that M is nonnegative definite by Bochner's theorem, and that the latter integrals are all finite as a consequence of (A4), as for instance, The first two integrals are positive according to Bochner's theorem.
A concrete model for V .A natural way to model the distribution of the heterogeneities is by a collection of spheres with radii {R i } centered at {x i } and with (complex-valued) contrasts {τ i }.We then consider for h the characteristic function of the unit sphere.Note that the elongated shape of the brine pockets could easily be implemented in this framework by setting for instance where R i is the length of the pockets along the principal axis, γ < 1 models the aspect ratio, and the ˜in the coordinates represents a (random) rotation of the axes modeling the orientation of the pockets.The locations of the centers are random and are assumed to form a Poisson point process.This random distribution of spatial points satisfies the following: for any bounded, open (or closed) set D ⊆ R d , the number of points both in the random collection and in D is a random number N (B) that follows the Poisson distribution with mean λ|B|, where |B| is the volume of the set B and λ is called the intensity of the Poisson point process.Also, given N (B), the points {x i } ∩ B are independent and uniformly distributed in B. We may also assume that the radii {R i } and contrasts {τ i } are random and are i.i.d random variables with appropriate distributions.Since the inhomogeneities have a maximal and a minimal radius, we can assume that R m ≤ R i ≤ R M .In the same way, we have 0 < ℜ(τ m ) ≤ ℜ(τ i ) ≤ ℜ(τ M ) for the air bubbles and the brine pockets, with a similar relation for the imaginary parts.A simple calculation shows that where τ is the average contrast and v the average volume of the heterogeneities.In the same way, the variance of W is It is not expected that the constitutive parameters of one type of heterogeneities vary too much from one another, i.e. |τ | 2 ∼ |τ | 2 , and as a consequence (Var(W )) Since, as mentioned in the introduction, the volume fraction of the inclusions is at least 10%, then (λv) 1/2 ≥ 0.3, and this shows that W and (Var(W )) 1/2 are of the same order.This fact can help us investigate, in experimental settings, the validity of the assumption that the homogenized coefficients are just the averages of the random coefficients.We will see in the next section that the assumption holds provided that the wavelength is sufficiently large, and that the random medium fluctuations are not too strong.In the data provided in [15, figure 4], the real part of a diagonal component of the effective permittivity tensor ranges, depending on the frequency, from about 3.5 to 4.5 for a temperature of -14 degrees Celsius.More precisely, it is about 3.75 in the range 800Mhz-1Ghz.This means that in such configurations, the effects of the inclusions are not negligible, but are not too strong either, and correct the background of pure ice (with permittivity 3.14) by about 20%.Based on the fact that (Var(W )) 1/2 ∼ W , it seems therefore reasonable to assume that the fluctuations around the expected value are not too large as well, and that the homogenized coefficients are well approximated by the expected value of the random permittivity.At higher temperatures, the experimental permittivity can increase up to 5.5, in which case the effective coefficients should not be obtained by a simple averaging, but rather by a much more involved analysis [15].
We recover the scaling V (x/β) as follows: set ℓ c = R M , and suppose the intensity λ satisfies λ = λ 0 β −d .Then, setting x = Hx in W (x) yields, with the change of variables Above, we used the fact that a Poisson point process {x i } with intensity λ is statistically equivalent to a Poisson point process {βx i } with intensity λβ −d .The standard properties of Poisson processes show that V is stationary and has finite range correlations, so that (A1) and (A4) are satisfied.Assumption (A3) is also trivially verified since the real and imaginary parts of τ i are positive.Regarding (A2), even though for each realization the number of points x i in a bounded domain is finite, Poisson point processes allow clustering.This implies that the constant C 0 in assumption (A2) is not uniform in the randomness and that (A2) does not hold as stated.This issue can be fixed with additional technicalities that we just sketch here: replace the constant Its probability can be estimated by considering for instance the probability to find N > f (β) points in a domain of size β d , and can be shown to be very small as β → 0. Define then Ṽ = V on the complementary of Ω β in Ω, and Ṽ = 0 on Ω β .Choosing an appropriate f (β) and using Ṽ to model the heterogeneities, our proofs can then be directly adapted to recover the corrector result of Theorem 3.1.We present our main results in the next section.
3 Main Results

The homogenization and corrector results
We suppose that the source f is a smooth function with bounded support.In the regime we consider here, the solution u β to (2.3) converges as β → 0 (in proper sense) to the solution u of the following homogenized equation where n is defined in (2.4).We write u = Gf for the solution to (3.1), and have u ∈ L 2 (R d ) thanks to the artificial absorption α.We are mostly interested in the corrector where the randomness appears in the source term.Above, k 2 σ = σk 2 , and χ is the characteristic function of the sea ice layer S. To study the higher order correctors in u β − u, we rewrite (3.2) in a form that serves as the starting point of an iteration scheme: Performing another step of iteration, we get Defining v β = −Gk 2 σ χq β u, our first result below shows that, in terms of β, u β − u − v β is of order β 2 for d = 3 , and of order β 2 | log β| when d = 2.We will see in our second result (Theorem 3.2) that v β is of order β d/2 , and as a consequence v β is the leading corrector as β → 0. Theorem 3.1.Let d = 2, 3. Assume that the random field q satisfies (A1)-(A4) and let K be a bounded set in the upper half-space where and C is independent of β, α, k, κ m and κ e .
In the theorem, the domain K can be seen as the location of the detector where measurements are performed.In (3.4), we kept track of the dependency of the constants with respect to σ, k, and the absorptions κ e and κ m , in order to quantify how β must relate to these parameters for the corrector result to hold.Relatively to the size of u (measured here by u L 2 (S) ), we need C β,k,σ,κe β 2 ≪ 1, and remark first that for the wave to propagate in the sea ice layer, κ e and κ m have to be sufficiently small and of order k −1 .Indeed, when the absorption κ e is small compared to n 2 e (this holds in most practical configurations, see [15]), the imaginary part of n e behaves like κ e , and we expect the wave energy to decrease exponentially with a factor proportional to kκ e |x|.This means that kκ e has to be of order one for the wave not to go extinct at a shallow depth in the sea ice.We write then κ e ∼ κ m ∼ k −1 , and suppose that the variance of V is not necessarily small, say σ ≥ 1.The condition C β,k,σ,κe β 2 ≪ 1 then becomes, to leading order, The second condition boils down to β 2 σ 3 k 5 ≪ 1 since when β satisfies the latter condition, then the one involving log(βk) is in turn verified.Recalling that β = εη, and that k = 2π/η, we find the conditions εσ ≪ η 5/4 when d = 3, and εσ 3/2 ≪ η 3/2 when d = 2.These conditions hold when the wavelength is sufficiently large and the fluctuations not too strong.The proof of Theorem 3.1 is fairly standard and follows the lines of those e.g. of [2,3,10].It is based on estimates of fourth-order moments of q and on estimates on the Green's function of (3.1).The main differences with these references are that we work in an unbounded domain, which brings in additional technical difficulties, and that we keep track of important constants, in particular of k ≫ 1, in order to obtain a corrector result uniform in the main parameters in the problem.
The second result is a precise description of the limiting distribution of the normalized homogenization error (u β − u)/ β d , in the limit as β → 0 keeping the other parameters fixed.For its statement, we will need the square root of M , given by Denote by {α ij } the entries of M 1/2 , and let where W (1) y and W (2) y are standard independent multiparameter (y-parameter) Wiener processes [19].We have then the following result.Theorem 3.2.Under the same assumptions as in Theorem 3.1, we have, for any bounded set K in the upper half space, where G is the Green's function associated to problem (3.1).
The result above is a convergence in distribution of functions in L 2 (Hilbert) spaces, see section 4.5 and [23] for an introduction on the subject.In Theorem 3.2, the Wiener integrals in the limiting corrector are simply, after integration in x against a test function, normal random variable with zero mean and appropriate variances.The first-order corrector v is therefore a mean-zero Gaussian field, whose information is contained in correlations of the form In the next section, we approximate these correlations in the high frequency limit η → 0, assuming the condition C β,k,σ,κe β 2 ≪ 1 holds for the corrector result to be true.

The high frequency limit
With Theorems 3.1 and 3.2 at hand, we approximate the wavefield u β by u β 0 = u + β d/2 v, where v is formally the solution to, recalling that k = 2π/η, where we have introduced The main tool to study the high frequency limit η → 0 of (3.7) is the Wigner transform, defined by, for a function w: See [20,12] for an introduction on Wigner transforms.Rigorous mathematical analyses with Wigner transforms are typically difficult and technically involved, all the more in domains with sharp interfaces, and are beyond the scope of this work.We hence decided to remain at a formal level in the present section.We refer to [11,21] for rigorous derivations.The correlation and an inverse Fourier transform yields Our next result provides us with an asymptotic description of E{W[v]} as η → 0. Before stating it, we need to introduce a few notations.For |ξ| ≤ k e , let k j (ξ) = (k 2 j − |ξ| 2 ) 1/2 , for j = {0, e, 2}.We recall that k 0 < k e < k 2 , and we set k 0 (ξ) = 0 when |ξ| > k 0 .The reflection-transmission amplitudes at the interfaces z = 0 (associated to j = 0) and z = −L (associated to j = 2) are defined by The boundary values of a function W(x, p) at these interfaces are denoted as follows: write where the upperscript ± refers to the upper/lower limit.We have then the following theorem.
, where o(η) represents a term that converges to zero as η → 0 in the distribution sense, and where W satisfies equipped with the following reflection-transmission conditions at the interface z = 0, and at z = −L, The (formal) proof is direct and follows the lines of [5].The interpretation of Theorem 3.3 is the following: as the homogenized solution u propagates in the sea ice layer, it interacts with heterogeneities at scales much smaller than the wavelength.This results in an isotropic radiation (modeled by the delta function δ(|p| 2 − k 2 e )) with intensity |u(x)| 2 , which then propagates according to ray theory, and is refracted at the interface between the air and sea ice and the interface between sea ice and seawater following Snell-Descartes laws.Note that, even though the refractive index is complex in the sea ice and seawater, we still obtain the usual Snell-Descartes laws.This is because the absorption is small compared to the real part of the index.When this is not the case, Snell-Descartes laws have to me modified, see [7].Moreover, we do not need to consider the critical case |k ⊥ | = k 0 in the boundary conditions since this situation corresponds to a set of measure zero in R d .There is no critical case at the bottom interface since the real part of the refractive index in seawater is larger than the effective index in sea ice, and therefore k 2 > k e .
Using the stationarity of the random medium, the correlations can be calculated as follows.Recalling the form of the source term (2.2), and indexing u β as u β t , we have, invoking ergodicity, From a practical viewpoint, what is measured is for some T 0 sufficiently large and a few frequencies ω in the bandwidth.The correlation C(x, y) is asymptotically modeled by The next step is the resolution of an inverse problem with data C(x, y) at the detector and model C 0 (x, y), with the goal of recovering some information about the sea ice, namely the thickness of the layer and the homogenized coefficients.This will be addressed in future works.The next section is devoted to generalizations of Theorems 3.1 and 3.2.

Generalizations
A multi-layer configuration can be considered as follows: suppose the sea ice consists of N L layers defined by and assume there are N i H different types of heterogeneities in the layer i.We then define the new random field by where the V ij are independent and have a typical scale much smaller than the wavelength.We set ℓ c to be the largest of these scales.Above, χ i is the characteristic function of L i .When the functions ℓ i vary slowly compared to the scale rescaled variables, then the situation is essentially similar to the flat interface case since the scales of variation of the surfaces and of the heterogeneities are well separated.It is then expected that, asymptotically in β, In the latter, u and G are the solution and the Green's function of the multi-layers multi-species homogenized Helmholtz equation, and dW ij y has the form (3.5), where the α coefficients are the entries of square root of the matrix M ij , which is defined in a similar way as in (2.5).
Non-flat interfaces can also be considered; suppose that each interface profile consists of a slowly varying part, i.e. varying on scales much larger compared to the wavelength, augmented with a highly oscillatory (possibly random) part.When the oscillations are much faster than the wavelength, homogenization theory predicts a limit model where the rough boundary is replaced by the slow varying profile that is surrounded by a boundary layer where an equation with effective coefficients is solved.Outside of this layer, the original equation is satisfied with appropriate continuity conditions.See e.g.[22] for more details and the references therein.In terms of (3.9), this means that u and G are replaced by some homogenized versions u e and G e .When the oscillations of the rough boundary are of the same order as the wavelength, the situation is more complicated.Nevertheless, if the amplitude of the fast fluctuations is sufficiently small, a transport model for u(x)u * (y) can be obtained provided the high frequency hypothesis holds in the layer of interest.See [4] for more details.One of our next objectives is to compare numerically these various boundary effects with the volume scattering modeled the corrector v.
The rest of the paper is devoted to the proofs of our main results.

Proofs of the homogenization and corrector results
We start with some preliminaries about Green's functions.

The Green's function of the Helmholtz equation
We derive some estimates for the Green's function G(x, y) of the homogenized equation, i.e. the solution to ∆ x G(x, y) where n is defined in (2.4), by exploiting the layered structure of the underlying medium and the physical absorption in the middle layer (sea ice).Note that the reciprocal property G(x, y) = G(y, x) holds.Let Φ = Φ(x, y) be the Green's function in free space, that is where n e is the homogenized complex refractive index of sea ice, see (2.4).The function Φ admits the following explicit formulas: where H (1) 0 is the zero order Hankel function of the first kind.We have the proposition below.
Proposition 4.1.Consider the Green's function G(x, x 0 ) solution to (4.1) and let K be a bounded set in the upper half space With (2.4), suppose that n e > n 0 and that n 2 > n 0 .Set p(x, x 0 ) = G(x, x 0 ) − Φ(x − x 0 ) for x 0 ∈ S. Then there exists a generic constant C independent of k and κ e (and α) such that we have the following L ∞ estimate sup and G verifies sup where H s is the usual Sobolev space.
The proof is postponed to section 4.6.When d = 2, we will use the result below, proved in section 4.7.1:Lemma 4.2.For d = 2, we have the estimate, where C > 0 does not depend on k and n e .
In the rest of the paper, we will use the notation R(x − y) := E{q(x)q * (y)}.
Note that, due to assumption (A4) and a similar relation to (2.6), R is integrable over R d .

Error estimates
The goal of this subsection is to prove Theorem 3.1.The strategy of proof is an adaptation of that of [3].The starting point is the relation.
With the goal of controlling the last term on the right above, we get a first direct estimate for u β − u with the lemma below.We recall that χ is the characteristic function of the sea ice domain S.
Lemma 4.3.Let g ∈ L 2 (S) and v satisfy Then, we have the estimate Proof.Multiplying (4.8) by v * , integrating over R d and taking the imaginary part leads to and it suffices to use the Cauchy-Schwarz inequality to conclude.
With v β = −Gk 2 σ χq β u, we then note the following: In view of Lemma (4.3) and assumption (A2), we find The next step is to estimate v β .We use for this the following lemma, whose proof is postponed to the end of the section.
With the shorthand u 2 , we then find which yields the homogenization result provided (1 satisfies the same bound; however, the same argument does not yield uniform in α control of u β − u − v β L 2 (K) , for a set K in the upper half space, because the absorption in the upper space is precisely α.We will use the above estimate together with the next lemma in order to control the higher order terms in (4.7).Fourth order moments of q β are needed in the proof, which is given in section 4.4.Below, L(E, F ) denotes the space of bounded operators from the Banach space E to the Banach space F .Lemma 4.5.Assume (A1)-(A4) hold.Then, for any bounded set K in the upper half-space, and k,σ,κe = Cσ 4 k d+4 /κ 2 e , and C k,κe = k 2d−4 /κ 2 e .
We can now estimate the last term in (4.7).We have, with (4.9), (A2), and the Cauchy-Schwarz inequality, In terms of β, this term is of order β d .Regarding the second term in the r.h.s of (4.7), we find, with (4.10), This completes the proof of Theorem 3.1 by inspection.

Proof of Lemma 4.4
From direct calculations, we have Integrating over x first and using Cauchy-Schwarz inequality, we find We can recast the remaining integral as where R β is a short-hand notation for R(•/β).Note that R integrable, and h is square integrable.We have then, using Young's inequality, and obtain finally Using the first estimate of G in (4.5), we get the desired result and the proof is complete.

Proof of Lemma 4.5
We start with the proof of (4.9).We have, after a Cauchy-Schwarz inequality, for all x ∈ K and all f ∈ L 2 (S), Taking expectations on both sides yields In the last line, we used Young's inequality.The second estimate of (4.5) and the symmetry of G then completes the proof of (4.9).Regarding (4.10), we need fourth order moments estimate for mixing random fields.We have, for q satisfying (A1)-(A4), where ρ = ̺ 1 2 (•/3) and α j ∈ {r, i} for the real and imaginary parts of q; see [16] for a proof.Note that ρ is bounded by one and integrable in R d according to (A4).Decomposing q into real and imaginary parts, there are 16 terms to estimate, all of which are treated in the same way with (4.11).We therefore pick one as an example.Direct computation shows Taking expectation on both sides yield Following (4.11), we need to control three integrals.The most singular one corresponds to the combination of arguments localizing the Green's functions on the diagonal, and is given by We then decompose G as in Proposition (4.1) into G(x, y) = Φ(x − y) + p(x, y), where Φ is the free space Green's function solution to (4.2) and p a regular function.We split the integral above further into two terms.The first one is obtained by replacing the first G by Φ above.Then, after a change of variables, this term becomes When d = 3, the last term is controlled by Here, we used the fact that ρ is integrable and bounded in R The contribution of p(y, z) is controlled as follows.We have, for all x ∈ K: We used Young's inequality in the last step.With (4.4) and (4.5), the term above is controlled by Cβ d k 3d/2−4 κ −2 e .Gathering all previous results, we find that the contribution associated with the first term in (4.11) is of order σ ) for d = 2.The other terms in (4.11) lead to weaker contributions in β.We take the following term as an example: after integrating in x, We recast the integral as Multiplying by k 8 σ sup y∈S G(•, y) 2 L 2 (K) , and accounting for (4.5), this term generates of contribution of order k 2d σ 4 β 2d κ −4 e , which is weaker than the other leading terms since k ≫ 1.By exactly the same method, we see that the third term on the right hand side of (4.11) leads to a contribution of order β 2d .Combining all the results above, we proved (4.10).

Asymptotics of the corrector
We address here the convergence of the corrector v β = −Gk 2 σ χq β u as β → 0, keeping the other parameters fixed.We introduce first some mathematical preliminaries.
Preliminaries.Let H be a Hilbert space, and X : (Ω, F , P) → H be a random variable with value in H. Then X induces a (Borel) probability measure on H as follows: for each Borel set B ∈ H, define λ X (B) := P(X ∈ B).
A family {X β : Ω → H} β∈(0,1) of random variables in H is said to converge in distribution to X : Ω → H, as β → 0, if the family of probability measures induced by {X β } on H, denoted by λ X β , converges to that induced by X, denoted by λ X .In other words, for any Borel measurable subset B ⊆ H, the real numbers λ X β (B) converges to λ X (B).An important and very useful criterion for the convergence in distribution of X β to X is given as follows.
We refer to [23] for its proof.
Proposition 4.6.Suppose H is a separable Hilbert space.Let •, • denote the inner product on H. Then X β converges in distribution to X, provided that (1) for each h ∈ H, the random variables X β , h converges in distribution to X, h .
The norm on H s (K) is given by the formula f 2 In view of the fact that the embedding H s (K) ֒→ L 2 (K) is compact, a very simple yet useful criterion for tightness of measures on L 2 is given as follows; see [18,Theorem A.2]. Proposition 4.7.Let {X β } β∈(0,1) be a family of random variables on the separable Hilbert space L 2 (K).Suppose that there exists C > 0 and s ∈ (0, 1), independent of β, such that Then the family {λ X β } β∈(0,1) of probability measures induced by {X β } is tight.
Next, we prove Theorem 3.2.We go back to the expansion formula (3.3) and realize that, according to Theorem 3.1, It follows that the limiting distribution in L 2 (K) of the normalized homogenization error In view of the criteria in Proposition 4.6 and Proposition 4.7, it suffices to prove the results in Lemma 4.8 below.Recall the notations in Theorem 3.2: M 1/2 = {α ij } is the square root of the matrix M defined in (2.5), and W y is defined by where W (2) y are standard independent multiparameter (y-parameter) Wiener processes, and W y := (W where G is the Green's function of the homogenized equation 3.1.We have then: Lemma 4.8.Under the same assumptions of Theorem 3.2, we have: (1) For each ϕ ∈ L 2 (K), (2) Let d = 2, 3 and s < (4 − d)/2.Then there exists constant C > 0 such that Proof.For item one, we compute Write then for simplicity m 0 (y) = u(y)m(y)χ(y).We address the convergence of ζ β := β −d/2 (v β , ϕ) by considering the real vector (ℜ(ζ β ), ℑ(ζ β )) T .For this, we consider random vector h β of the form, with Q β (y) = (q r ( y β ), q i ( y β )) T , where A = {a ij } is matrix whose entries are real L 2 (S) functions supported in S.These are integrals involving the product of an L 2 (S) function and the random oscillatory function q β α , α = {r, i}, which is a rescaled version of the short range correlated random field q α (x, ω).Adapting the results of [2,3] to our case, we find that where M is defined in (2.5).With the multiparameter Wiener process W y = (W y , W y ) T defined earlier, we find that the distribution on the right hand side is realized by The left hand side can be written as a sum of several Wiener integrals in terms of the independent components W y ; those integrals are simply random normal variables with zero mean and they have appropriate joint variances.
To verify (4.12), it suffices to write ζ β as We then obtain (4.12) by inspection.
Next, we show that the measures generated by v β /β d/2 are tight in L 2 (K), by proving (4.13).Here, K is any bounded set in the upper half space.By definition, where the second term is the H s seminorm given by Taking expectations, we have Then, with the Cauchy-Schwarz inequality, sup With (4.6) and plugging the latter estimate in (4.14), we have By similar arguments as in the proof of Lemma 4.4, we get to conclude that Together with Lemma 4.4 which controls the L 2 norm of v β , this completes the proof.
An immediate result of the lemma is that Since β − d 2 (u β − u) has the same limiting distribution, the proof of Theorem 3.2 is complete.

Proof of Proposition 4.1
Write x = (x ′ , z) with x ′ = (x 1 , x 2 ) when d = 3, and x ′ = x 1 when d = 2. Taking the Fourier transform of (4.1) in the x ′ variable (with dual variable ξ ∈ R d−1 ), we find where x 0 = (x ′ 0 , z 0 ) and the complex number k 2 (z, ξ) is defined by We write the dependencies of Ĝ as Ĝ(z, z 0 , ξ, x ′ 0 ).Although it is possible to obtain exact expressions for Ĝ, we follow a different route that seems somewhat more direct and would apply if the refractive index in the sea ice layer depends also on z.We start with the case z 0 ∈ (−L, 0).
L 2 estimates.Inequality (4.19) allows us to estimate p for |ξ| ≤ kn 0 in L 2 as follows: with the change of variables ξ → kξ in the r.h.s, we obtain with We only treat the term in n 0 in the r.h.s of (4.22) since the other one is handled similarly.Some direct algebra involving Taylor expansions shows that there exists C > 0 that is independent of ξ and α such that With κ α = κ e + α, we are then left with estimating the integral where C is independent of κ e and n 0 .Above, we used the fact that n 0 < n e .This gives us a bound for p in L 2 in the region where |ξ| ≤ kn 0 .Assume now without lack of generality that n e > n 2 (the opposite case n e < n 2 can be handled similarly).For |ξ| ≥ kn e , we combine (4.20) and (4.21) to arrive at The first term in the r.h.s is easily taken care of by (4.23) and the same argument as before.
For the second one, with the notation κ 2,α = κ 2 + α, we need to estimate the integral: .
Using the change of variable |ξ| → κ α r + n e , we see that the above is bounded by for any δ > 0. It remains to treat the case kn 0 ≤ |ξ| ≤ kn e .We start by using (4.21) to obtain, after multiplying by k Using (4.20) and Young's inequality, we get which, together with (4.21), leads to Integrating over kn 0 ≤ |ξ| ≤ kn e , and using that (k 2 n 2 e − |ξ| 2 )/|k 1 | 2 ≤ 1 together with (4.23), we find Collecting results, we have proven (4.3).We turn now to the L ∞ estimates.
L ∞ estimates.We control p L ∞ by p L 1 .For z ∈ [−L, 0], the solution to (4.16) reads and as a consequence |p(z)| ≤ |p(0)| + |p(−L)|.As before, assuming without lack of generality that n e > n 2 , we split the integration into three parts, |ξ| ≤ kn 0 , kn 0 ≤ |ξ| ≤ kn e , and |ξ| > kn e .We estimate the part |ξ| ≤ kn 0 using (4.19).We only treat p(0) and the term involving k 0 in the r.h.s since all other calculations are similar.Using (4.23) and the fact that n 0 < n e , we then find For the part |ξ| ≥ kn e , we only treat as above the k 0 part of the r.h.s in p(0) and write for some R > 0 that will be defined below.The first term is easily estimated using (4.23): For the second term, we need to refine estimate (4.23) for large |ξ|, and find, for appropriate C R and R sufficiently large, With the change of variables |ξ| → κ α r + n e in the second term of (4.25), we find, using (4.26), The last piece kn 0 ≤ |ξ| ≤ kn e is estimated by using (4.24).Since (n 2 e − |ξ| 2 )/|n e | ≤ (n 2 e − |ξ| 2 ) 1/2 , the k 0 part in the r.h.s in p(0) is bounded by This proves the first estimate in (4.4).The second one follows from the observation that, for z > 0 and z 0 ∈ (−L, 0), |p(z)| = |p(0)e ik 0 (ξ)z | ≤ |p(0)|.
We turn now to the case z 0 > 0 and prove the second estimate in (4.5).The calculations are very similar and we only provide the main ideas.
Estimates for z 0 in the upper half-space.We start with (4.15) with now z 0 > 0, and derive first boundary conditions at x = 0 and x = −L.The one at x = −L is the same as when z 0 ∈ (−L, 0).For the one at x = 0, we notice that ∂ z Ĝ has a jump of value e ix ′ 0 •ξ at z = z 0 .Accounting for this, solving (4.15) for z > z 0 and z ∈ (0, z 0 ), and eliminating the unknown reflection coefficient at x = 0 yields the boundary condition In the region z ∈ (−L, 0), Ĝ satisfies the Helmholtz equation with wavenumber k 1 and has the boundary conditions specified above.We can then estimate Ĝ in L 2 by reproducing the steps for the analysis of p in the previous section.In particular, we replace, in (4.17 The last piece kn 0 ≤ |ξ| ≤ kn e is estimated by Gathering previous estimates, we obtain the second part of (4.5).
We continue the proof of lemma by estimating the free Green's function Φ ξ .
Estimates on Φ ξ .In order to obtain the first estimate in (4.5), we compute With the change of variables |ξ| → κ α r + n e , we find .
It is direct to see that the above integral is finite for r > 0 and then bounded by Cκ −1/2 e .We focus therefore on the piece r ∈ (−n e /κ e , 0).Since the function x → .If not, we estimate the integral by, using (4.27) and the decrease of Collecting estimates, and using κ e ≤ 1, we finally arrive at sup Together with (4.3), this yields the first estimate in (4.5).
Estimate in H s (K).We write G(x, y) = Φ(x − y) + p(x, y) and start with Φ, for which a direct computation shows that it does not belong to the Sobolev space H 1 (R d ).Nevertheless, we can check that for any s < (4 − d)/2, we have Φ ∈ H s (R d ).Indeed, since we have from (4.2) that The zero order Hankel function H  The reflection-transmission boundary conditions are obtained for instance by adapting the results of [24] to our layered case.This ends the proof.

Conclusion
We have derived in this work, under appropriate conditions on the wavelength, the typical length scale of the fluctuations and their variances, an asymptotic model for the first-order corrector to the homogenization limit.We have considered a simplified model where the inhomogeneities occupy a single layer with flat interfaces, and where the propagation of waves is described by a random Helmholtz equation.The corrector takes the simple form of a Gaussian field whose statistics depends on the homogenized solution and the homogenized Green's function.In addition, when the thickness of layer is sufficiently large compared to the wavelength, we have obtained a transport model for the correlations of the correctors.Finally, we have explained how to generalize our results to more complicated settings of multi-layers with rough boundaries.
This work is first step towards the understanding of the true physical problem and it leaves many questions open.We describe next some of those problems.Firstly, what is the form of the corrector when the fluctuations of the media have stronger strength?In this case, the homogenized model involves the resolution of a cell problem and the characterization of the corrector becomes much more difficult.Secondly, how significant are boundary effects?It is unclear whether the main contribution to the backscattered signal comes from volumic scattering by the heterogeneities in the sea ice or from surface scattering from rough interfaces, in particular the sea ice / sea water interface that is the most relevant for our purpose.Finally, how all of this translates to the true physical situation that involves electromagnetic waves and Maxwell's equations?These questions will be addressed in future works.