Self-averaging of kinetic models for waves in random media

Kinetic equations are often appropriate to model the energy density of high frequency waves propagating in highly heterogeneous media. The limitations of the kinetic model are quantified by the statistical instability of the wave energy density, i.e., by its sensitivity to changes in the realization of the underlying heterogeneous medium modeled as a random medium. In the simplified It\^o-Schr\"odinger regime of wave propagation, we obtain optimal estimates for the statistical instability of the wave energy density for different configurations of the source terms and the domains over which the energy density is measured. We show that the energy density is asymptotically statistically stable (self-averaging) in many configurations. In the case of highly localized source terms, we obtain an explicit asymptotic expression for the scintillation function in the high frequency limit.


Introduction
Let us consider the following scalar wave equation for the pressure potential p(τ, x, t): 1 c 2 (x, t) where τ is time, (x, t) ∈ R d ×R denote the spatial variables, ∆ x is the Laplace operator in the transverse variables x, and c(x, t) is the local sound speed.
Our objective is to understand the properties of p(τ, x, t) when c(x, t) is a highly oscillatory random field and the initial conditions for p(τ, x, t) oscillate at the same frequency. The analysis of high frequency waves in random media based on (1) is extremely complicated and still not totally established mathematically. Since the wave field is oscillatory, its (weak) limit typically misses most of the energy of the wave field p. Kinetic models are then used to capture the energy density of the wave fields; see e.g. [6,7,12,17] for rigorous results, [3,19] for more formal derivations, and [11,14,16,20] for references in the physical literature.
The validity of the kinetic model is limited by its statistical instability, namely by its variability when the realization of the underlying random medium is changed. In many situations, the energy density is self-averaging [2,6,7], which means that the energy density measured (averaged) on a (sufficiently large) domain is asymptotically, as the frequency goes to infinity, independent of the realization of the random medium. The above results often require that the domain of measurements be of size independent of the wavelength and that the source term for the kinetic model be sufficiently smooth.
In this paper, we are interested in the statistical stability of such kinetic models in a very simplified regime of wave propagation, namely the Itô-Schrödinger regime. The latter regime arises when the wave field is a very narrow beam propagating in the direction t and the sound speed c(x, t) oscillates more rapidly in the direction t than it does in other directions. Such assumptions are valid in somewhat restrictive practical settings. However, this regime of wave propagation is relatively simple to analyze mathematically and provides interesting qualitative answers regarding the statistical stability of more general kinetic models.
The validity of kinetic models has been analyzed numerically in several settings [8,9,10], with quite good agreements with the energy density given by wave equations of the form (1). Such kinetic models may then be used to solve inverse problems, where constitutive parameters in the transport equation modeling e.g. buried inclusions or statistics of the random medium, are reconstructed from available boundary measurements. We refer the reader to [9,10] for reconstructions based on synthetic (numerical) data and to [5] for kinetic reconstructions from experimental data in the micro-wave regime; see also [4] for a review on the use of kinetic models in the imaging of buried inclusions. These studies show that the kinetic models perform relatively well. Their limitations are almost entirely caused by our lack of knowledge of the random medium, which generates some statistical instabilities in the measurements. Understanding these instabilities will allow us to improve on the reconstructions and to have a better understanding of the maximal resolution that can be achieved.
Itô-Schrödinger regime. In the Itô-Schrödinger regime, we introduce ψ(x, t; κ) as where c 0 is the background sound speed, assumed to be constant. Thus ψ represents waves at position (x, t) propagating with frequency ω = c 0 |κ|. After appropriate scalings and simplifications, the wave field ψ satisfies the following Itô-Schrödinger stochastic partial differential equation: Since κ plays no significant role in the sequel, we set it to κ = 1. Here, B(x, dt) is the standard Wiener measure, whose statistics are described by where E is mathematical expectation with respect to the measure of an abstract probability space on which B(x, dt) is defined and t ∧ t ′ = min(t, t ′ ). We shall not justify (3) from (1). See [1] for a justification in one dimension of space and [2] for the scaling arguments leading to (3). For our purposes, ψ η (x, t) satisfies a wave equation with highly oscillatory coefficients oscillating at a frequency inversely proportional to the small parameter η ≪ 1. We assume that ψ η (x, 0) also oscillates at a frequency comparable to η −1 and are interested in the properties of the wave field as η → 0. Because the field oscillates rapidly, its weak limit is of little interest. A more interesting quantity is the energy density of the waves |ψ η | 2 (x, t), or the probability density in the context of quantum waves. Because the energy density does not satisfy a closed-form equation, it is more convenient to analyze energy densities by introducing the following Wigner transform of the wave field: where ψ η denotes complex conjugation of ψ. Note that R d W η (t, x, k)dk = |ψ η (x, t)| 2 by inverse Fourier transform so that W η may be seen as a phase space (microlocal) decomposition of the energy density. Let ψ η (x, 0) be a sequence of functions uniformly bounded in L 2 (R d ), η-oscillatory, and compact at infinity in the sense of [13], i.e., such that for every continuous compactly supported function ϕ on R d , we have: A practical sufficient condition is that ψ η (x, 0) is compactly supported and η∇ψ η (x, 0) is square integrable with L 2 (R d )-norm bounded independently of η. Then, we have the following convergence result [13,15]: The Wigner transform W η (0, x, k) converges, after possible extraction of subsequences, in the space of distributions D ′ (R 2d ) to a Radon measure W 0 (0, x, k), and moreover, we have In other words, the limiting Wigner transform captures all the energy of the incident wave field ψ η in the limit η → 0.
Kinetic Model. Upon using the Itô formula, we obtain that the average Wigner transform a η (t, x, k) = E{W η (t, x, k)}, solves the following kinetic equation where we assume that ψ η (x, 0), whence a η (0, x, k), is deterministic; see e.g. [2] for the details of the derivation. We have defined R 0 = R(0) andR(k) as the Fourier transform of R(x), with the convention that Since R(x) is a correlation function,R(k) is non-negative by Bochner's theorem. For the rest of the paper, we assume thatR(k) Note that R 2d a η (t, x, k)dxdk is independent of time so that the total energy of the initial condition is preserved by the transport evolution.
Scintillation. The validity of the kinetic model (8) to describe the ensemble averaging of the phase space energy density of the wave field is trivial in the Itô-Schrödinger regime: the kinetic model (8) is here exact for all η ≥ 0, unlike what happens in other regimes of wave propagation [6,7,19]. It remains however to understand how stable it is. In other words, how good an approximation is a η (t, x, k) of the random field W η (t, x, k). A natural object in the study of the statistical stability of W η is the following covariance function: We refer to this function as the scintillation function, in analogy to how stars are perceived to twinkle because the realization of the atmosphere changes in time.
We shall see that the size of the scintillation function crucially depends on the smoothness of the initial conditions ψ η (x, 0) and a η (0, x, k) and on the support of the domain over which the energy density is averaged. The effect of the averaging will be quantified by measuring J η in appropriate (weak) norms.
One of the main advantages of the Itô-Schrödinger regime of wave propagation is that J η (t, x, k, y, p) satisfies a closed form equation. Another application of the Itô formula [2] shows that J η is the solution of the following kinetic equation: with vanishing initial conditions J η (0, x, k, y, p) = 0, where In the absence of the operator K η , the variables (x, k) and (y, p) remain uncoupled in (11) and the scintillation vanishes. Scintillation is created as the waves propagate through the random medium with a rate of creation proportional to K η a η ⊗ a η . Notice that K η involves a highly oscillatory integral. Outside of the diagonal x = y, this oscillatory integral is small, whereas in the vicinity of the diagonal x = y, it is not. We thus observe that K η h is small when h is smooth and large when part of h is concentrated near x = y.
Outline. The rest of the paper is structured as follows. The main results of the paper are summarized in section 2. We obtain estimates for J η in various norms, and in the specific case of initial conditions for a η of the form a η (0, x, k) = δ(x)f (k), show that η −1 J η converges to a measure J solving an explicit kinetic equation. Section 3 presents stability estimates for the scintillation operator K η defined in (12) and for the kinetic equations (8) and (11). The proof of the stability estimates for J η are given in section 4 whereas the proof of convergence of η −1 J η when a η (0, x, k) = δ(x)f (k) is given in section 5.

