Corrections to Einstein's relation for Brownian motion in a tilted periodic potential

J. C. Latorre, G. A. Pavliotis, P. R. Kramer

math-ph cond-mat.stat-mech physics.comp-ph

Introduction

Brownian motion in periodic potentials is one of standard models in condensed matter physics. Applications include the modeling of Josephson junctions, polymer dynamics, superionic conduction, dielectric relaxation, plasma physics and surface diffusion . A detailed discussion and extensive bibliography can be found in .

The goal of this paper is to study Brownian motion in a tilted periodic potential for arbitrary values of the drift and of the tilt (external forcing). The dynamics of the Brownian particle is governed by the Langevin equation

where $V$ is a periodic potential with period $L$ (in each direction), $\boldsymbol{F}$ denotes a constant external force, so that the effective potential is

The main objective is the calculation of the drift and diffusion coefficients which are defined as

Here $\langle\cdot\rangle$ denotes ensemble average and $\otimes$ stands for the tensor product. Explicit formulas for these coefficients are available only in the overdamped limit and mostly in one dimension. An exact analytical formula for the effective velocity of an overdamped Brownian particle moving in a one dimensional tilted periodic potential was obtained many years ago by Stratonovich (,[30, Ch. 9]). A corresponding analytical formula for the diffusion coefficient was obtained and analyzed more recently , and verified in an experimental realization of the model involving rotating optical tweezers . Simpler algebraic formulas were deduced from these for the special case of piecewise linear potentials in . Only potentials with very specific geometries can lead to analytical formulas in dimension higher than one . A wealth of information on the problem of Brownian motion in a tilted periodic potential in one dimension can be found in [27, Ch. 11].

It is well known that the equilibrium diffusion coefficient (i.e., the diffusion coefficient in the absence of an external drift) and the drift, or, rather, the mobility are related through the famous Einstein formula:

The validity of this formula has been proved rigorously for several models , including that of a Brownian particle in a tilted periodic potential . Formulas of the form (5) can be understood in the more general framework of linear response theory and of the Green-Kubo formalism . A recent rigorous analysis of the Green-Kubo formalism for the calculation of the shear viscosity can be found in .

The main goal of the present paper is to investigate the validity and usefulness of corrections to linear response theory. In particular, we calculate all terms in the power series expansions (with respect to the forcing $F$ ) for the drift and diffusion coefficients and we use these in order to calculate corrections to Einstein’s formula (5). Our analysis is based on the formalism of averaging and homogenization . From this formalism we know that both drift and diffusion coefficients can be expressed in terms of the solution of appropriate Poisson equations (14)-(17). Details are presented in the next section.

For simplicity of notation and presentation, we will restrict our calculations for corrections to linear response theory to the one dimensional case $d=1$ , and hence hereafter drop vector and tensor notation. Completely analogous formulas are applicable in multiple dimensions. We present our results in detail in Section 3, and summarize them here rather imprecisely:

The drift and diffusion coefficients admit the asymptotic expansions

The coefficients $V_{j},\,\Sigma_{nj},\,\Xi_{nj};n,\,j=1,2,\dots$ can be computed in terms of solutions to Poisson equations for the generator of the equilibrium dynamics $F=0$ . In particular, the higher order corrections to the drift and diffusion coefficients are not compatible with an extension of the Einstein relation (5) beyond the linear response regime $F\rightarrow 0$ .

Thus, it is possible, in principle, to calculate the drift and diffusion coefficients of the nonequilibrium dynamics (1) in terms of the equilibrium dynamics

for at least a finite interval of values of $F$ .

The validity and usefulness of the power series expansions (6) is tested by performing numerical experiments. For the calculation of drift and diffusion coefficients we need to solve Poisson equations of the form

with $\mathcal{L}$ being the generator of the Markov process $\{q,\,p\}$ with $p=\dot{q}$ . We solve equations of the form (9) using a spectral method that is an extension of Risken’s continued fraction expansion method . By comparing the results obtained using our spectral method with results obtained from (the computationally more expensive) Monte Carlo simulations, we demonstrate that our method performs very well.

The rest of the paper is organized as follows. In Section 2 we present the formulas for the drift and diffusion coefficients obtained using homogenization theory. In Section 3 we calculate the power series expansions for the drift and the diffusion coefficient. In Section 4 we present results of numerical simulations on the calculation of $U$ and $D$ . Section 5 summarizes our conclusions. The details of the spectral method for the solution to the Poisson equation are presented in Appendix A. Some discussion of how our formulas relate to an alternative approach developed in can be found in Appendix B.

