CONSTRUCTION OF RESPONSE FUNCTIONS IN FORCED STRONGLY DISSIPATIVE SYSTEMS

We study the existence of quasi--periodic solutions of 
the equation 
\[ 
e \ddot x + \dot x + e g(x) = e f(\omega t)\ , 
\] 
where $x: \mathbb{R} \rightarrow \mathbb{R}$ is 
the unknown 
and we are 
given $g:\mathbb{R} \rightarrow \mathbb{R}$, $f: \mathbb{T}^d \rightarrow \mathbb{R}$, 
$\omega \in \mathbb{R}^d$ (without loss of generality we can assume that $\omega\cdot k\not=0$ 
for any $k \in \mathbb{Z}^d\backslash\{0\}$). 
We assume that there is a $c_0\in \mathbb{R}$ 
such that $g(c_0) = \hat f_0$ (where $\hat f_0$ denotes the 
average of $f$) and $g'(c_0) \ne 0$. Special cases of this 
equation, for example when $g(x)=x^2$, are called the ``varactor 
problem'' in the literature. 
   
We show that if $f$, $g$ are analytic, and $\omega$ satisfies some 
very mild irrationality conditions, there are families of 
quasi--periodic solutions with frequency $\omega$. These families 
depend analytically on $e$, when $e$ ranges over a complex 
domain that includes cones or parabolic domains based at the 
origin. 
   
 
The irrationality conditions required in this paper are very weak. 
They allow that the small denominators $|\omega \cdot k|^{-1}$ grow 
exponentially with $k$. In the case that $f$ is 
a trigonometric polynomial, we do not need 
any condition on $|\omega \cdot k|$. 
This answers a delicate question raised in 
[8]. 
   
We also consider the periodic case, when $\omega$ is just a number 
($d = 1$). We obtain that there are solutions that depend 
analytically in a domain which is a disk removing countably many 
disjoint disks. This shows that in this case there is no Stokes 
phenomenon (different resummations on different sectors) for the 
asymptotic series. 
   
 
The approach we use is to reduce the problem to a fixed point 
theorem. This approach also yields results in the case that $g$ is 
a finitely differentiable function; it provides also very 
effective numerical algorithms and we discuss how they can be 
implemented.


Introduction
The goal of this paper is to study a differential equation describing a strongly dissipative dynamical system with a forcing function: where ε > 0 is a small parameter, g : R → R, f : T d ≡ [0, 1) d → R and ω ∈ R d are given. Note that we allow d = 1, which corresponds to the case of a periodic forcing. We will refer to all functions of the form of f as quasi-periodic, so that for us periodic is a particular case of quasi-periodic. We will mainly use the equivalent form of (1.1) εẍ +ẋ + εg(x) = εf (ωt) . (1.2) Using the terminology of electronic engineering, equation (1.2) with g(x) = x a and a ∈ [1.5, 2.5] is called the varactor equation and it has been considered several times in the literature ( [10], [11], [8], [9]). Notice that from the mathematical point of view, (1.2) is a singular perturbation problem because the small parameter affects the term with derivatives of the largest order.
In this paper we establish existence and local uniqueness of quasiperiodic solutions of frequency ω for all ε in a range around zero.
Under appropriate non-degeneracy assumptions on g and (rather weak) non-resonance conditions on ω, we prove the existence and uniqueness of families of quasi-periodic solutions with frequency ω, depending analytically on ε, when ε ranges over complex domains accumulating at the origin and including the physical values corresponding to ε real and positive. Since the argument covers already the physically relevant values of ε, we have not optimized the analyticity domain.
We stress that the non-degeneracy and Diophantine conditions required by our approach are very mild. For example, it suffices to assume that g (c 0 ) = 0 for c 0 such that g(c 0 ) equals the average of f . The non-resonance conditions we assume allow that |ω · k| −1 ≈ α exp(β|k|), k ∈ Z d \ {0}, for some constants α, β, and even that ω · k vanishes in some case (see Theorem 3 for precise statements). Note that the nonresonance condition above -which is not the most general condition we can deal with -is more general than Diophantine and even than the so-called Brjuno condition 1 ; this solves some questions raised in [8], [9]. In the case that d = 1 -the periodic case -we also obtain that the solution can be defined for complex values of ε in a disk with a countable number of smaller disjoint disks removed. In this domain we can also prove local uniqueness. Note that this domain is simply connected and it includes loops that enclose the origin in any arbitrarily small ball centered at ε = 0. This shows that the resummations of expansions of solutions do not present the Stokes phenomenon, namely one cannot obtain different resummations in different sectors.
We also note that the method is constructive and it leads to easy and fast numerical algorithms (see Section 8.1). The main theorems in this paper are in an a-posteriori format (that is, given a function that solves the equation approximately, we conclude that there is a true solution nearby), which allows us to be sure that the numerically computed solutions are close to true solutions (see Section 8.1). The method we present also yields results in the case that f and g are only finitely differentiable (see Section 8.2).
The method of proof in this paper consists first in constructing an approximate solution to finite order in ε and, then, in transforming (1.2) into a fixed point equation, which we can solve thanks to the fact that the strong dissipation leads to very few small divisors. Therefore, the leading assumption of the main theorem (see Theorem 3) is that there exists a perturbative solution to a finite order. This assumption on the existence of solutions to finite order is the only source of small divisor equations. The approximate solution is used only as a starting point of an iterative method that does not involve any small divisors.
In this paper we have decided to use only the elementary fixed point contraction theorem and not to use techniques based on quadratically convergent methods that also seem promising. Remark 1. The equation (1.2) is conformally symplectic; the paper [3] develops a Newton's method for maximal dimensional conformally symplectic systems. The theory of [3] would produce solutions in which the frequency of the forcing is supplemented by a natural frequency of , α n (ω) = inf{|2πω · k| : k ∈ Z such that 0 < |k| ≤ 2 n }.
the system. A Newton's method that applies to periodically forced systems and which does not depend on geometric considerations can be found in [14].
1.0.1. Relation with the literature. The method used in the present paper is very different from that used in [8], [10], [11], which is based on the construction of the full perturbative series and on showing that these series can be summed using Borel transforms. The construction of the perturbative expansions is affected by small divisors, because the expansions are computed from ε = 0, which is precisely the value that makes the equation singular. Furthermore, the Borel summability involves obtaining estimates in complex domains of ε, which are larger than the physical domain ε ∈ R. These two effects lead to assumptions in the small denominators, which are not present in the method used in this paper. We also point out that since the equations are non-linear, it is not clear that the Borel summation of a solution in the sense of formal power series is a solution in the standard sense of differential equations.
The method presented here goes in the opposite order than the methods based on resummation. It first produces an analytic function that solves the desired equation in the standard sense of differential equations. For us, the investigation of whether this function is obtainable from the asymptotic expansion is an afterthought that is only considered once the function is constructed and it has been shown to satisfy the differential equation.
We provide examples (see Example 12) in which the perturbative series does not exist, but we can produce solutions through the fixed point methods.
To provide some possible contact with the summability methods, we have shown that our solutions are defined in a domain which includes circles tangent to the imaginary axis at the origin. This is one of the hypotheses of a celebrated theorem on summability (see [17]). The main theorem of [17] shows that the analytic functions defined in domains containing circles tangent to the imaginary axis and whose asymptotic expansion satisfies factorial bounds are Borel summable (that is, they are determined by their asymptotic expansion and they can be reconstructed from their asymptotic expansion by the Borel transform).
In the same vein of providing contact with summability methods, in Theorem 6 and in Section 7 we have studied the periodic case and showed that the solutions of (1.2) can be defined and are locally unique for ε lying in a disk centered at the origin from which one removes a sequence of disjoint disks centered at the imaginary axis. In the periodic case, it was shown in [10] that the factorial bounds in the coefficients hold; hence, we can conclude that the solution is recovered from the Borel transform. Notice also that, since the domain we establish contains loops that enclose the origin within any small disk, we obtain that Stokes phenomena (the resummation in two sectors disagree) do not happen because the analytic function constructed here allows to perform an analytic continuation from one sector to the other.
1.0.2. Organization of the paper. This paper is organized as follows. In Section 2 we formulate the problem and we state the main result. Some preliminaries, such as the definition of function spaces and norms, are provided in Section 3. The construction of an approximate solution is given in Section 4. The proof of Theorem 3 is presented in Section 5 with the existence part given in Section 6. In Section 7 we provide the proof of Theorem 6 which considers the periodic case d = 1 and establishes a large domain. In Section 8.1 we describe how to turn the present proof into effective numerical algorithms. In Section 8 we present some results that follow from the same circle of ideas, such as results for finitely differentiable problems.