Main results
Let ψ η (x, 0) be a sequence of η−oscillatory, compact at infinity, functions uniformly bounded in L 2 (R d ). This is the case of interest for us here, where we can define the Wigner transform (5) and pass to the high frequency limit η → 0 while still ensuring that energy is conserved as in (6). We are interested in quantifying the statistical stability of the Wigner transform W η (t, x, k) and do so by analyzing the scintillation function J η defined in (10).
We present two results. The first result proposes an upper bound for J η in different norms and for different initial conditions ψ η (x, 0). The second result analyzes the convergence properties of J η as η → 0 for initial conditions of the form a η (0, x, k) = δ(x)f (k), which correspond to localized sources at position x = 0 radiating energy smoothly in wavenumber k. In this context, we will show that J η is of order O(η) and will obtain the limit of η −1 J η as η → 0.
Some typical initial conditions. Let us consider initial conditions ψ η (x, 0) oscillating at frequencies of order η −1 and with a spatial support of size η α for 0 ≤ α ≤ 1. The parameter α quantifies the macroscopic concentration of the initial condition.
The simplest example is a modulated plane wave of the form: where χ(x) is a smooth compactly supported function on R d . The direction of propagation is given by k 0 . Note that the above sequence of initial conditions is indeed uniformly bounded in L 2 (R d ), compact at infinity, and η-oscillatory.
As another example of initial conditions, we consider where J 0 is the zero-th order Bessel function of the first kind. Such an initial condition is supported in the Fourier domain in the vicinity of wavenumbers k such that |k| = |k 0 | so that ψ (2) η emits radiation isotropically at wavenumber |k 0 |; see [8,9] for more details. We again verify that the above sequence of initial conditions is indeed uniformly bounded in L 2 (R d ), compact at infinity, and η-oscillatory. For this, we use that J 0 (z) = 2 πz cos(z − π 4 ) + O(z −3/2 ). Domain of measurements. For the above initial conditions for ψ η , we are interested in the corresponding Wigner transform W η (t, x, k) and scintillation function J η . It turns out that J η is itself oscillatory so that its size depends on the scale at which it is measured. In order to capture this scale, we introduce a test function ϕ ∈ S(R 2d ), a fixed wavenumber k 1 ∈ R d , and define We then denote by ·, · the duality product S ′ (R n )-S(R n ) for n = 2d or n = 4d and want to quantify W η , ϕ η,s 1 ,s 2 , the energy density averaged over a domain (in the phase space) of width η s 1 in space and η s 2 in wavenumbers. By using the Chebyshev inequality, we obtain the following estimate on the probability that W η deviate from its ensemble average a η : Here, a ⊗ a(x, k, y, p) = a(x, k)a(y, p). In other words, when the above right-hand side converges to 0, then we find that W η (t), ϕ η,s 1 ,s 2 converges in probability to 0, which implies that W η (t) converges weakly and in probability to 0. The measured energy density is thus asymptotically statistically stable. A very relevant practical question pertains to the largest values of s 1 and s 2 that can be chosen so that the Wigner transform is still statistically stable in the limit η → 0. We are now ready to state our main theorem on this issue.
Bounds for the scintillation function. For any ϕ(x, k) ∈ L 2 (R 2d ), let F x ϕ(u, k) and F k ϕ(x, ξ) be the Fourier transforms of ϕ in the first variable only and in the second variable only, respectively. We also denote by a b the inequality a ≤ Cb, where C > 0 is some universal constant. Then we have the following result: Theorem 2.1 Let ψ η (x, 0) be a sequence of functions uniformly bounded in L 2 (R d ), compact at infinity, and η-oscillatory. Let a η (0, x, k) be the corresponding sequence of Wigner transforms given by (5). We assume that F x a η (0) and F k a η (0) are integrable functions and that for some α ∈ R and β ∈ R. Then we find that Of interest here is the following corollary: Corollary 2.2 Let ψ η (0) be given by one of the expressions in (13) or (14).
We can deduce the following results from the above corollary. In what follows, we consider that averaging takes place over a large domain of wavenumbers so that s 2 = 0, as e.g., in spatial measurements of the physical energy density.
Support of the sources. Let us assume that the spatial support of the domain of measurements is large so that s 1 = 0 as well. Then we find that In other words, the scintillation is of order O(η d ) when α = 0, which corresponds to a large support of the initial source term. This corresponds to the ideal case where the scintillation is smallest. In such a setting, we obtain that W η − a η , ϕ is of order η d 2 . This is the most stable situation. For a very narrow support of the initial source term comparable to the correlation length of the medium, namely when α = 1, we obtain that the scintillation is of order O(η) so that W η − a η , ϕ is now of order η 1 2 . We thus obtain statistical stability of the energy density generated by a very localized source term whose radiation pattern in k is smooth, although the statistical instability is much larger than in the case α = 0. We know that for sources that are highly localized both in space and in wavenumbers, the scintillation does not converge to 0 and the energy density is not asymptotically statistically stable; see [2]. Such highly localized initial conditions would correspond to a choice α = β = 1 in Theorem 2.1. We will confirm in the next theorem that the order O(η) above is optimal.
Small domain of measurements. Conversely, we can consider the case of a source term with a large support, which corresponds to α = 0, and a very small measurement domain. In this setting, we find that This means that the energy density becomes asymptotically statistically stable as soon as it is measured over an area that is large compared to the correlation length of the medium. This is an optimal result of self-averaging as we cannot expect the energy density to be statistically stable point-wise, or when averaged over sub-wavelength domains. The above result, which is based on estimating K η in (12) in appropriate norms, improves on estimates obtained in [2,18].
We can also consider intermediate situations where both the source and the measurement domain have small support. In that case, the optimal estimate for the scintillation depends on whether α < s 1 or s 1 < α. These results are in fact optimal when the source term and the domain of measurements are located at the same place. Such a geometry explains why we do not obtain scintillation proportional to η d(1−α−s 1 ) . We should obtain better estimates when the domain of measurements and the source term are not centered around the same point, though this cannot be inferred from our current results.
Convergence of scintillation. Let us consider the case of initial conditions of the form (13) or (14) with α = 1, i.e., for tightly localized source terms, in (transverse) dimension d ≥ 2. The Wigner transform of such source terms converges in the limit η → 0 to a distribution of the form δ(x)f (k), where f (k) is a smooth function when χ(x) is smooth [15]. We consider the kinetic equations with such initial conditions and obtain the following result. Theorem 2.3 Let J η be the solution of (11) with the initial condition in (8) given by a η (0, x, k) = δ(x)f (k) for some smooth function f (k) in dimension d ≥ 2. Then η −1 J η (t) converges in the space of distributions uniformly in time to the limit J(t), which solves the following kinetic equation where The above theorem should be interpreted as follows. As time propagates, the transport ballistic part a 0 (t, x, k) = e −R 0 t δ(x − tk)f (k) creates some instabilities, which converge after appropriate scaling to J 0 (t). The scintillation thus generated is then transported by the transport equation (21). We also observe that the error estimate of order O(η) in (19) with α = 1 is optimal.