Formulas for the Drift and Diffusion Coefficients

We start by writing (1) as a first order system, in $d=1$ dimension:

The process $\{q,\,p\}$ is a Markov process with generator

The Fokker-Planck operator, i.e. the $L^{2}$ –adjoint of the generator, is

The potential function $V$ has period $L$ . We can use homogenization theory to prove that the rescaled process

The diffusion coefficient is given by the formula

where $\phi$ is the solution of the Poisson equation

Equations (14) and (17) are equipped with periodic boundary conditions in $q$ and suitable integrability conditions . Formula (16b), which shows that the effective diffusion tensor is positive semidefinite, follows from (16a) after an integration by parts.

It is possible to prove that both $U$ and $D$ are analytic functions of the forcing $F$ . This has been proved for the drift in (in fact, in this paper the analyticity of the drift with respect to the forcing is proved for several models including systems of coupled Fokker–Planck equations). A similar analysis can be used to prove the analytic dependence of $D$ on $F$ .

Corrections to Linear Response Theory

In this section we solve perturbatively equations (14) and (17) in one dimension, in order to obtain the power series expansions (6). Calculations of this form are quite standard when investigating the effect of colored noise on the drift and diffusion coefficients, e.g. . Related calculations have presented recently in .

The main result of this section is a precise formulation of Proposition 1.1. To state the result, we need to introduce some notation. We denote by $H_{0}$ the Hamiltonian of the unperturbed (equilibrium) dynamics (1):

denote the reversible and irreversible parts, respectively. The operators $\mathcal{A}$ and $\mathcal{S}$ are antisymmetric and symmetric, respectively, in ${\mathcal{H}}$ . We introduce now the creation and annihilation operators

These two operators are ${\mathcal{H}}$ -adjoint:

The drift and diffusion coefficients admit the asymptotic expansions

where $f_{j},\,\phi_{j},\;j=0,\dots$ are solutions to the (adjoint) Poisson equations

We have used the notation $\hat{\mathcal{L}}_{0}=-\mathcal{A}+\gamma\mathcal{S}$ to denote the ${\mathcal{H}}$ -adjoint of the generator $\mathcal{L}_{0}$ .

Using the notation that we have introduced in this section, Einstein’s formula (linear response theory) can be written in the form

However, formula (22a) shows that it is not true that a similar simple relation holds for higher order terms in the expansions for the drift and the diffusion coefficients. In particular, it is not true that

but instead there is a non-trivial correction to (25) that is given by the second term on the right hand side of (22a). As an example, we present the formula for the diffusion coefficient that is valid up to $\mathcal{O}(F^{3})$ :

Notice that the calculation of the next two terms in the expansion for the diffusion coefficient requires the solution of an additional Poisson equation, in order to compute $\phi_{1}$ , as well as the calculation of three additional integrals.

Similarly, it is not true that the Einstein relation (5) can be extended away from $F=0$ in the form:

because of the presence of correction terms in Eq. (22b). This issue is investigated numerically in the next section, see Figures 4 and 5. The relation Eq. (27) was indeed hypothesized in , but showed through analytical and numerical studies that while it seems qualitatively correct, and is quantitatively correct in the three limits $F\rightarrow 0$ , $F\rightarrow\infty$ , and $\beta\rightarrow\infty$ , it is not quantitatively accurate for general parameter values. Our results in Proposition 3.1 give quantitative formulas for this discrepancy, for example, through third order:

The violation of the Einstein relation for $F\neq 0$ in the model under consideration, and other nonequilibrium steady-state models, was recently analyzed by from a different nonperturbative perspective, expressing the correction terms with respect to various time-correlation functions of the dynamics. But as we discuss in Appendix B, our framework based on perturbation expansions of the equations from homogenization theory appear to yield more easily computable expressions. We remark also that have examined deviations from the Einstein relation in the context of stochastic tracer dynamics in a random environment.

Proof of Prop. 3.1. We start with the analysis of the stationary Fokker-Planck equation (14). We set

We substitute (30) into (14) and use the symmetry and antisymmetry of $\mathcal{S}$ and $\mathcal{A}$ , respectively as well as equation (20) to obtain