Formulation of the problem and statement of the main result
We proceed to set up the notation, formulate the problem and perform some manipulations transforming it into an equivalent form. In this section, we will not indicate the spaces where the problems are defined; this will be done later in Section 3, but it will be trivial to show that the elementary manipulations done here apply to these spaces (and indeed in almost any reasonable space).

2.1.
Quasi-periodic solutions of equation (1.2). For the functions f and g defining the problem (1.2), we will assume that: • H1 there exists c 0 ∈ R such that g(c 0 ) =f 0 , g (c 0 ) = 0, wherê f 0 denotes the average of f ; • H2 f , g are analytic functions. The first part of H1 holds if g has a range which covers all the reals. The second part is the nontrivial content of the assumption, that at the point c 0 the derivative of g is not zero. Note that hypothesis H1 is satisfied by Our aim is to find the so-called response solution, namely a solution x ε (t) of (1.2) of the form for a suitable function X ε : T d → R, indexed by the small parameter ε, and c 0 ∈ R.
Note that, of course, what we are seeking is a quasi-periodic solution, namely a function X ε : T d → R. We have just found notationally convenient to segregate the average c 0 because, as we will see later, the average needs to be dealt with in a different way. We are particularly interested in showing that the function ε → X ε is analytic, when ε ranges over a domain.
We emphasize that we are not assuming that ω is rationally independent. In particular, for us, periodic solutions are a special case of quasi-periodic solutions.
In terms of the hull function X ε = X ε (θ) parameterized by θ ∈ T d , the equation (1.2) reads as Note that the hull function X ε , being a function from the torus, is a periodic function in the variable θ. On the other hand, the response solution x ε (t) = c 0 + X ε (ωt) is a quasi-periodic solution.
We will introduce the notation Note that M ε is an operator acting on periodic functions of θ ∈ T d . For convenience, we also introduce the notation . With the notations above, (2.1) can be rewritten as: As we will see later, under appropriate hypotheses on ω, M ε will be invertible, so that (2.3) can be transformed into a fixed point problem.
We will show that M −1 ε is defined in a suitable complex domain accumulating at ε = 0, but its norm will grow when ε approaches 0. Nevertheless, if we start with an approximate enough solution, we will be able to transform the problem further (as in the implicit function theorem), so that we obtain a solution of (2.4) and we will show that it is analytic in ε in a large complex domain that accumulates to ε = 0. This method also works for finite differentiable functions, even if in this case the conclusions are only that the solution is defined for a range of real ε accumulating at ε = 0 and that it depends on ε in a differentiable way.

Approximate formal solutions.
We define an approximate formal solution as follows.
Definition 2. We say that X ≤N ε (θ) = N j=1 ε j X j (θ), where all the X j are analytic functions, is an approximate solution of (2.1) to order N , whenever it is The existence of solutions to arbitrarily high orders is discussed extensively in [8]. We will revisit some of these results in Section 4. As mentioned in the introduction, we will not need to study convergence or re-summation of these series expansions. For our method, just the existence of a finite number of terms is enough. Hence, the non-resonance assumptions we will need are much weaker than the assumptions adopted in [8].