Functional setting and stability estimates
In preparation for the proof of the theorems and the corollary presented in the preceding section, we prove here some stability results for the transport equations (8) and (11) and for the scintillation operator K η .
We denote by F the operator of Fourier transform with respect to all variables of the function on which it applies. For 1 ≤ p ≤ ∞, we introduce X p as the subspace of tempered distributions in S ′ (R 4d ) such that for 1 ≤ p < ∞ and We also define Y p as the subspace of tempered distributions in S ′ (R 2d ) such that for 1 ≤ p < ∞ and Finally, we define Y as the subspace of tempered distributions in S ′ (R 2d ) such that g Y = sup Morally (though this is inexact), the space X 1 corresponds to scintillation functions that are integrable in one spatial variable (bounded in the corresponding dual variable v) and bounded in another spatial variable (integrable in the corresponding dual variable u). It is this boundedness that allows us to obtain the result (20) in the presence of small domains of measurements. In contrast, X ∞ corresponds to scintillation functions that are integrable in both spatial variables (bounded in u and v), which allows us to get the result (19).
The above spaces are well-adapted to the estimation of the scintillation operator K η . More precisely, we have the following result: (ii) Let µ ∈ Y p and ν ∈ Y . Then Proof. With obvious notation, we recast K η = ǫ i ,ǫ j ǫ i ǫ j K ij η . Let h ∈ X p . Then we have so that using the Hölder inequality with 1 = 1 This proves (i). Let now h := µ⊗ν. Upon performing the change of variables w → ηw, we have which concludes our proof. We need stability estimates for the kinetic equations. We start with the first kinetic equation: andR non-negative. Then we have: Let S = 0 and let a 0 (t, x, p) := a 0 (x − tp, p)e −R 0 t be the ballistic part of a. Then, assuming that F k a 0 ∈ L 1 (R 2d ), we have the following estimate for all t > 0: Proof. The proof is a direct application of the integral formulation of (30), where G t is the free transport semigroup given by G t a(x, p) := a(x − tp, p).
The operators Q and G t are both continuous in Y p . Indeed, for ϕ ∈ Y p , we have: Standard fixed point techniques then provide existence and uniqueness results for (30). When S = 0, estimate (31) follows from the maximum principle and the observation that a 0 Yp is a majorizing solution to (30). When a 0 = 0, (31) is an application of the Gronwall lemma. For S = 0, we have the following Neumann series expansion in terms of multiple scattering: a n (t) = t 0 e −R 0 (t−s) G t−s Qa n−1 (s)ds, with the ballistic part a 0 (t, x, p) := e −R 0 t a 0 (x − tp, p). By induction, we find the following expression for the Fourier transform of a n : The change of variable k + tu → ξ yields Summing over n ≥ 1 gives the result. The last lemma deals with the fourth-order transport equation (11): 3 Assume a 0 ∈ X p and S ∈ L 1 ((0, T ), X p ), for T > 0 and 1 ≤ p ≤ ∞. Then, the above system admits a unique solution in C 0 ([0, T ], X p ) such that: Proof. The result stems from the integral formulation of (11) given by where G 2 t is the semigroup defined as G 2 t a(x, p, y, q) := a(x − tp, p, y − tq, q).
From Lemma 3.1, we know K η is continuous in X p , and so are G 2 t and Q 2 since for ϕ ∈ X p . Existence and uniqueness follow as before from standard fixed point theorems while estimate (34) stems from separate applications of the maximum principle and the Gronwall lemma.