We look for a solution to (31) as a power series expansion in $F$ :

We substitute the expansion for $f$ into (31) to obtain the sequence of equation

The null space of $\hat{\mathcal{L}}_{0}$ , as well as its $\mathcal{H}$ -adjoint $\mathcal{L}_{0}$ is one-dimensional and consists of constants. Consequently, the solvability condition for equations of the form (35) is that

Provided that the solvability condition (36) is satisfied, the Poisson equation (35) has a unique mean zero solution, $\langle\psi,{\bf 1}\rangle_{\beta}=0$ . We correspondingly define the operator $(-\hat{\mathcal{L}}_{0})^{-1}$ on the subspace of functions satisfying (36) to be this unique mean zero solution.

From the first equation in (34) and the normalization condition we deduce that

The properties of the operators $a^{\pm}$ immediately yield that the solvability condition is satisfied for all equations for $f_{j},\,j=1,\dots$ :

The solution of Equations (34b) can be written as

where $\hat{\mathcal{K}}$ is the $\mathcal{H}$ -adjoint of $\mathcal{K}:=a^{-}\left(-\mathcal{L}_{0}^{-1}\right)$ . Consequently,

Thus, we have obtained a power series expansion for the invariant distribution in powers of $F$ :

from which we immediately deduce the expansion for the effective drift:

Now we proceed with the analysis of the Poisson equation (17) which, in view of (40), (30), (32), and (37) , can be written as

with $V_{j}$ given by (41). The generator of the perturbed dynamics is

where $\mathcal{L}_{0}$ is given by (19). We look for a solution of (42) in the form of a power series expansion in $F$ :

We substitute this expansion into Equation (42) to obtain the sequence of equations (recalling from Eq. (37) that $f_{0}=1$ ):

Equation (43a) is precisely the Poisson equation for the unperturbed dynamics $F=0$ . Now we show that the solvability condition (36) is satisfied for equations (43b). We need to show that

The solvability condition (44) is satisfied for all $j\geq 1$ , and moreover the relation

Our strategy pivots on the observation that if we can establish (45) for $m=0,1,2,\ldots,M$ , then the solvability condition (44) follows for $j=1,2,\ldots,M+1$ :

using Eq. (45) with $m=j-1$ in the penultimate equality.

Since, by the above argument and the induction hypothesis, the solvability condition for (43b) is satisfied for $j=1,2,\ldots,M+1$ , we can write

where the second sum of constants is included to meet the side condition in Eq. (17), as we have defined $(-\mathcal{L}_{0})^{-1}$ to yield a mean zero solution. But the operator $a^{-}$ will kill these constants, and therefore, for $k=M+1,M+2,\ldots$ , we can write:

Now we derive (22a). Using the centering condition in (17) we have that the diffusion coefficient is given by

We can also alternatively restructure this expansion as follows, using the relations (44) and (45):

establishing the statement (22b) in the proposition. In the last equality, we used that $a^{-}f_{0}=0$ . ∎

Numerical Simulations

In this section we present results of numerical simulations that corroborate the theoretical results presented in the previous section. The calculation of the drift and diffusion coefficients is based on the numerical solution of the hypoelliptic boundary value problems (14) and (17) as well as the calculation of the integrals (15) and (16a). Both PDEs are solved using a spectral method that relies on the expansion of the solution of the stationary Fokker-Planck and the Poisson equations in a Fourier-Hermite expansion. This method is adapted from Risken’s continued fraction expansion method ; see also . This method was used previously in the study of the diffusion coefficient for a Brownian particle in a periodic potential in . Details about the numerical method can be found in Appendix A.

In all the numerical experiments we use a cosine potential, $V(q)=V_{0}\cos(\omega_{1}q)$ , with $\omega_{1}=2\pi/L$ . As a first test for the validity of our numerical method, in Figures 1, 2 and 3 we compare the results obtained from the solution of the two PDEs with results obtained using Monte Carlo simulations. In particular, in Figures 1 and 2 we reproduce the results reported in for $\gamma=0.01$ and go beyond this for larger values of $\gamma$ . In all the Monte Carlo simulations reported in this paper we take a sufficiently large number of realizations, a sufficiently small time step and sufficiently long paths so that the results of the simulations are very accurate.In fact, in all the figures where the results of Monte Carlo simulations are presented, we also include the error bars. However, they are so small that they are barely visible. Details on the values of the parameters used in the simulations can be found in the caption figures. In Figure 3 we present results for $U$ and $D$ for larger values of $\gamma$ .