2.3.
Formulation of the main result. The main result of this paper is the following Theorem that shows that the existence of formal solutions to a low order implies the existence of a solution analytic in a domain.
The solution X ε is locally unique and it is asymptotic to X ≤N ε in the sense that for ε ∈ Ω σ,B we have: , for some positive constant C, where m > d and || · || ρ,m is a norm on spaces of analytic functions that will be introduced in Section 3.
In Section 3 we will be more precise on the domains of analyticity, namely on the definition of the norm · ρ,m (see Definition 9) and on the relations between B, σ, ρ, but this will require to introduce more notations and in particular to formulate the results in terms of some well defined spaces (see Section 3).
Note that the statement of Theorem 3 does not need to assume any Diophantine property on ω. Therefore, in particular, Theorem 3 works for any frequency, including ω rational. Of course, to apply the result to concrete systems, we need to show that there exist asymptotic solutions to a finite order bigger than 1 and this entails some (very mild) small divisor conditions. We can also take advantage of the fact that the perturbation has a special structure (e.g., it is a trigonometric polynomial).
Remark 4. We define conical domains Υ δ,σ as Note that the domain of analyticity in ε, θ includes sets of the form Υ δ,σ × T d ρ for 0 < δ, σ, ρ 1. In the case of the conical domains, we can give a proof that provides more information on the rates of convergence of the iterative process used in the proof of the theorem and as a basis of the algorithms. Furthermore, to obtain the existence and uniqueness, we need only to assume that there is an approximate solution to order N ≥ 1 (see Section 6.2), instead than requiring that N ≥ 2 as in the parabolic domain. Hence, we have included Section 6.1 with a proof of these results. We hope that this will illustrate also the fact that the method presented here is very flexible and admits many variants.
Remark 5. We note that the domain Ω σ,B includes disks tangent to the imaginary axis. This is one of the hypotheses of the main theorem in [17].
The result in [17] shows that if a function X ε defined for ε ∈ Ω σ,B admits a formal expansion, for some positive constants C, γ, then the function X ε can be obtained by Borel summation from the expansion.
Notice that the domain Ω σ,B does not include the imaginary axis, but that it is tangent to it.
For sectors that include the imaginary axis, there is a much better known result provided by Watson theorem ( [19], [12]). The width π of the sectors is crucial because of the well known examples exp(−ε 2 ), which have identical expansions in sectors of width close to π.

2.4.
Results in the periodic case, d = 1. In the periodic case d = 1, we can prove the existence and analyticity of the response solution in a domain which is obtained removing a sequence of disjoint disks from a small disk centered at the origin. where where R is chosen sufficiently small Then, we can find a solution X ε of (2.1) defined for ε ∈ D R . For a fixed ε, then X ≤N ε (·) ∈ A ρ (ρ > 0 independent of ε) and it depends analytically on ε. The solution X ε is asymptotic to X ≤N ε and it is locally unique in a neighborhood of the asymptotic solution.
Note that the domain D R in (2.7) is connected and it includes loops enclosing the origin, so that when we have a function defined on it, it cannot present the Stokes phenomena (different resummations in different sectors).
Remark 7. Note that, there are at most two k for which (ωk) 2 = g (c 0 ). For most values of ω, there will be no such k. In such a cases, ε k is not defined or is infinity. Theorem 6 remains valid if in the resonant case, we set ε k = ∞ and therefore B r k (ε k ) ∩ B R (0) = ∅.
To avoid complicating the statement of the Theorem 6 we adopt the convention ε k = ∞, when k is resonant, which makes the argument go through.

2.4.1.
Some indications on the proofs. The construction of the asymptotic expansions in the assumption H3 of Theorem 3 will be undertaken in Section 4. This will require several hypotheses of Diophantine properties of ω and of non-degeneracy of the problem.
We will show also that the formal power series expansion is uniquely determined by the fact that it satisfies the equation (1.2) in the sense of power series expansions. Hence, the final result of this paper will be that the hypotheses discussed in Section 4 which imply H3, also imply the existence of an analytic function in a domain. The reason why we have formulated the main result in this way is because there are several sets of hypotheses that lead to H3. Another reason to split the results in two steps is that the proof of Theorem 3 is based on very different techniques, than the verification that the other hypotheses imply H3.
As we will see later, the Diophantine assumptions we make are much weaker than the so called Brjuno conditions. Thus, we settle a question raised in [8] on whether it is possible to obtain solutions when the frequency has Diophantine properties weaker than Brjuno.
Since, for some technical purposes, it seems desirable to establish larger domains of analyticity in ε, we start by giving the proof in conical sectors (see Section 6.1); then we extend it to parabolic domains requiring slightly stronger hypotheses, but obtaining a larger analyticity domain (see Section 6.2). The results presented are far from optimal, but we hope to convey that the method presented here is flexible. One motivation to obtain larger analyticity domains is to study the summability properties of the expansions.
We will also state some results for finitely differentiable functions. In such a case, we will only conclude that the response functions can be constructed for ε ∈ (0, ε 0 ] for some 0 < ε 0 1 and that they are finitely differentiable (see Section 8.2).
The methods also work in cases that the x is higher dimensional. We discuss this very briefly in Section 8.3.
We note that Theorem 3 will be proved by a contraction mapping argument, so that if we give an approximate solution (e.g. a formal power series or a numerical computation) that solves the invariance equation approximately, then there is a true solution close to the approximate one. This is often called an "a-posteriori" formulation of KAM theory, which can be used to validate approximately computed solutions.