Estimates for the scintillation
We are now ready to prove Theorem 2.1 and Corollary 2.2.
Proof [Theorem 2.1]. According to Lemma 3.3, the fourth-order transport equation (21) is stable in X p , so that we have the following estimate, uniformly on [0, T ], Provided that a η belongs to Y ∩ Y p , then K η a η ⊗ a η is small in X p . Indeed, item (ii) of Lemma 3.1 yields for s ∈ [0, T ] and 1 ≤ p ≤ ∞ that: First, we control the Y norm by the Y 1 norm since Y 1 ⊂ Y . Lemma 3.2 shows that the radiative transfer equation (8) is stable in Y r , for 1 ≤ r ≤ ∞, so that we just need to estimate the initial condition a η0 (x, p) := a η (0, x, p) in these Y r norms. Denoting by F x a η0 (u, p) the Fourier transform of a η0 with respect to the spatial variable x only, we obtain, so that the assumption of the theorem gives Moreover, defining ψ η0 (·) := ψ η (·, 0), we have the relation from which it follows, using the Cauchy-Schwarz inequality, that We have thus obtained that for all s ∈ [0, T ], which yields by interpolation, for 1 ≤ p ≤ ∞, This induces a first estimate for J η , which is not optimal for initial conditions with small support when α − β > 0. The stability of the transport equation (8) in Y p is not sufficient to deal with such irregular initial conditions. Rather, we need to separate the ballistic part from the scattering part in the kinetic equation to obtain sharper estimates and thus introduce: where a 0 η (t, x, p) = e −R 0 t a η0 (x − tp, p) is the ballistic part and a s η satisfies ∂a s Since the Fourier transform of a 0 η is given by e −R 0 t F a η0 (u, k+tu), its Y norm can be estimated for t ∈ (0, T ] as: Now, Lemma 3.2 and estimate (32) imply that: so that the time singularity of a s η is weaker than that of a 0 η . Thus, for 1 ≤ p ≤ ∞, For short times, we then use estimate (36) since it is independent of s and for longer times, we use the above estimate. We thus write: Setting t 0 (η) = η α−β when α > β above and using (36) and (35), we find, for We conclude by using the Parseval-Plancherel equality which yields, for t ∈ [0, T ], , It remains to verify the scaling properties: .
We conclude the proof of the theorem by choosing p = ∞ or p = 1 in the above estimates.
Proof [Corollary 2.2]. We simply need to estimate F x a 0η and F k a 0η in it follows that: It suffices to set β = 1 − α in Theorem 2.1 to conclude the proof of the corollary.