We emphasize the fact that the spectral method enables us to calculate the drift and diffusion coefficients very accurately for a very wide range of values of the friction coefficient $\gamma$ as well as the forcing $F$ . As expected, the numerical method becomes computationally more expensive as $\gamma$ decreases, since more Hermite and Fourier modes are needed for the accurate calculation of the diffusion coefficient. We note also that, in two and higher dimensions, the underdamped regime requires appropriate preconditioning for the efficient solution of the resulting linear algebraic problem.

Now we turn our attention to the numerical study of formulas (21) and (27). In Figure 4 we have calculated numerically the effective drift $U$ using (15), and we have also calculated numerically the coefficients $V_{n}$ in (21). For this we need to solve the Poisson equations (24a), where the generator of the unperturbed dynamics, i.e. with $F=0$ , appears. We can see that as we increase the number of terms in the power series expansion, the series converges to the value of $u$ computed from solving the stationary Fokker-Planck equation (14) and computing the integral in (15). We stress that, using the expansion (21) we can calculate the nonequilibrium drift for arbitrary values of the external forcing $F$ using only information from the equilibrium dynamics.

Finally, we investigate the overdamped limit. The drift and diffusion coefficients of an overdamped particle moving in a one dimensional periodic potential under constant external force can be computed analytically in terms of quadratures. The exact formula for the effective drift is computed in (,[30, Ch. 9]), whereas the exact formula for the diffusion coefficient can be found in and .

Expanding (14) and (17) in inverse powers of $\gamma$ we obtain

The generator is posed on $[0,L]$ equipped with periodic boundary conditions. Similarly, the diffusion coefficient is given by

on $[0,L]$ with periodic boundary conditions. Higher order corrections in (46) and (47) can be obtained through the solution of further auxiliary Poisson equations. As shown in Figure 7 , the overdamped formulas for the drift and diffusion coefficients offer a very accurate approximation even for moderately high values of the friction coefficient, uniformly in $F$ .

Conclusions

Using the framework of homogenization theory and multiscale analysis, we have developed a systematic expansion of the effective drift and effective diffusivity for the nonequilibrium dynamics of a particle in a tilted periodic potential. The coefficients in this expansion relate the nonequilibrium transport coefficients to statistical averages involving the equilibrium dynamics (with no imposed tilt), computed through the solutions of boundary value problems for deterministic partial differential equations of hypoelliptic type. The expansions give a detailed description of how Einstein’s relation between the diffusivity and mobility of a particle is violated in higher orders with respect to the perturbation from equilibrium. Our theoretical results were confirmed by numerical simulations based on a new efficient spectral method for the solution of Poisson equations for the generator of the Langevin dynamics.

Our method of analysis can be readily extended, with suitable elaboration of notation, to multiple dimensions. Other substantial directions for future research include the application of the homogenization procedure to multiscale and locally periodic potentials, as well as to time-dependent external forcing. This last setting could have particular relevance to the study of stochastic resonance phenomena.

Acknowledgments This work is supported by the DFG Research Center Matheon “Mathematics for Key Technologies” (FZT86) in Berlin. The research of G.A.P. is partially supported by the EPSRC, Grant No. EP/H034587. PRK wishes to thank the Zentrum für Interdisziplinäre Forschung (ZiF) for its hospitality and support during its “Stochastic Dynamics: Mathematical Theory and Applications” program, during which part of this research was completed. The research of J.C.L. was partially supported by NSF DMS-0449717.

Appendix A Numerical Algorithm

In this appendix we present a numerical approach for solving the 1-dimensional stationary Fokker-Planck equation (14) together with the cell problem (17) for computing $V$ and $D$ via (15) and (16a). This numerical method is based on the approach by and consists in a spectral decomposition of the solution of (14) and (17) in terms of Hermite polynomials and Fourier series, followed by a recursive method to solve the resulting system of algebraic equations. Since this approach is presented in for finding numerically $\rho_{\beta}(q,p)$ and $U$ , we will focus on the computation of $D$ via the solution for the auxiliary field $\phi(q,p)$ in equation ((17)) and equation (16a).

