Perturbation and Numerical Methods for Computing the Minimal Average Energy

We investigate the differentiability of minimal average energy associated to the functionals $S_\ep (u) = \int_{\mathbb{R}^d} (1/2)|\nabla u|^2 + \ep V(x,u)\, dx$, using numerical and perturbative methods. We use the Sobolev gradient descent method as a numerical tool to compute solutions of the Euler-Lagrange equations with some periodicity conditions; this is the cell problem in homogenization. We use these solutions to determine the average minimal energy as a function of the slope. We also obtain a representation of the solutions to the Euler-Lagrange equations as a Lindstedt series in the perturbation parameter $\ep$, and use this to confirm our numerical results. Additionally, we prove convergence of the Lindstedt series.


Introduction
Let d ∈ N be fixed, and let V : R d × R → R, be periodic under integer translations. That is V(x + k, y + l) = V(x, y) for all (k, l) ∈ Z d × Z, where (x, y) = (x 1 , . . . , x d , y) ∈ R d × R. Furthermore, assume V is analytic. We consider the formal variational problem where ε is a small parameter, so that S ε is a small perturbation of S 0 (u) = R d 1 2 |∇u| 2 dx. We call u ∈ H 1 loc (R d ) a minimizer (or minimal solution) of S ε if for all φ ∈ H 1 comp (R d ) Any minimizer of (1) must solve the Euler-Lagrange equation (2) − ∆u + εV y (x, u) = 0, and by standard elliptic regularity theory will be at least as regular as V, so analytic in our case.
Definition 1.1. Following [Mos86], we say a continuous function u is non-selfintersecting if the graph of u does not intersect integer translates of itself. That is, if u ∈ C 0 (R d , R) and ∀(k, l) ∈ Z d × Z (3) u(x + k) + l − u(x) > 0, or < 0, or ≡ 0, where the three alternatives are independent of x.
This is also referred to in the literature as the Birkhoff property. From a geometric viewpoint, this means that the graph of u projects into T d+1 = R d+1 /Z d+1 without intersecting itself, unless it coincides exactly.
If u is non-selfintersecting, then there is a rotation vector, ω ∈ R d , associated to u such that Mos86]. A function u satisfying (4) is called plane-like because its graph is at a bounded distance from the a hyperplane in R d+1 with normal vector (ω, −1). We denote the set of all non-selfintersecting minimizers of S ε with rotation vector ω by M ω (S ε ), and more briefly as M ω , where dependence on S ε is understood. As shown in [Mos86], M ω is nonempty for all ω ∈ R d .
The rational dependency of ω can play a role in the structure of M ω . The case wherē ω = (ω, −1) is rationally independent was studied in [Ban87b], and the rationally dependent case in [Ban89].
We define the minimal average energy A ε : R d → R by where u ∈ M ω , B R = {x ∈ R d : |x| ≤ R}. It was shown in [Sen91] that A ε is well-defined, that is, the limit exists and is independent of the choice of minimizer u ∈ M ω . Furthermore, A ε is convex, so that one-sided derivatives of A ε exist at each ω ∈ R d , [Sen91]. In fact, for ε = 0, any u ∈ M ω has the form u(x) = ω · x + α for some α ∈ R (see [Mos86]). Thus, A 0 (ω) = 1 2 |ω| 2 is smooth. However, for typical V the differentiability of A ε breaks down when ε > 0, and for large enough ε the set of points where A ε is not differentiable will be dense in R d (see [Ban87a]). Because A ε has one-sided derivatives the graph of A ε will have "corners" (i.e. the jump in the gradient of A ε in a direction e j ) at points of nondifferentiability. In [Sen95], page 356 there is a formula for the one-sided directional derivative of A ε involving special types of minimizers of S ε , described in Section 3.
In this paper, we present a numerical approach to computing solutions of (2) via a gradient descent method known as the Sobolev gradient, as explained in Section 2. We use this to compute A ε and D e j A ε (ω) + D −e j A ε (ω), that is the jump in the gradient of A ε in the direction e j . We will sometimes refer to a jump in the gradient as a "corner". In Section 3 we present a perturbation method for finding the minimizers of S ε needed to apply the formula provided in [Sen95] to compute D e j A ε (ω) + D −e j A ε (ω).
Throughout the paper, we will deal with the case ω ∈ 1 N Z d . There are two good reasons for this. The first is that minimizers of (1) with irrational rotation vectors can be obtained as limits of sequences of minimizers with rational rotation vectors. The second reason is that for rational ω, the partial differential equations we need to solve are well-defined on the torus T d . That is, they have periodic boundary conditions. This allows the use of Fourier transforms in the numerical method. It also makes each equation in (16) of the form ∆φ = g, which can be solved for φ provided g has average zero.
In this setting, for a fixed ω ∈ 1 N Z d , the formal functional in (1) can be replaced by the reduced functional where z is NZ d -periodic. Considering the limits of solutions u N (x) = ω N x + z N (x), ω N ∈ 1 N Z d as N → ∞, we see that this problem is related to periodic homogenization. Rescaling so that (6) is defined on the unit cube, the process of letting N → ∞ is the same as letting the x dependence oscillate very rapidly.