Some spaces of analytic functions
To implement the program outlined above, we need to make precise the spaces in which we are considering the problem. Note that the discussion in Section 5.1 just requires that the space is characterized by Fourier coefficients -so that indeed M ε is invertible -and that the composition operator X → G • X is smooth, when X is given the appropriate topology. With a bit of foresight for further applications, we would also like that the spaces can be easily adapted to higher dimensions or to the case that G is only finitely differentiable.
Given a function f : T d → R, we denote its Fourier expansion as If the function f is analytic, it satisfies the Cauchy bounds for some positive constant C, where ρ > 0 is a lower bound on the radius of analyticity around any point and |k| = |k 1 | + · · · + |k d |.
Definition 9. Given ρ > 0, m ∈ N, we define for a function X : We denote by A ρ,m the Banach space of functions X, such that ||X|| ρ,m is finite endowed with this norm.
Remark 10. Note that for a multi-index, we use |k| to denote the 1 norm of k (the sum of the coefficients). The choice of a norm of k is important for the exponential term in the norm and different equivalent norms for k lead to different spaces of analytic functions. In contrast, for ρ = 0, any choice of the norm of k leads to an equivalent norm in the space of functions.
The spaces A ρ,m are very used in analysis; an important property is that the norm can be expressed as an integral. In fact, by Parseval identity, we have that ||X|| 2 ρ,m is comparable to integrals up to a positive constant, say C: where the integral is understood with respect to the 2d-dimensional Lebesgue integral in T d ρ . Note that for ρ > 0 the functions in A ρ,m are analytic in the interior of T d ρ and extend to Sobolev functions on the boundary of T d ρ . In particular, when m > d they extend to continuous functions on the boundary of T d ρ and the A ρ,m -norm is stronger than the C 0 (T d ρ )-norm. For our purposes, it will be enough to remark that when m > d it is very easy to see that the space A ρ,m can be identified with the closed space of H m (T d ρ ) -the standard Sobolev space -consisting of analytic functions (notice that the Sobolev embedding theorem shows that the limit in H m is also a limit in C 0 and it is very elementary that the uniform limit of analytic functions is analytic. Variants of this are also true in lower regularities, but we will only need the elementary results.) When ρ = 0 we recover the standard Sobolev spaces, since A 0,m = H m (T d ). Therefore, by Sobolev embedding theorem, we have that A 0,m embeds into the space of continuous functions for m > d (note that we are using 2d-dimensional integrals, [18]). However, for ρ = 0 the integrals we are considering in (3.3) are d-dimensional Lebesgue integrals, rather than the 2d-dimensional Lebesgue integrals that we were considering when ρ > 0. This is a rather substantial change in the notation, which causes that the m index needed for the Banach algebra changes by a factor of 2.

Construction of the approximate solution
In this section we show how to construct approximate solutions and to obtain some preliminary bounds. This is, of course, the realm of the standard formal perturbative theory. We note that there is a compromise between the complication of the non-linearity and the small divisor conditions that are needed. We present several results in this direction. We note that the conditions needed are very mild. In particular, they are milder than the Brjuno condition, for example Lemma 11 allows for rational frequencies. Of course, Theorem 3 shows that once we have these approximate solutions, we have a true solution.

Trigonometric polynomials.
Lemma 11. Assume that • H3.a f is a trigonometric polynomial of order K and ω is nonresonant up to the order N K, i.e. ω · k = 0 for all k ∈ Z d \{0} with |k| ≤ N K, where N is the order of the approximate solution. Then, there exists an approximate solution to order N ≥ 1 of (2.1) (compare with H3).
Proof. Let us expand X ε (θ) in Fourier-Taylor series as (4.1) we show that we can define an approximate solution X ≤N ε (θ) of (2.1) as a finite truncation of (4.1) up to the order N . Taking into account To first order in ε we obtain: for k = 0 we have g(c 0 ) =f 0 providing whileX 10 is unknown and it will be determined at next order. For k = 0, |k| ≤ N K, we obtain so that X 1 is well defined, due to the non-resonance condition up to the order N K. One can proceed recursively up to the order N by defininĝ X N k as the solution of the equation (4.5) can be solved provided that ω · k = 0 for all the k for which (S N ) k = 0. Finally, it is well known (see e.g. [6]) and easy to see (using recursion relations for the expansions of trigonometric functions of the composition) that the assumption that f is a trigonometric polynomial of order K implies that S N is a trigonometric polynomial of order N K. Hence, under the assumption H3.a, we can perform the step N − 1 times and obtain a solution to order N .
Note that, the way that we obtainX N,k in (4.5) is to divide by ω · k. Hence, the size of these small divisors is crucial. In the periodic case (d = 1), we have |ωk| ≥ |ω| > 0 when k = 0; hence, in the periodic case, the divisors are not small. When d ≥ 2, one sees that ω · k might become arbitrarily small for large enough k. Indeed, the results are stronger in the periodic case (see Theorem 6 and Section 7).