The cell problem for the auxiliary field $\phi(q,p)$ can be written in terms of the infinitesimal generator of the Ornstein-Uhlenbeck (OU) process as introduced in Section 3,

where $\phi_{n}(q)$ is a series of functions to be determined. $H_{n}(p)$ are rescaled Hermite polynomials

which are the eigenfunctions of the operator $\mathcal{S}$ ,

Also, these rescaled Hermite polynomials are orthonormal with respect to the unperturbed stationary distribution:

Upon substituting (50) into ((17)), projecting against $H_{0}$ , $H_{1}$ , and $H_{n}$ for $n\geq 2$ respectively , and using the orthonormality property of the Hermite polynomials, we obtain the following infinite system of ordinary differential equations for $\phi_{n}(q)$ ,

A.2 Spectral decomposition.

Since the solution to the cell problem must be periodic in $q$ , the auxiliary functions $\phi_{n}(q)$ must also be periodic. It is natural then to express these functions in terms of their Fourier series,

For simplicity, we will focus now on the simplest periodic potential, namely,

although more complex potentials can be studied. In terms of this potential, the equations take the following form,

We now proceed to describe the numerical algorithm for computing $D$ . In order to solve (17) in its spectral representation (A.2), we approximate $\phi_{n}(q)$ by a Galerkin truncation of the Fourier series after the $M$ th term,

The infinite system of algebraic equations (A.2) becomes then an infinite, tri-diagonal system of equations expressed as follows. By explicitly writing the real and imaginary parts of $\Phi_{n}^{j}=\xi_{n}^{j}+i\eta_{n}^{j}$ and using the fact that the solution must be real-valued (which implies that $\xi_{n}^{-j}=\xi_{n}^{j}$ and $\eta_{n}^{-j}=-\eta_{n}^{j}$ ) we form the vectors,

This representation leads to the following system of equations,

These matrices are given, for $n=0,1,\ldots$ , by,

where $\boldsymbol{I}_{k}$ is the $k$ x $k$ identity matrix. For $n=1,2,\ldots$ we have,

where $[\boldsymbol{B}]_{k}$ represents the $k$ th element of the vector $\boldsymbol{B}$ (respectively for $\boldsymbol{A}$ .) In order to solve the infinite system of algebraic equations, we impose some boundary condition of the form $\boldsymbol{\Phi}_{N+1}=\boldsymbol{S}_{N}\boldsymbol{\Phi}_{N}$ , for large $N$ . Tested boundary conditions include $\boldsymbol{S}_{N}=\boldsymbol{0}$ (Dirichlet boundary condition), which we employed in the simulations in Section 4 , and $\boldsymbol{S}_{N}=\boldsymbol{I}_{2M+1}$ (Neumann boundary condition). Defining matrices $\{\boldsymbol{S}_{n}\}_{n=0}^{N-1}$ recursively downwards from $n=N-1$ by

we can check by induction (again downwards) that for $n=1,\ldots,N$ ,

Indeed, this relation is already in force for $n=N$ , and assuming it to be true for some $n=m\geq 2$ , from Eq. (52c) we find:

so that Eq. (53) holds for $n=m-1$ as well. Turning now to Eqs. (52a) and (52b), we have

from which we find by solving for $\boldsymbol{\Phi}_{1}$ in terms of $\boldsymbol{\Phi}_{0}$ :

Substituting this expression into Eq. (54b), we finally obtain a closed equation for $\boldsymbol{\Phi}_{0}$ :

The matrix $\boldsymbol{Q}_{0}^{-}\boldsymbol{S}_{0}$ will have one null eigenvalue (corresponding to the null space of $\mathcal{L}$ ). One can verify, by considering the analogous numerical solution scheme for $\rho_{\beta}$ and $U$ , presented in , that the right hand side of Eq. (55) satisfies the solvability condition that it be orthogonal to the left eigenvector of $\boldsymbol{Q}_{0}^{-}\boldsymbol{S}_{0}$ with zero eigenvalue. A unique solution for $\boldsymbol{\Phi}_{0}$ is then obtained by discretization of the auxiliary condition in Eq. (17). In particular, representing the solution to the stationary Fokker-Planck (14) by a Hermite polynomial expansion

and approximating the functions $R_{n}(q)$ by a finite Fourier series, with coefficients organized into vectors $\boldsymbol{R}_{n}$ analogously to Eq. (50), this auxiliary condition reads:

This then determines, with Eq. (55), $\boldsymbol{\Phi}_{0}$ from which he remaining $\{\boldsymbol{\Phi}_{n}\}_{n=1}^{N}$ are found recursively using the matrices $\boldsymbol{S}_{n}$ and the relations (53). Once the vectors $\boldsymbol{\Phi}_{i}$ are found, $D$ is easily computed by replacing the proposed solution for $\rho_{\beta}(q,p)$ and $\phi(q,p)$ in (16a) and using the Hermite polynomial properties to obtain:

Appendix B Alternative Approach to Obtaining Corrections to Einstein’s Formula

The relation between the diffusivity and mobility is expressed in as follows (in our notation):

where $\langle g;h\rangle=\langle gh\rangle-\langle g\rangle\langle h\rangle$ and $\langle\cdot\rangle$ denotes an average over the stochastic noise (and possibly random initial conditions). The correction term was studied in on the model system (1) as well as other non-equilibrium systems through direct numerical simulation of the governing dynamical equations and Monte Carlo estimation of the statistical average. We can express Eq. (56) in terms of deterministic operators through the following formal manipulations. First, we re-express

which avoids the complication of working with the nonperiodic variable $q(t)$ . We then have:

Now, thanks to the large factor of $T$ in the denominator, we may neglect initial transients and evaluate the statistical average in the nonequilibrium steady state, i.e., with single-time statistics governed by the invariant density $\rho_{\beta}$ , the solution of the stationary Fokker-Planck equation (14). We then express the two-time correlation function formally using the evolution operator $e^{\mathcal{L}\Delta t}$ , where $\mathcal{L}$ denotes the generator of the Langevin dynamics, for the forward-in-time variable, and the projection operator

A somewhat more compact formula can be obtained by defining the $L^{2}(\rho)$ adjoint of $\mathcal{L}$ , which can be computed as:

where $\hat{\mathcal{L}}_{0}$ is defined at the end of Proposition 3.1, so that

Inspecting expression Eq. (57) for the correction to the Einstein relation, we see that beyond computing $\rho_{\beta}$ as the stationary solution of the Fokker-Planck equation (14), we must solve a Poisson equation of the form (17) as well as a second Poisson equation of the form

In the expression Eq. (59), we must solve a stationary Fokker-Planck equation (14), a Poisson equation of the form (17), as well as an adjoint-Poisson equation of the form

In both cases, it seems that an additional equation would need to be solved beyond the stationary Fokker-Planck equation (14) and a single Poisson equation (17) necessary in the homogenization approach. On the other hand, computing the mobility $\frac{dU}{dF}$ at general values of the tilt $F$ from the nonperturbative homogenization equations would require a differentiation between different values of $F$ . The direct formula (56) would generally of course need to be evaluated through Monte Carlo averages involving a large number of sample trajectories run for sufficiently long.

The perturbation theory with respect to $F$ developed for the homogenization equations in Section 3 has the virtue of allowing the simultaneous numerical computation of the diffusivity and drift for a range of values of tilt $F$ , rather than one value at a time. One could introduce similar perturbation expansions with respect to tilt $F$ into the formulas (57) and (59). We attempted to examine whether this would give equivalent results, but found this effort frustrating. On the one hand, computing Eq. (57) perturbatively would introduce the perturbative series solution to a second Poisson equation completely absent from the homogenization theory, so it would be difficult to relate the results. Expression (59) has more promise because to leading order, $\hat{\mathcal{L}}^{-1}$ is identical to the simple operator $\hat{\mathcal{L}}_{0}^{-1}$ , which is just a time reversal of the operator $\mathcal{L}_{0}^{-1}$ . However, implementing the perturbation expansion on Eq. (59), even to first order, produced considerably more unwieldy equations than emerged from the homogenization equations, and again how to relate the resulting expressions was unclear. The main complication is the propagation of the perturbation expansion (32) for the invariant density through the adjoint operator $\hat{\mathcal{L}}$ (58). Perhaps a more clever analysis would provide a linkage between the formula for the correction (56) to the Einstein relation from and the perturbative expansion we have developed in Proposition 3.1, but it appears that computations are considerably simpler by conducting the perturbation expansion on the homogenization equations as we have done in Section 3.