Numerics
In this section we investigate numerically the size of the corners of A ε as ε varies. Since our problem comes from a variational principle, steepest descent methods are a natural approach. We consider a fixed ω ∈ 1 N Z d and seek solutions of (2) of the form u(x) = ω · x + z(x), where z is N-periodic. To find such a function u, we solve for z, and set u(x) = ω· x+z(x). In this setting, we see that z is a critical point of the reduced variational problem We recall that the gradient of S ε,N with respect to a Hilbert space, H, is the unique element of h ∈ H such that DS ε,N (z)η = g, η H for all η ∈ H. For instance, the L 2 -gradient of S ε,N at z is the unique element of L 2 = H 0 , which we will write as Considering gradients with respect to other inner products has been a fruitful endeavor (see [Neu10]). In our particular case, by considering the gradient of S ε,N in H β for β ∈ (0, 1] we avoid the stiffness of the problem that appears in the H 0 case. This is explained in more detail in the remark following the derivation of ∇ β S ε,N (z).
The standard inner product on For γ > 0 we define the inner product u, v β = (γI − ∆) β u, v H 0 , which determines a norm on H β that is equivalent to the standard norm. For more details on this see the introduction of [BLV11]. In that paper, the authors show that the descent equation ∂ t u = −∇ β S ε,N (u) satisfies a comparison principle, and preserves the class of non-selfintersecting functions. In this way, they obtain critical points of S ε,N with rotation vector ω by choosing initial condition u(x, 0) = ω · x, equivalently, z(x, 0) = u(x, 0) − ω · x = 0.
We note that different choices of γ result in different inner products, and therefore different gradients of S ε,N . Thus, ∇ β S ε,N (u) depends not only on the choice of Hilbert space, but also on the choice of inner product on that space. We calculate the H β -gradient as follows: Thus, our steepest descent equation in H β , ∂ t z = −∇ β S ε,N (z), becomes z : [0, N] d → R with periodic boundary conditions. If the solution z(x, t) of (8) approaches a critical point of S ε,N , that is z(x, t) → z c (x) as t → ∞, then z c will solve (γ + ∆) 1−β z c = (γ + ∆) −β (γz c − V y (x, ω · x + z c )), which reduces to −∆z c + εV y (x, ω · x + z c ) = 0. Then u c = ω · x + z c will solve (2) and sup x |u(x) − ω · x| < ∞.
Remark 1. The stiffness of (8) is significantly reduced for values of β ≈ 1. For instance, the L 2 gradient descent equation would be ∂ t u = ∆u − εV y (x, u). If we set G(k) = F {−εV y (x, u)}(k), then the descent equation in the Fourier domain becomes the system ∂ tûk = |k| 2û k + G(k) for each k. As |k| grows, these equations become increasingly stiff (the eigenvalues of the linearization have a large spread). However, the H β gradient descent equation in the Fourier domain becomes ∂ tûk = (γ + |k| 2 ) 1−βû k + (γ + |k| 2 ) −β G(k). So for β ≈ 1 the stiffness is greatly reduced.
2.1. Implementation of Sobolev gradient descent. We fix ω ∈ 1 N Z d and seek z : [0, ∞) × [0, N] d → R solving (7) with periodic boundary conditions. Thus, (7) can be rewritten in the frequency domain via the Fourier transform. The main benefit of this is that the pseudo-differential operators in (7) simplify greatly in the frequency domain. However, the composition V y (x, ω · x + z(x)) is complicated by the Fourier transform. We take advantage of the simplicity of the operators in the frequency domain in our numerical scheme, and pay the price of computing an inverse fast Fourier transform at each time-step. This is explained in the remark following equation (11).
We now develop the details of the implementation for d = 2. We know that z(x, t) = u(x, t) − w · x will be N−periodic if we have Nω ∈ Z 2 and if z(x, 0) = u(x, 0) − ω · x is an N−periodic function. So we set z 0 (x) = z(x, 0) = 0, and consider the equation satisfied bŷ z, the Fourier transform of z. We get (9)ẑ t (t, ξ) = −(γ + |ξ| 2 ) 1−βẑ (t, ξ) where we have also used F [·] to denote the Fourier transform. We now choose a time step ∆t and set t n = n∆t. We also break up the domain [0, N] 2 into m 2 discrete points and represent z(t n , x) as an m × m array z n (i, j). So, representing the discrete fast Fourier transform as fft[·] (and withẑ n = fft[z n ]) equation (3) is approximated as This is a quasi-implicit method, since nearly all of the right-hand-side is evaluated at the later time t n+1 . Only the nonlinear term is evaluated at time t n , so we can easily solve for z n+1 : After the n-th step, we haveẑ n . In order to compute V y (x, ω· x+z n ) we apply the inverse fast Fourier transform to get z n = ifft[ẑ n ], which requires order m 2 log(m) operations because z n is an m × m array (or a vector of length m 2 ). With z n computed, we evaluate V y (x, ω · x + z n ), requiring m 2 operations, and then we transform it with FFT, again costing O(m 2 log(m)) operations.
Ifẑ n and fft[V y (x, ω · x + z n )] were represented as vectors of length m 2 , then this same procedure would amount to multiplying fft[V y (x, ω · x + z n )] by the m 2 × m 2 diagonal matrix representing ∆t(γ + ξ 2 1 + ξ 2 2 ) −β , addingẑ n , and then multiplying the result by the m 2 × m 2 diagonal matrix representing (1 + ∆t(γ + ξ 2 This method would also work in a setting where the gradient part of the energy functional S ε is replaced by the fractional laplacian. We could use all of the same techniques above for solving the gradient descent for The critical points of (12) will solve the Euler-Lagrange equation Using the metric on H βδ given by v, Because z is periodic and the operator (−∆) δ is diagonal in the Fourier coefficients, implementing (13) numerically is the same as described above, except with powers of |ξ| 2δ in place of |ξ| 2 .
3. Perturbation method for computing minimizers and the jumps in DA ε (ω) 3.1. Foliations and laminations of minimizers.
For each ω ∈ R d the set M(ω) is closed and totally ordered. The closedness is a classical argument (see [Mor73]). The total order, meaning that if u, v ∈ M(ω) then either u > v or u < v or u ≡ v, is a consequence of the maximum principle for elliptic partial differential equations. Let x = (x 1 , . . . , x d ), and we write (x, x d+1 ) ∈ R d+1 . Let ω ∈ R d , then the total order of M(ω) means that for a given (x, x d+1 ) ∈ R d+1 there is at most one u ∈ M(ω) such that x d+1 = u(x). That is, each point in R d+1 belongs to the graph of at most one u ∈ M(ω). If for all (x, x d+1 ) ∈ R d+1 there exists u ∈ M(ω) with x d+1 = u(x), then we say M(ω) is a foliation of R d+1 . Because of the non-selfintersection property (3), such a foliation projects into T d+1 . It can happen that there are points (x, x d+1 ) ∈ R d+1 for which there does not exist any u ∈ M(ω) with x d+1 = u(x). In this case we say M(ω) is a lamination of R d+1 (and projects to a lamination of T d+1 ). For this reason, a lamination is sometimes referred to as a "foliation with gaps".
Ifω is rationally dependent, and if M(ω) defines a lamination, then there are minimizers u ∈ M ω whose graphs lie in the gaps of M(ω), [Ban89]. In addition, if we choose a direction then there is a u ∈ M ω such that u is asymptotic to some u + ∈ M(ω) in the direction β and asymptotic to some u − ∈ M(ω) in the direction −β. For details of the asymptotic behavior, see [Sen95], page 350. Such a u is said to be heteroclinic in the direction β. A formula for the one-sided directional derivative of A ε at a point ω is given on page 356 of [Sen95], and involves integrating the action over the gaps defined by the elements of M(ω) and the heteroclinics between them. We will use perturbation methods to calculate the gap borders and the heteroclinics lying in the gaps.
3.2. Lindstedt series for solutions. We seek plane-like solutions u ε (x) of that can be expanded as u ε (x) = u 0 + εu 1 + ε 2 u 2 + . . .. The series j≥0 ε j u j shall be referred to as the Lindstedt series for the solution u ε . Substituting the series into equation (14) and matching powers we arrive at the following equations for each order of ε ε 0 : ∆u 0 = 0 In order that u ε be a plane-like solution, we require u 0 to be affine and u j be periodic for each j > 1. Using the notation [·] j to refer to the j-th coefficient of the power series in ε, we write We will also write u < j for the first j terms in the Lindstedt series: u < j = u 0 + . . . + ε j−1 u j−1 . The j-th order equation in the list (15) has the form The zeroth order equation is satisfied by any affine function, so we take u 0 = ω · x + α, and at this point we are free to choose α as we like. To solve the j-th order equation, we must have that This compatibility condition is what forces specific choices of α once ω has been fixed.
We also note that the solution of (16) is determined only up to an additive constant, which will be chosen so that the equation of the following step has a solution. That is, the average of u j is chosen so that equation for u j+1 is solvable.
3.3. Existence of the Lindstedt series to all orders. We consider ω ∈ 1 N Z d fixed, and seek u ε solving (14) such that u ε (x) − ω · x ∈ L ∞ (NT d ). We must have u 0 (x) = ω · x + α to solve ∆u 0 = 0. The choice of α is free, and each value of α ∈ [0, 1) will result in a different set of equations for the u j with j ≥ 1. In this section, we show that there are at least two choices of α ∈ [0, 1) such that (16) has a solution for each j ≥ 0.
For any such choice of α, NT d V y (x, u 0 ) dx = 0, and therefore there exists a family of periodic solutions, u 1 (x) = u * 1 (x) + λ, of ∆u 1 = V y (x, u 0 ), differing only by an additive constant. We will write u * 1 for the member of the family with average zero.
has a solution for all j ≥ 1.
3.4. Convergence of the Lindstedt series. We use a Newton method to produce a sequence of functions U n (x, ε) that are analytic in ε, and converge uniformly to a solution of −∆u ε + εV y (x, u ε ) = 0 for ε ∈ B δ (0) ⊂ C for small enough δ > 0. Thus, we produce an ε-analytic function, u ε , solving (14) for ε ∈ B δ (0), so the Taylor series of u ε must coincide with the Lindstedt series, proving the convergence of the Lindstedt series. For similar convergence results, but for the case of Diophantine frequencies, see [CF96]. Lemma 3.3. Let u j ∈ H m+2 (T d ) and ε ∈ C and define u(ε, x) = j∈N ε j u j (x), where the series is convergent in H m+2 for |ε| < r, for some r > 0 and m > d/2. Define F : C → H m (T d ) by F(ε; u) = ∆u(ε, x) + εV(x, u(ε, x)) Then the derivative of F with respect to ε is From the Gagliardo-Nirenberg inequality [Nir59], we have that (D ε u(ε, x)) 2 ∈ H m because m > d/2, and that there is a constant depending on m and d, such that (D ε u(ε, x)) 2 H m ≤ C(m, d) D ε u(ε, x) 2 H m . Hence, Thus, F is differentiable and D ε F(ε; u) = G(ε; u) We define F ε (u) = −∆u + εV y (x, u). It was shown in the Section 3.3 that for any fixed M ∈ N, we can solve the first M equations from (15) for u <M ε , and that u <M ε will solve (14) up to order ε M . That is, for a sufficiently large M to be determined later, and for j ≥ 1 we define Let m > d/2 be fixed, and note that for any M > 0, u <M ε ∈ H m by the regularity theory for elliptic PDEs. If U n (x, ε) is analytic in ε and is H m+2 in x, then F ε (U n ) = −∆U n +εV y (x, U n ) is analytic in ε and H m in x by the result in Lemma 3.3. We have DF ε (U n )η = −∆η + εV yy (x, U n )η, and we need to consider carefully the behavior of DF ε (U n ) −1 . To simplify notation, we define L n ε : H m+2 (NT d ) → H m (NT d ) as L n ε = DF ε (U n ) = −∆ + εV yy (x, U n ). L n ε is a small perturbation of −∆ : H m+2 → H m . Now, −∆ maps the codimension one subspace H m+2 /R of its domain to the codimension one subspace H m /R of its range in a bounded, invertible way. But it has the simple eigenvalue λ = 0, with eigenspace spanned by the constant functions.
Lemma 3.4. Let P 0 : H m → H m /R be the orthogonal projection onto H m functions with zero average. That is, if f ∈ H m , then P 0 f = f − NT d f dx. Let −∆ 0 : H m+2 /R → H m /R be the restriction of the laplacian, and let Y be the image of H m+2 /R under −∆ 0 + εV yy (x, U n ). Suppose there are constants q, c 1 , δ 0 > 0 such that Proof. Let Q = −∆ 0 + P 0 εV yy (x, U n ), then Q : H m+2 /R → H m /R, and for small ε will have a bounded inverse, Q −1 : H m /R → H m+2 /R. We have that −∆ 0 + εV yy (x, U n ) = −∆ 0 + P 0 εV yy (x, U n ) + P ⊥ 0 εV yy (x, U n ) maps H m+2 /R into H m and is a small perturbation of Q, and there is a c 2 > 0 such that − ∆ 0 η + εV yy (x, U n )η H m ≥ c 2 η H m+2 for all η ∈ H m+2 /R. Thus, Y is a codimension one linear space isomorphic to H m /R lying inside H m , and Y ⊥ = {g ∈ H m : NT d (−∆ 0 η+εV yy (x, U n )η)g dx = 0, ∀η ∈ H m+2 /R} is one dimensional.
The condition (19) implies that the image of constant functions under L n ε does not lie entirely in Y, but has a component in the Y ⊥ direction. Thus the image of H m+2 under L n ε is H m , and for small ε, L n ε will be invertible. Let P Y and P ⊥ Y denote the orthogonal projections of H m onto Y and Y ⊥ . We let g ∈ Y ⊥ be a unit vector such that P ⊥ Y ξ = ξ, g H m g for ξ ∈ H m , which is possible because Y ⊥ is one dimensional. Let ξ ∈ H m , so ξ = L n ε η for some η ∈ H m+2 . We write η = η 1 + η 0 with η 1 ∈ H m+2 /R, and η 0 ∈ R, and ξ = ξ Y + ξ ⊥ with ξ y ∈ Y and ξ ⊥ ∈ Y ⊥ . Writing ξ in terms of η 1 , η 0 we have ∆η + εV yy (x, U n )η = ∆ 0 η 1 + εV yy (x, U n )η 1 + εV yy (x, U n )η 0 = ∆ 0 η 1 + εV yy (x, U n )η 1 + P Y εV yy (x, U n )η 0 + P ⊥ Y εV yy (x, U n )η 0 . The term ∆ 0 η 1 + εV yy (x, U n )η 1 + P Y εV yy (x, U n )η 0 ∈ Y, so the component of L n ε η in Y ⊥ is ξ ⊥ = L n ε η, g H m g = εV yy (x, U n )η 0 , g H m g. Thus, η 0 = ξ ⊥ ( εV yy (x, U n ), g H m ) −1 and η 1 is given by Hence, Recall that (L n ε ) −1 is a compact operator by the regularity theory for elliptic PDE. In particular, the eigenvalues of L n ε are isolated from the spectrum, and if λ n is an eigenvalue of L n ε , then the resolvent of L n ε , written as R(L n ε , ζ) = (L n ε − ζI) −1 , has the representation where Q n is bounded, and P n is the spectral projection on the the λ n eigenspace: and Γ is a closed curve enclosing λ n but no other point of the spectrum (see [Kat66]). The principal eigenvalue, λ 0 (ε), of L n ε is simple because L n ε is an elliptic operator, [Eva98]. This means that λ 0 (ε) is analytic in a neighborhood of ε = 0, [Kat66]. In the iteration process, (L n ε ) −1 = R(L n ε , ζ) = ∞ j=0 λ j n Q j+1 n + 1 λ n P n , will act on F ε (U n ). At each step, we need the function U n to be analytic in ε, so we need (L n ε ) −1 F ε (U n ) to be analytic. Proposition 1. Assuming condition (19) holds, there is a choice of M ∈ N such that, if U 0 = u <M , then U n will be analytic in ε for all n ≥ 0, and the Newton method (18) converges uniformly in ε in a neighborhood of ε = 0.
Proof. In the iteration process, (L n ε ) −1 = R(L n ε , ζ) = ∞ j=0 λ j n Q j+1 n + 1 λ n P n , will act on F ε (U n ). As explained in the previous paragraph, λ 0 (ε) is analytic in a neighborhood of ε = 0, and will have a zero of order p ∈ N at ε = 0. From the result of Lemma 3.4 we may assume p < q.
We will use induction on n to show the analyticity of U n+1 . F ε (u <M ) has a zero of order M at ε = 0 by construction of u <M , and we take M > 2q. By the expansion in (20) and the result from Lemma 3.3, for ζ = 0 and λ n = λ 0 (ε), (L 0 ε is a small perturbation of L 0 ε , and λ 1 (ε) will have a zero of order p. Thus, we have established the first step of the induction on n, because U 1 is analytic in ε, F ε (U 1 ) has a zero of order M 1 > 2q, and the principal eigenvalue has a zero of order p < q at ε = 0. Now assume that U n is analytic in ε and F ε (U n ) has a zero of order M n > 2q at (21) we have that U n+1 is analytic, and F ε (U n+1 ) has a zero of order at least 2q at ε = 0. Just as for the n = 0 case, U n+1 − U n H m ≤ c|ε| q so λ n+1 (ε) has a zero of order p < q at ε = 0, and the induction is complete.
3.5. Connecting orbits and corners of the energy. The connecting orbits, or heteroclinic orbits, exist only when M(ω) fails to produce a foliation of T d+1 , as described at the beginning of Section 3. Each u ∈ M(ω) satisfies ∆u = εV y (x, u) because it is a minimizer of S ε . The order-zero approximation to u has the form u 0 = ω · x + α, for α ∈ [0, 1). This is a continuous family, whose graphs foliate T d+1 .
From Lemma 3.1 we know that if V satisfies the twist condition (17) then there are at least two choices of α ∈ [0, 1) such that the Lindstedt equations (15) can be solved to all orders. Recall that the function Φ 1 , defined in Lemma 3.1, has at least two zeros in [0, 1). If the zero set of Φ 1 is not all of [0, 1), then a choice of α must be made in order to solve ∆u 1 = V y (x, ω · x + α) for periodic u 1 . Thus, the order-ε approximations, given by u 0 + εu 1 are not a continuous family indexed by α ∈ [0, 1), but rather by a strict subset of [0, 1). The graphs of the functions in t his family no longer foliate T d+1 .
If at the j-th step we find that Φ j (α) 0 for some value of α then we make the ansatz that the heteroclinic orbit will be of the form u h (x) = ω · x + α(ε j/2 x) + εu 1 (x) + . . .. And α(ε j/2 x) → α ± as x → ±∞ where α ± satisfy Φ j (α ± ) = 0. Then ∆α = O(ε j ) and the j-th order equation will be ∆α + ∆u j = [V y (x, u < j )] j−1 . Thus, α will solve a PDE of the form ∆α = f (x, α) with boundary conditions α(ε j/2 x) → α ± as x → ±∞, where That is, f is the term that keeps [V y (x, u < j )] j−1 from having a zero average for any α. This is carried out in detail in Section 4 for a specific example.
3.6. Energies involving the fractional Laplacian. The existence of minimizers in the case ω ∈ 1 N Z d has been established for energy functionals involving fractional powers of the laplacian, [LV09], [BLV11]. As described in the remark at the end of Section 2.1, the Euler-Lagrange equation for the functional Much of what has been described so far regarding solutions of ∆u = V y (x, u) carries over to this case. However, the analogous properties of the associated minimal average energy A δ ε that are presented for δ = 1 in [Sen91] and [Sen95] need to be established if one desires to investigate the differentiability properties of A δ ε in this case of non-local energy. This remains an interesting challenge.

An example of the Lindstedt series
Consider the potential V : R d+1 → R given by V(x, y) = sin(2πk · x) cos(2πy), with k ∈ Z d . For a fixed ω ∈ R d , We compute the jumps in the derivative of the average energy functional, A ε (ω) defined in (5), using asymptotic expansions of the connecting orbits.
Let u ε (x) solve ∆u ε = εV(x, u ε ). Writing u ε = j ε j u ( j) we see that ∆u 0 = 0, and thus u 0 = ω · x + α, ω ∈ R d , α ∈ R. We are considering ω fixed, but the choice of α is free and we have a family of solutions parametrized by α ∈ R.
However, in order to calculate u 1 we must solve ∆u 1 = V y (x, u 0 ) for a periodic function u 1 . The average of V y (x, u 0 ) depends on ω, so in this example we will choose ω = k so that the average of V y (x, u 0 ) is nonzero for some choices of α. Thus, the requirement that V y (x, u 0 ) = −2π sin(2πk · x) sin(2πk · x + 2πα) has zero average forces cos(2πα) = 0. We say that u ε has a first order resonance if α is restricted in solving for the order ε 1 term in the Lindstedt series.
We can calculate the first two terms in the Lindstedt series for each admissible value of α = 1/4, 3/4. We have ∆u 1 = ∓π sin(4πkx) so that the two solutions are: The superscripts c and m indicate that the solution for α = 1/4 is a minimizer of the energy functional and the solution for α = 3/4 is a critical point, but not a minimizer. We restrict our attention to minimizers, where A ε can be considered as a function of the rotation vector ω.
To compute the jump in the gradient of A ε at ω = k using the formula in [Sen95] we need to find minimizers of A ε that are heteroclinic between u m ε and u m ε + 1. We search for an asymptotic solution to ∆u = εV y (x, u) that has the form u h ε (x) = kx+α( √ εx)+εu 1 (x)+o(ε). The order ε 1 equation is We want to choose α so that it is essentially one-dimensional (i.e. for some η ∈ R d , we want α to be a function of η · x). It should also satisfy ∆α = −π cos(2πα) and we can eliminate those two terms from the equation above. At this point, any choice of η is reasonable, and in the end we will chose η to be in the direction in which we want to differentiate A ε (ω) (thus, η will typically be a standard basis vector).
Lettingη denote 1 |η| η, we have Recall that in the expression for u, α is evaluated at z = √ εx. Expanding sin(2πα( √ ε x)) and cos(2πα( √ ε x)) in Taylor series, for this choice of α, we have Then equation (22) becomes ∆u 1 = π sin(4πkx) + O( √ ε ) and we get for u h ε the expression: To be precise, we should write u h ε,η to show the dependence of u h ε on η. It is important to notice that u h ε,η is heteroclinic from u m ε + 1 to u m ε in the direction −η. Similarly, u h ε,−η is heteroclinic from u m ε + 1 to u m ε in the direction η, and from u m ε to u m ε + 1 in the direction −η.

4.1.
Computing the gradient of A ε (ω). Here we restrict η to be a standard basis vector e j . To ease notation we will write M(x) = u m ε (x) and H e j (x) = u h ε,e j (x). A formula for the derivative D e j A ε (ω) is found in [Sen95], and in our case becomes Thus, the derivative of A ε (ω) in the direction e j is difference in the energies of the minimizer defining the top of each connected component of gap and the minimizing heteroclinic in the direction e j . Note that |∇(M + 1)| 2 + εV(x, M + 1) = |∇M| 2 + εV(x, M).

4.2.
Two dimensional case. We now focus on the two dimensional case so that we can compare with the numerical computations, which were done for the two dimensional case. The computations in higher dimensions are impractical for us at the moment. In our example, ∇H η = k + η √ 2ε/ cosh( √ 2επηx) − kε cos(4πkx)/4|k| 2 . If we select η = e 1 and compute D e 1 S (k) we have The heteroclinic solution H −e 1 has the same expression as above, with a sign change on the second term. The |k| 2 terms from |∇H ±e 1 | and |∇M| = |k| 2 + O(ε) will cancel when computing D ±e 1 A ε (ω). When computing the jump in the derivative, that is the sum D e 1 A ε (ω) + D −e 1 A ε (ω), the terms ± k 1 √ 2ε cosh( √ 2ε πx 1 ) will cancel, and we are left with two terms of the form ε cosh 2 ( √ 2ε πx 1 ) . These two terms then contribute for a total of is not easy analytically, but numerically we find that they also yield 2 √ 2 π ε 1/2 . Thus, the jump in the derivative is D e 1 A ε (k) + D −e 1 A ε (k) = 4 √ 2 π ε 1/2 ≈ 1.8ε 1/2 , which agrees well with the numerical computations.

4.3.
Comparison with results from the numerical computations. We use the Sobolev gradient method to compute the minimizers u ε for several values of ε, where each u ε has rotation vector ω. We then repeat this computation for the same values of ε, but for rotation vectors ω ± ∆ωe j . With these minimizers we can compute A ε (ω) for each value of ε and each rotation vector. Taking differences over the ω-variable, we can approximate the derivative of A ε (ω) with respect to ω.
The plots in Figures 1 and 2 are for three cases of potential function. In each case we examine a first-order resonance, so the choice of α in the Lindstedt series is made when solving the order ε equation. This also means the twist condition (17) is satisfied.  Figure 1. Logarithm of the jump in D e 1 A ε (ω) against log(ε) in the two dimensional example: V(x, u) = ε sin(2πkx) cos(2πu), where k = (2, 3). This is a "first order resonance" so ω = k. 256 modes were used for each Fourier frequency direction ξ 1 , ξ 2 . V(x, u) = ε sin(2πk 1 x 1 ) sin(2πk 2 x 2 ) cos(2πu) on the left and V(x, u) = ε 2 sin(2πkx)(cos(2πu) + sin(2πu)) on the right. k = (2, 1) in each case, and ω = k. 256 modes were used for each Fourier frequency direction ξ 1 , ξ 2 . Figure 1 is log-log plot of the jump in the e 1 direction of the gradient of A ε (ω) versus ε. The dotted line is the logarithms of the numerically computed points plotted against log(ε). The solid line is the fit J(ε) = log( 4 √ 2 π ε 1/2 ), which was derived in Section 4.2. The same relation is found for the resonance at ω = −k. Figure 2 has the same type of plot for different choices of potential function. In 2(a) the potential is V(x, u) = ε sin(2πk 1 x 1 ) sin(2πk 2 x 2 ) cos(2πu), and the choice of resonance was ω = k = (2, 1). In this case, the fit is J(ε) = log( 4 π ε 1/2 ). The same relation is found at ω = (±k 1 , ±k 2 ) and at ω = (±k 2 , ±k 1 ).
Higher order resonances can be computed, for instance at ω = 2k in each of the previous examples. However, the behavior of the jump in DA ε (ω) behaves like ε to a power greater than one. To get reasonable accuracy one needs to take many more Fourier modes in the numerical approximation, which is impractical for us at the moment in the case of two spatial dimensions.