4.2.
Some examples of analytic solutions without asymptotic expansions. The construction of Lemma 11 leads to the following example that shows that there are situations where we can apply Theorem 3 to construct solutions in the physical domain ε > 0, but nevertheless the formal series expansions cannot be defined beyond a finite order.
Example 12. Given g satisfying g (c 0 ) = 0 (c 0 as before) and given integers N 0 > 1, K > 0, consider a frequency ω such that ω · k = 0 for all |k| < N 0 K, but such that ω · k = 0 for some |k| = N 0 K. Then, we can find trigonometric polynomials f of degree K, such that it is possible to find approximate solutions up to order N 0 − 1, but there is no approximate solution up to order N 0 .
Clearly, for the trigonometric polynomials in Example 12 we cannot find any approximate solution to order N 0 -a fortiori, we cannot find an asymptotic expansion to all orders. Nevertheless, Theorem 3 shows that we can obtain an analytic solution defined in a large set (e.g. a parabolic domain).
In fact, we recall that, due to the non-resonance condition of Example 12, Lemma 11 allows to find solutions up to the order N 0 − 1, no matter what the trigonometric polynomial is.
Hence, to verify the claim in Example 12, we just need to verify that it is possible to construct polynomials in such a way that it is impossible to find the coefficient X N 0 . This follows easily because of the fact that when we look at the equations of order N 0 we have to solve the equations (4.5).
The key observation is that, if there is a k with |k| = N K 0 , the coefficient (S N 0 ) k will be an algebraic expression of the coefficients of f which are not identically zero. This can be easily seen because if we take for example k = (N 0 K, 0, . . . , 0) in the expression of (S N 0 ) k , there is one term of the form A(f (K,0,...,0) ) N 0 where A is a non-zero factor.
Remark 13. Note that the above argument shows that if we fix a resonance of order N , the only requirement on the trigonometric polynomial to exhibit the phenomenon described in Example 12 is that the coefficient in the expansion at order N corresponding to the resonance does not vanish. Since the coefficient of the perturbative expansion is a polynomial in the coefficients of f , and we have shown that the polynomial is not trivial, the set of zeros will be a codimension 1 variety in the space of coefficients.
Therefore, given a resonant frequency as in the assumptions of Example 12, if we consider the space of coefficients of the trigonometric polynomial, there will not be a perturbative expansion to all orders, except if the coefficients lie on a codimension 1 variety.

Approximate solutions for general analytic forcing functions.
Eliminating the assumption that f is a trigonometric function, we obtain the following result. Note that the Diophantine conditions that we assume are much weaker than the Brjuno condition.
With some stronger assumptions on the Diophantine properties of the mapping, we obtain the existence of solutions to all orders. The assumptions are still weaker than the Diophantine conditions. Proposition 16. Let f ∈ A ρ,m with zero average. Assume that, for some 0 < η < ρ, α > 0, we have |ω · k| −1 ≤ α exp(2πη|k|). Then, there is a unique solution u of with u having zero average. Furthermore, we have that for some constant C > 0.
Proof. Due to the definition of the norm in (3.2), we obtain: for some constant C > 0.
The following result is also proved by the same method, but it is very standard. We recall that we say that ω is Diophantine when for some ν ≥ 1, τ ≥ 1.
Proposition 17. Let f ∈ A ρ,m (where ρ ≥ 0) with zero average. Assume that ω satisfies (4.7). Then, there is a unique solution u of (4.6) with u having zero average. Furthermore, we have that for some constant C > 0.
Note that Proposition 17 remains valid in the case ρ = 0.

Proof of Lemma 14.
To prove Lemma 14 we argue that, as in Lemma 11, X 1 is given by (4.2) as In a similar way, we can proceed as in Lemma 11 to define the solution up to the order N , which is given as the solution of the equation for a suitable function T N (θ) = k∈Z dTN k e 2πik·θ defined in terms of c 0 , X 1 , ..., X N −1 and such that the right hand side of (4.8) has zero average, thanks to the choice ofX N −1,0 . Then, we have: The hypothesis H3b.2 guarantees that we can apply Proposition 16 to get estimates of the solutions X j , j = 1, ..., N , on domains of size ρ − jη/N , so that at each step we can take a domain loss equal to any number bigger than η/N . The result of Lemma 15 is very similar to those stated in [8], [10], [11], except that the Diophantine conditions considered here, namely H3b.2, are weaker than the Brjuno conditions considered in [11] and a fortiori than the Diophantine conditions considered in [10]. The question of whether it is possible to weaken the Diophantine conditions was raised in [8, Section 6].
In Section 8.2 we will present a result very similar to Theorem 3as well as analogues of Lemmas 11 and 14 -for finitely differentiable problems.

Proof of Theorem 3: estimates on
We will prove Theorem 3 by formulating the problem as a fixed point theorem in Banach spaces. Of course, we will have to specify the Banach spaces, which will make Theorem 3 more quantitative than hitherto stated.

Preliminaries and formal considerations.
The key observation is that for ε = 0, the operator M ε defined in (2.2) is boundedly invertible from some spaces of analytic functions to themselves. Note that in Fourier coefficients the operator M ε is diagonal: where m ε (t) = −εt 2 + εg (c 0 ) + ıt . Therefore, we can obtain estimates on the operator M −1 ε (in many Banach spaces of functions, whose norm is determined by the Fourier coefficients) by estimating Of course, the estimates for Γ will be singular as ε tends to zero, but it will be apparent that if ε converges to zero away from the imaginary axis, one can obtain estimates (see Section 5.3). Moreover, we note that if we consider a Banach space which is also a Banach algebra, the operator X → G • X will be Lipschitz and, given that G is quadratic, we have that the Lipschitz constant will be small in a ball of small radius. Hence, if we consider X in a space of functions which is at the same time determined by the Fourier coefficients and a Banach algebra under multiplication, we will show that the fixed point problem has a solution for ε in a suitable domain, which we will make explicit. Note also that in the periodic case we can obtain much better estimates for M ε . In the periodic case we have that ω · k is bounded away from zero, so that it suffices to study As we will see, in the periodic case taking the sup over 2πωk will be much smaller than taking the sup over the whole real line. Of course, when ω is higher dimensional, since ω · k is dense on the real line, it is the same to take the sup over ω · k or over the whole real line. In Sections 5.2, 5.3 we will present the estimates that allow to apply the contraction argument and to obtain response solutions of (1.2).

Estimates on multiplication and composition operators.
We have the following well known result (see, e.g., [18]).  is analytic as a mapping from A ρ,m to itself. If the mapping X = X ε depends analytically on the parameter ε, then Λ G depends analytically on the parameter ε. Furthermore, we have We note that, because G is analytic, we have that G(t) = n G n t n and that we have for some η > 0 the estimate |G n | ≤ Cη n . We see that, by the Banach algebra property of the norm, one has ||G n X n || ≤ Cη n ||X|| n . Hence, when ||X|| is sufficiently small, the series n G n X n converges uniformly.