Convergence of the scintillation
We now prove the announced convergence result. We first observe that the existence results obtained in Lemmas 3.2 and 3.3 hold when the spaces Y p and X p are replaced by the spaces of bounded measures M(R 2d ) and M(R 4d ), respectively or by the spaces of continuous functions C 0 (R 2d ) and C 0 (R 4d ), respectively. We recall that d ≥ 2 here. Proof [Theorem 2.3]. The scintillation function satisfies the following transport equation in integral form (37) We recast this, with obvious notation, as We denote by T 2 the formal limit operator of T 2η defined as The source contribution. We verify that where the ballistic part is given by Indeed, we know from Lemma 3.1 that and from (32) in Lemma 3.2 that That a Y∞ is bounded comes from the stability of the transport equation in Y ∞ established in Lemma 3.2. The term K η a 0 ⊗ (a − a 0 ) is treated similarly. Let us define We find that J 0 Indeed, we deduce from (39) and the stability of Up to a smaller-order error term in the space of distributions, we may thus replace J 0 η by J 00 η in the sequel since the transport equation (37) is stable in X ∞ . Now, calculations with K η replaced by duds. Upon sending η → 0, we find in the limit that in the space of bounded measures M(R 4d ). After accounting for all four terms in the definition of K η and using the fact thatR(u) =R(−u), we find that the limit of η −1 J 0 η (t) is given by: This gives us the source term (22) in the transport equation (21).
Kinetic equation for the scintillation. We have shown that η −1 J 0 η converged to J 0 . It remains to obtain convergence of the whole sequence η −1 J η . Let φ(t, x, p, y, q) be a a smooth function on [0, T ] × R 4d . Then we have by integration on the latter space that and equivalently that with We have shown that the difference between the source terms η −1 J 0 η and J 0 converges to 0 as a distribution and has a negligible effect on η −1 J η . So we can replace the initial condition for the error term by J 0 and look at the problemJ η = T 2ηJη + J 0 , whereJ η is now of order O(1). We observe that J 0 (t) = e −2R 0 t δ(x − tp)δ(y − tq)H(p, q), where H(p, q) is a smooth function.
Let now J 1 η =J η − J 0 be the solution of We recall that T 2η J(t) = t 0 e −2R 0 (t−s) G 2 t−s (Q 2 + K η )J(s)ds, so that T 2η J 0 = T 2 J 0 + J 2 η , where J 2 η is given by a bounded operator in M(R 4d ) applied to K η J 0 . The latter is given by plus similar contributions. Because H is a smooth function, this term converges to 0 in M(R 4d ) as η → 0. This shows that J 2 η converges to 0 as η → 0.
The other contribution, T 2 J 0 , involves a bounded operator applied to Q 2 J 0 , which is equal to For f , whence H, and R sufficiently smooth, the above function is bounded in C 0 (R 4d ). The function is not bounded uniformly in time, however, and we split the contribution J 0 (t) into J 0 δ (t) = J 0 χ (0,δ) (t) and J 0 χ (δ,T ) (t), which we still denote by J 0 (t). The source term T 2 J 0 δ generates a small contribution, which goes to 0 as δ goes to 0 in the sense of distributions since the term in (45) is bounded in e.g. L 1 (R 4d ) uniformly in time so that after time integration in (38), the contribution is bounded by O(δ) → 0. The remaining contribution is bounded in the uniform norm uniformly in time with bound inversely proportional to δ 2d .
We now have a problem of the form J 1 η = T 2η J 1 η + T 2 J 0 , where T 2 J 0 is uniformly bounded in the uniform norm by O(δ −2d ). Weakly, this means that (J 1 η , φ) = (J 1 η , T * 2η φ) + (T 2 J 0 , φ), where φ(t, x, p, y, q) is a smooth function. The solution J 1 η is bounded in C 0 (R 4d ) uniformly in η by stability of the fourth-order transport equation in the uniform norm. There is therefore a subsequence that converges weak * in L ∞ (R 4d ) to a limit J 1 ∈ L ∞ (R 4d ).
The above convergence to the limiting transport equation holds for every cutoff δ. Thus, by stability of the limiting transport equation, we can remove the cut-off in δ and obtain that weakly in the space of distributions. The above integral equation admits a unique solution, which shows that the whole sequence η −1 J η converges to J solution of: This completes the proof of Theorem 2.3.