5.3.
Estimates on the operator M ε . We start with some elementary considerations on the computation of the supremum involved in the definition of Γ(ε) in (5.1). Note that we just need to estimate the minimum of the modulus of a quadratic complex polynomial in t. This can be readily transformed into computing the minimum of a quartic real polynomial, which is a rather straightforward albeit tedious problem. We want to pay, of course, special attention to the case when ε is small.
We note that, once we obtain estimates on the infimum of the multiplier Γ, we obtain estimates for M −1 ε in all the spaces whose norms are obtained from sizes of the Fourier coefficients. This includes, of course, the spaces A ρ,m for all values of ρ ≥ 0, m ∈ N and also the Sobolev spaces with ρ = 0, which we will use for the results on finite regularity in Section 8.2.
Whenever | Im (ε)| ≤ α|ε|, α 1, we can use first order perturbation theory for the value of the minimum and conclude that Γ(ε) ≤ C|εA| −1 on a conical sector that includes the real axis.

5.3.2.
Estimates on the imaginary axis. Note that, when ε = is, For fixed s, the equation S(is, t) = 0 is a quadratic equation in t so that its discriminant is 1 − 4(−s)sg (c 0 ). When g (c 0 ) ≥ 0, for any s we can find real roots for t. When g (c 0 ) < 0, real roots exist for small enough s.
Therefore, when ε = is, the operator M ε is unbounded and the elementary fixed point argument used in this paper does not apply (of course, more sophisticated methods such as KAM theory ( [3]) may work under stronger resonance assumptions).

Estimates in a parabolic domain.
To study the analytic properties of the function X ε (in particular, whether it can be reconstructed by Borel summation of its asymptotic expansion), it will be interesting to study its properties in domains which are tangent to the imaginary axis. We fix with r ∈ R and B > 0 large enough, say B > B 0 for some B 0 ∈ R + . We aim to prove the following result.
Proposition 20. When ε is on the the parabola (5.2), we have: for a suitable positive constant C.
To prove (5.3), we start by writing S as The first term vanishes at t = ± √ A. We can define three regions in t, namely and the complement of I + ∪ I − . If t belongs to the complement, we have the estimate for as suitable positive constant C s . Then we have In the intervals I − and I + , we can bound the first term in S(t) from below by zero and we have Taking the infimum we have that Γ(ε) satisfies the bound (5.4) for any t.

Proof of Theorem 3: the existence part
The proof of the existence of the solution in Theorem 3 relies on a fixed point theorem on Banach spaces. We recall that we proved in Lemma 11 and Lemma 14 the existence of an approximate solution X ≤N ε as in (2.5). We can write (2.4) as We define the operator T ε , acting on functions analytic in ε and taking values in A ρ,m , as The equation (6.1), equivalent to (1.2), is a fixed point equation for T ε . To study this equation, we need to find the domain of T ε . Then, we will show that the operator T ε maps the domain into itself and that it has a Lipschitz constant smaller than 1.
We will present two versions of the arguments. One in cones including the physical sector and another argument including parabolic domains. Of course, the parabolic domain is more general, but we will present both, because they involve different considerations and in the cone case we will obtain faster convergence. The cone case will also generalize better to finite differentiable functions.
6.1. Proof of Theorem 3 in conical sectors. In this section we provide a proof of a slight modification of Theorem 3, since we will consider domains which are just sectors. We will need to assume that the approximate solution is a solution at least to order 1 and we will obtain explicit estimates on the convergence.
For the subsequent estimates, it will be important for us that the domain Υ δ,σ (see its definition in (2.6)) does not include the origin. Hence, we can bound from above the singular factors by σ −1 , since σ is the minimum distance to the origin in the domain. Also, we see that the maximum distance to the origin can be bounded from above by a constant times σ. As we will see later, this will allow us to show that the bad factors are dominated by the good factors. It will be quite important that in the domain we consider, the distance is bounded away from zero and that the maximum and the minimum distances are comparable.
Remark 21. Note that, once we consider the operator T ε acting on spaces of analytic functions in ε, we will obtain automatically that the solution depends analytically on parameters. On the other hand, we note that the way that we prove the contraction in the sup-norm is to prove estimates for a fixed ε and take the supremum. Therefore, the argument we present also contains as a particular case the convergence for each fixed value of ε. In the numerical considerations (or when we consider only the real values), it is advantageous to consider just fixed values of ε and to obtain only pointwise convergence. This will also play a role in establishing uniqueness and analytic continuation (see Corollary 22).
Due to the Banach algebra property of A ρ,m and to the fact that G (0) = 0, we have that for X ε , Y ε in a ball of radius α (α < α 0 ) centered at the origin of A ρ,m,δ,σ , the operator T ε is Lipschitz with constant K 1 α, namely for some constant K 1 > 0. The reason is that, as we have shown, we can bound ||M −1 ε || by a constant and that the Lipschitz constant of the composition with G on the right is bounded by a constant times α.
Assume that we are given an approximate solution X ≤N ε of (6.1) such that X ≤N ε ρ,m ≤ K 2 σ for some K 2 > 0. Consider a ball of radius β around X ≤N ε ; then, we have that T ε is contractive, if K 1 (K 2 σ + β) < 1 , provided K 2 σ + β ≤ α . Moreover, we know that the approximate solution satisfies . This is because we know that X ε solves the equation (2.1) up to O(ε 2 ) and we obtain the fixed point equation just by multiplying by M −1 ε , whose norm is bounded by |ε| −1 .
The ball of radius β in A ρ,m centered around X ε gets mapped into itself if K 3 σ + K 1 (K 2 σ + β)β ≤ β . To conclude the argument it suffices to show that given K 1 , K 2 , K 3 , the properties of the operator and the approximate solution, for all 0 < σ ≤ 1 (the size of the annulus in ε-space), we can choose α as the radius of a ball centered at zero in function space in which we are doing the a-priori estimates, β as the radius of a ball centered around the approximate solution, in such a way that we can fulfill the following conditions: 3) The second condition in (6.3) implies that the ball B β (X ε ) = {Y ε ∈ A ρ,m : X ε − Y ε ρ,m ≤ β} and T ε (B β (X ε )) are contained in the ball in which we have a-priori estimates for the Lipschitz constant in the ball B α (0); the third condition implies that T ε (B β (X ε )) ⊂ B β (X ε ) and the first one that T ε is a contraction on B α (0), a fortiori in B β (X ε ). Now we verify that indeed, for all σ > 0, we can find the desired α, β. To simplify the computations, we can take β = 2K 3 σ, so that (6.3) becomes which can be fulfilled by a proper choice of the constants. This implies that the operator T ε admits a fixed point, which is the solution of the equation (1.2) in the domain Υ δ,σ . Notice that one of the conclusions of the theorem is that we get local uniqueness for all the solutions. In particular, we obtain the following result.
Corollary 22. If we fix ε and denote by X ≤N ε an approximate solution, then for any θ ∈ T d Furthermore, we observe that since the contraction result is true for all σ sufficiently small, we can obtain the existence and local uniqueness for σ and 0.9σ. The pointwise uniqueness in Corollary 22 shows that the two functions produced in Υ δ,σ and Υ δ,0.9σ have to agree in the intersection of the two domains. Hence, they are analytic continuation of each other. Therefore, we can obtain a solution of (2.1) defined for ε ∈ ∪ 0<σ<R Υ δ,σ for some R > 0. This solution is locally unique because it can be obtained as the pointwise limit of the repeated iteration of T ε starting with the approximate solution.
Notice also that, if we started with an approximate solution satisfying the equation to a higher order, the true solution would differ by the same order.
6.2. Proof of Theorem 3 in parabolic domains. The proof is very similar to the proof of Theorem 3 in Section 6.1, but the argument is quantitatively different.
We will consider T ε defined in the space A ρ,m,B,σ consisting of the analytic functions of ε taking values on A ρ,m and with ε ranging on the domain Ω σ,B .
Following the argument before, but taking into account that, for In the ball B β (X ≤N ε ) the map T ε will be a contraction, if the upper bound on the Lipschitz constant is smaller than 1, namely Provided, of course, that As before, we want to show that, given K 1 , K 2 , α 0 , for all B > B 0 , σ < σ 0 , it is possible to find β > 0 such that (6.4), (6.5), (6.6) are satisfied. We note that if we take β = Aσ for some A > 0, the conditions (6.4), (6.5),(6.6) become For example, by choosing A ≥ 10K 2 , the first and third conditions in (6.7) are satisfied by taking B sufficiently large, while the second condition is fulfilled by taking σ small enough. The rest of the proof does not need any change from the proof in the sectorial domains.
7. Proof of Theorem 6 7.1. Some auxiliary estimates. To prove that (6.1) is a contraction, we will find it convenient to estimate A motivation for the study of this function is that ||M −1 ε || ≤ Γ per (ε) and that the Lipschitz constant of T ε , defined in (6.2) in a neighborhood of the approximate solution, is approximately |ε| 2 ||M −1 ε || (see the end of the argument in Section 7.1.1).
We consider k fixed and we study each of the functions |ε| 2 |m ε (2πωk)| −1 ; the desired estimates (7.1) follow by taking the supremum over k. The goal is to show that if we fix a certain level η ∈ R + , we have the bounds |ε| 2 |m ε (2πωk)| −1 ≤ η except in a small ball centered at a point in the imaginary axis or outside of a large ball.
To simplify the typography, we write We have: The function ϕ k (t) has a singularity at t = |A k |/|B k |; having fixed a level set η, the solutions of ϕ k (t) = η are determined by solving the pair of quadratic equations: Each of the above equations admit two roots t ± , one close to the singularity, say t − , and the other one far away, say t + . In particular, the first equation in (7.2) has the roots where t − can be expanded as where the last term denotes the displacement from the singularity and it is of order k −4 . An identical analysis, changing the sign of η yields the result for the second equation.
In conclusion we have that, for some constant C > 0 and for |ε| ≤ 1, we have ϕ k (|ε|) ≤ η except in the annulus: Now, we proceed to study the function |ε| 2 |m ε (2πωk)| −1 more carefully in the annulus (7.3). We denote ε k = A k /B k and we note that, for large k, ε k ≈ ik −1 ω −1 . We just observe that Then, if |ε − ε k | ≥ |k| −3 and ε is contained in the annulus (7.3), we have: Therefore, for any given η > 0, we can ensure that for large |k| the function |ε| 2 |m ε (2πωk)| −1 ≤ η for all |ε| ≤ 1 except for those in a small disk. For any η > 0, we choose R > 0 sufficiently small in order to obtain for all ε such that |ε| ≤ R, |ε − A k /B k | ≥ |k| −3 . We note that the domain in which we obtain the bound (7.4) is a disk from which we have removed a countable collection of disjoint sets. 7.1.1. The contraction argument in the periodic case. We note that in the annulus U σ ≡ {σ ≤ |ε| ≤ 2σ}, since the values of |ε| are comparable to σ, we have that σ 2 ||M −1 ε || is comparable to the function Ψ we estimated in Section 7.1.
Given any η, we have shown that for σ small enough, on the set U σ ≡ U σ \ B |k| −3 (A k /B k ) we have that Ψ(ε) ≤ η. As before, we will consider T ε defined on spaces of analytic functions of ε, θ, with ε ranging overŨ σ and we endow them with the norm of the sup in ε of the A ρ,mnorm of the function of θ obtained fixing ε.
As before, we note that starting with an approximate solution X ≤N ε restricted toŨ σ , we have that ||X ≤N ε || ρ,m ≤ K 2 σ for some K 2 > 0. Therefore, using that inŨ σ , |ε| can be bounded from above and below by σ times a constant, we obtain that the Lipschitz constant of T ε in a ball of radius β = Aσ around the approximate solution is bounded from above by σ|(K 2 σ + β)K 1 | ||M −1 ε || for some K 1 > 0, provided that K 2 σ + Aσ ≤ α 0 . (7.5) We will therefore obtain that T ε is a contraction in the ball when, because the left hand side of (7.6) is an obvious upper bound on the bound of the Lipschitz constant of T ε . Using that X ≤N ε satisfies the equation at least to order 1, we obtain that ||T ε (X ≤N ε ) − X ≤N ε || ρ,m ≤ K 3 σ. Therefore, we can ensure that the ball (7.7) Now, we can choose A, η, σ in such a way as to have the three conditions (7.5), (7.6), (7.7) are satisfied. This shows that, for σ small enough we obtain a solution which is again unique. Therefore, arguing as before, we obtain an analytic function defined for ε ∈ ∪ 0<σ<RŨσ , which coincides with D R in (2.7).

Some further results
8.1. The numerical implementation. The treatment of the problem we have presented is very well suited for the numerical calculation of the invariant objects. One possibility is to break the calculation in two stages following the rigorous treatment presented below.
A) In a first stage, one computes a few terms of the perturbation expansion as in the proofs of Lemma 11 and Lemma 14.
We note that the operations required are very similar to the standard calculation of the Lindstedt series in celestial mechanics, which have been performed in computers since the early 60's ( [2], [15]). One just needs to develop a toolkit to manipulate Fourier-Taylor series (arithmetic operations plus composition with some elementary functions). Some modern surveys with accessible code can be found in [13], [7]. B) In a second stage, for a fixed value of ε, it is convenient to revert to the fixed point method presented in (6.1). As shown in Section 6, iterating the right hand side of (6.1) will yield a contraction for small enough ε. The results in Section 6 also show that if we take as initial points the evaluation of the polynomial approximations obtained in part A), the iterations converge.
The implementation of (6.1) entails just composing with G and applying the operator M −1 ε . The composition is fast if we discretize X ε in a grid of points. On the other hand, the application of the operator M −1 ε is diagonal in Fourier space. Of course, one can pass from the representation in discrete points to the Fourier representation using FFT, which is a fast algorithm which has optimized implementations in practically every computer. Hence, it is very easy to implement a contraction algorithm for small ε.
Note also that the implementation of the iterations for different values of ε starting from the same polynomial guess is very easy to parallelize.
Note that the above algorithm is very well suited for the case when ε is very small. This is the regime when the equations are very stiff and conventional numerical integrators (e.g. Runge-Kutta, Adam-Bashforth, etc.) have difficulties in producing results.
One problem we have not considered is how to obtain numerical results in the optimal domain; on the other hand, we note that for the values of ε such that our operator is not a contraction, standard numerical integrators will work very well.

8.2.
The finite differentiable case. The same strategy presented here applies also to the case that g and f are finitely differentiable (but with sufficiently high derivatives). In this case, of course, we consider only that ε ∈ R.
As in the analytic case, we proceed to study first the existence of finite order approximations, then we reduce to a fixed point problem.
One important technical lemma is the following result.
Similarly, if we have g ∈ C m+ +1 , we obtain that there are expansions in the composition to order with a remainder of order + 1.
Proposition 23 is a well known result. A concise proof can be found in [4] (see also [18]).
There is a very extensive literature on this problem leading to many variations (sometimes C g is called a Nemitskii operator [1], [16]).
In analogy with the analytic case, we have the following result.
Theorem 24. With the same notations as in Theorem 3, assume that m > d/2. Assume furthermore that there exists an approximate solution to order N ≥ 2 and that g ∈ C N τ +m+2 . Then, for all |ε| sufficiently small there exists a solution X ε ∈ A ρ,m of (6.1) (and hence of (2.1)). Furthermore, the mapping sending ε to X ε is C N τ +m+2 , when X ε is given the A 0,m topology.
The proof of Theorem 24 goes along the lines of that of Theorem 3 and we sketch here the main steps. 8.2.1. Existence of finite order approximations. The existence of approximate solutions can be obtained by following the formal power series solutions. There are, as before, trade-offs between the assumptions on small divisors and the assumptions needed on the regularity of f .
As in Lemma 11, we point out that if f is a trigonometric polynomial of order K and g ∈ C N +1 , we can compute the expansion to order N in A 0,m provided that ω · k = 0 for 0 < |k| ≤ N K.
Also, in the case that ω satisfies (4.7) and that g ∈ C N τ +m+2 , we can find solutions to order N . To verify the latter statement, we just note that we can perform the expansions to order N and that, to match the solutions, we can use use Proposition 17.
8.2.2. The fixed point argument. As in the analytic case, we argue that the operator at the right hand side of (6.1) is a contraction in A 0,m for sufficiently small ε and that it maps a ball around the approximate solution into itself. The fixed point theorem provides a solution of (6.1) as claimed in Theorem 24. 8.3. Higher dimensional results. The method of proof adapts very easily to the case that x is multidimensional x : R → R and, of course, f is also a function taking values in R and g : R → R .
It suffices to assume that there exists a c 0 such that g(c 0 ) = f and that Dg(c 0 ) is invertible and, for convenience, diagonalizable. Notice that after a linear change of variables, we can, of course, assume that Dg(c 0 ) is diagonal. Therefore, the linear part of the problem becomes just uncoupled copies of the problem we have studied, so that the same estimates we have obtained for M −1 ε apply also to the multi-dimensional case. The estimates on the nonlinear part G remain the same. Hence, the fixed point argument remains identical.
The existence of formal approximate solutions is again copies of the argument before, so that the results remain true only with a change of notation.