On the stability of periodic orbits in delay equations with large delay

We prove a necessary and sufficient criterion for the exponential stability of periodic solutions of delay differential equations with large delay. We show that for sufficiently large delay the Floquet spectrum near criticality is characterized by a set of curves, which we call asymptotic continuous spectrum, that is independent on the delay.


1.
Introduction. Delay-differential equations (DDEs) are similar to ordinary differential equations (ODEs) except that the right-hand side may depend on the past. For example, they could be of the forṁ where x(t) is a vector in R n and the delay τ > 0 decides how far one looks into the past. When studying DDEs as dynamical systems one notices that equilibria do not depend on the delay τ . More precisely, their location and number is independent of τ . However, their stability changes significantly when one varies τ , an effect that is well known and of practical importance in engineering and control [6,12]. Exponential stability is given by the spectrum of the linearization of the DDE in its equilibrium. This spectrum, in turn, can be expressed as roots of an analytic function (a polynomial of exponentials λ → exp(−λτ )). Lichtner et al [5] classified rigorously which types of limits this spectrum can have as τ approaches infinity. Roughly speaking, for sufficiently large τ all except maximally n eigenvalues form bands near the imaginary axis. After rescaling their real part by 1/τ one finds that these bands converge to curves, called asymptotic continuous spectrum. They are given as root curves of parametric polynomials, and are, thus, much easier to compute than the eigenvalues of the singularly perturbed large-delay problem. Of practical relevance are then stability criteria based entirely on the asymptotic spectra that guarantee the stability of an equilibrium for sufficiently large delays τ . This paper gives a similar result for periodic orbits of (1). In contrast to equilibria, periodic orbits change as the delay τ varies, such that the statement about "independence" of the delay (or, rather, re-appearance) has to be formulated more carefully. Let us look at a two-dimensional example to illustrate the observation made by Yanchuk & Perlikowski [16]: where we fix α ≈ −0.1 and vary τ . System (2)-(3) consists of the normal form for the supercritical Hopf bifurcation with an additional delayed term y(t − τ ) in the second equation, which breaks the rotational symmetry of the instantaneous terms. Numerically, one can observe that the system has a family of periodic orbits for delays between τ ≈ 0.7 and τ ≈ 1.3. The period and the phase portrait projected onto the (x, y)-plane of these orbits are shown in Figure 1(a) and (b). The family repeats, the orbits keeping their shape, for every integer N if we change the horizontal axis in Figure 1(a) to N T (τ ) + τ . The transformation N T (τ ) + τ is not just a parallel shift, since the dependence T (τ ) is, generically, nontrivial (in the example dT / dτ > 0). This leads to an overlapping of the families for large τ as shown in Figure 1(c) and (d).
Specifically, let us consider a regular periodic orbit x 0 , existing for a fixed value τ 0 of τ (regularity means that the unit Floquet multiplier of x 0 is simple, see Section 2). Then x 0 persists for τ in a small neighborhood of τ 0 such that we have a branch of periodic orbits depending locally on the parameter τ close to τ 0 . Let U be a small neighborhood of this orbit along the branch. Along this branch the period T of the orbit is a smooth function T (τ ) of τ . In Figure 1(a) an example orbit x 0 is indicated by a grey square at τ 0 = 1.081 and its neighborhood U is highlighted by a slight thickening of the line. If we assume that T (τ 0 ) = 0 (which is a genericity condition) then T (τ ) will have a slope uniformly bounded away from zero for all τ in the neighborhood corresponding to U. This implies that the images of the neighborhood U are stretched proportionally to N under the transformation τ → N T (τ ) + τ . Thus, for any given large τ a large number of periodic orbits from U coexist, and the number of coexisting orbits is proportional to τ (see the overlapping images of U in Figure 1(c) and (d)).
The next question to ask is: which of those many coexisting periodic orbits for large delays are dynamically stable? Is it possible to derive sharp stability criteria for the large-delay orbits based on quantities independent of the delay? More precisely, what is the stability of a given periodic orbit x 0 (such as the one indicated by a grey square in Figure 1(a) for delays τ = N T (τ 0 ) + τ 0 as N goes to infinity (this would be the sequence of orbits indicated by grey squares in Figure 1(c,d)).
Our example suggests that the stability for τ = N T (τ 0 ) + τ 0 where N → ∞ can indeed be determined from spectral properties of the periodic orbit at the small delay τ 0 . Figure 1(e) shows the spectrum of the periodic orbit x 0 highlighted by a grey square in Figure 1(a) for N = 25. One can see that, first, its Floquet exponents form bands, and, second, there is a large number of Floquet exponents very close to the imaginary axis (note the scale of the horizontal axis in Figure 1(e)). Figure 1(f) illustrates one of the results of this paper: after rescaling the real part, Floquet exponents converge to curves for N → ∞. These curves, the asymptotic continuous spectrum are computable by solving regular periodic boundary value problems parametrized by ω, the vertical axis in Figure 1(f). Since the original nonlinear DDE is autonomous, one of the curves of the asymptotic continuous spectrum touches the imaginary axis. Similar to the equilibrium case, we establish that the asymptotic continuous spectrum (together with the strongly unstable spectrum, see Section 2) determines the stability of the periodic orbit x 0 for sufficiently large N .
We have colored the branch in Figure 1(a) already according to the conclusions from the asymptotic spectra. The part that is displayed as a black curve in Figure 1(a) is the part of the branch that will be exponentially stable as N tends to infinity, whereas the grey part will be exponentially unstable, having a large number of weakly unstable Floquet exponents.
In short, the observation by Yanchuk & Perlikowski [16] implies that, if we find a periodic orbit x 0 of the DDE (1) for a fixed small delay τ 0 , then (under some genericity conditions) the DDE (1) has a large number of similar orbits coexisiting for any sufficiently large delay τ . This paper establishes a sharp critierion (again under some genericity conditions) determining if all of these coexisting periodic orbits are dynamically stable. It does so by providing formulas for the Floquet exponents of the periodic orbit with delay τ = N T (τ 0 ) + τ 0 that are independent of N but valid asymptotically for large N . For the example DDE (2)-(3), our results imply that the system has a large number of coexisting stable periodic orbits for every sufficiently large delay τ (as suggested by the Figures 1(c,d) where the black parts of the branches are stable).
Section 2 gives a non-technical overview of the results gradually developed and proven in the later sections. One central part of our paper is the construction of a characteristic function the roots of which are the Floquet exponents of the periodic orbit for delay τ 0 + N T (τ 0 ), and for which we can study the limit N → ∞. This construction is given in Section 3. The existence of this function h permits us to follow the approach from [5] and to extend their techniques to the case of periodic solutions. In the following sections 4 and 5 we describe two parts of the Floquet spectrum that show a different scaling behavior for large N (and, thus, τ ). The strongly unstable spectrum, converging to a finite number of asymptotic Floquet exponents which are determined by the instantaneous terms, is investigated in Section 4. Then, in Section 5, we analyze the Floquet exponents given by the asymptotic continuous spectrum shown in Fig 1(f). Based on these results, we can then prove in Section 6 our main result, a criterion about asymptotic stability based on the location of the asymptotic continuous and strongly unstable spectrum.
In contrast to the case of spectra at equilibria, where the asymptotic continuous and strongly unstable spectrum in many cases can be calculated explicitly (see [5]), the corresponding parts of the Floquet spectrum of a periodic orbit x 0 , can typically only be computed numerically. This limitation is not specific to our results but applies equally to most stability results based on the Floquet spectrum for ODEs. Even the periodic orbit itself is in most cases only computable with numerical methods. An exception is the case, where the periodic orbit is at the same time a symmetry orbit of the system. In this case, examples of asymptotic continuous and strongly unstable Floquet spectrum have been calculated explicitly in [15,17].
Conversely, the results presented in this paper provide an approach to approximate Floquet spectra numerically for large delays τ . If one uses numerical methods on problems with large delays, one faces the difficulty that the size of the matrix arising in the discretized eigenvalue problem grows not only with the desired accuracy (which is natural) but also with τ , even if the period of the orbit remains bounded. This is the case for the numerical methods used in DDE-Biftool [7]. This increase is to be expected because the number of Floquet exponents close to the imaginary axis increases with τ (see Theorem 6). In contrast to this, the asymptotic spectra can be computed with the same numerical method as in DDE-Biftool and a matrix size that is uniformly bounded for large τ . This paper does not discuss the details of the numerical computation of asymptotic spectra. However, our construction of the characteristic function h is uniform with respect to τ and could in principle be implemented numerically. In practice, it is better to apply the same analysis as is done in this paper to the large discretized eigenvalue problem.
In Fig. 1(f) we demonstrated the good agreement between original Floquet exponents and asymptotic spectra already in the situation of a moderately large delay, where the original Floquet exponents are still computable. Moreover, our stability criterion, which can be reliably verfied in this way, is valid independently on the actual value of the large delay.

Basic concepts and overview of the results.
PERIODIC ORBITS AND LARGE DELAY 5 2.1. Periodic orbits, stability, and Floquet exponents. Let x * be a periodic orbit of the n-dimensional autonomous nonlinear delay differential equation (DDE) that is, x * (t) satisfies (4) for all times t and has period T : x * (t) = x * (t + T ) for all t ∈ R. Without loss of generality we may assume that T = 1 (this can be achieved by a rescaling of time and the delay τ ). If the period T is equal to 1 then x * is also a periodic orbit ofẋ where N is a natural number and τ ∈ [0, 1). We denote the restriction of the periodic function We are concerned with the question whether the periodic orbit x * is stable or unstable for sufficiently large N in the following sense: Definition 1 (Exponential orbital stability and instability). Let X(t; ·) be the semiflow on C([−τ − N, 0]; R n ) induced by DDE (5). The periodic orbit x * is called exponentially orbitally stable if there exists a decay rate γ > 0 such that all initial history segments x 0 in a neighborhood of x * satisfy for some time shift t 0 and some constant C ≥ 0.
Similarly, x * is called exponentially unstable if there exists a growth rate γ > 0, a neighborhood N of x * and a constant C > 0 such that one can find initial history segments x 0 = x * arbitrarily close to x * that satisfy X(n; x 0 ) − X(n; x * ) ∞ ≥ C exp(γn) x 0 − x * ∞ > 0 for all n ∈ N as long as X(n; x 0 ) stays in the neighborhood N . This is the standard definition for stability of periodic orbits used also for ODEs except that the phase space is C([−τ − N, 0]; R n ) instead of R n . The notation · ∞ refers to the usual maximum norm in C([−τ − N, 0]; R n ).
Textbook theory of delay equations reduces the stability problem to the problem of finding eigenvalues of the linear map M N : which is given as the time-1 map of the linear DDĖ where the time-dependent n × n-matrices A(t) ∈ R n×n and B(t) ∈ R n×n are the partial derivatives of f in x * : [2]. If the right-hand side f of the nonlinear problem (4) is smooth in its arguments then the matrices A and B are also smooth periodic function of time t in the interval [0, 1].
Since the time-1 map (M N ) N +1 is compact the spectral theory for compact operators and the polynomial spectral mapping theorem imply that the spectrum σ(M N ) consists of a sequence of eigenvalues of finite multiplicity accumulating only at zero (zero is the only element of σ(M N ) that is not in the point spectrum). Also, λ = 1 is always an eigenvalue of M N becauseẋ * (t) satisfies (6) and has period 1.
The periodic orbit x * is exponentially (orbitally) stable if and only if (1) the eigenvalue 1 of M N is algebraically simple, and (2) all other eigenvalues of M N have modulus less than 1.
Similarly, it is exponentially unstable if at least one eigenvalue has modulus greater than 1.
Thus, the stability of iterations of M N is determined by its eigenvalues. We also use the term orbitally stable for the map M N , meaning that M N satisfies both of the above conditions.
The state space of DDE (6) is the function space C([−τ − N, 0]; R n ). Thus, initial value problems for (6) require specifying an infinite-dimensional initial condition. Similarly, one expects that a boundary value problem for a DDE requires the specification of infinitely many boundary conditions. However, periodic boundary value problems are easier to formulate: for example, a solution of the periodic boundary value problem for the general DDE (4) for period T = 1 is simply a function x ∈ C 1 ([−1, 0]; R n ) satisfyinġ and k ∈ Z is the integer part of τ − t. Since x is continuous and satisfies the periodicity condition (8) it can be extended continuously to a continuous function on the whole real line by defining x(s) = x(s mod[−1,0] ). Consequently, the right-hand side of (7) is continuous for all t ∈ [−1, 0], which guarantees that x can really satisfy the differential equation pointwise and is an element of C 1 . The solution x then automatically satisfiesẋ(−1) =ẋ(0), and, thus, by induction is as smooth as the right-hand side f . In this respect, periodic boundary value problems for DDEs are similar to boundary value problems for ODEs. In Section 3 we will reduce linear periodic boundary value problems for DDEs to low-dimensional linear systems of algebraic equations. Floquet exponents of M N can be found as those complex numbers µ for which the periodic boundary value probleṁ has a nontrivial solution y ∈ C 1 ([−1, 0], C n ) [2]. Note that in (9) we use the delay τ ∈ [0, 1) to calculate the Floquet exponents of M N , i.e. for a periodic orbit of system (6) with delay τ + N . Only the factor exp(−N µ) in front of B(t) accounts for the large delay whereas we have just extended y periodically for arguments less than −1. If (9) has a non-trivial solution y(t) for µ ∈ C then it has the non-trivial solution y(t) exp(2πikt) for µ + 2kπi for any integer k. Hence, we choose the Floquet exponent µ such that its imaginary part is between [−π, π).

2.2.
Asymptotic spectra for N → ∞. The set of Floquet exponents, Σ N , forms a discrete subset of the complex plane, the point spectrum of exponents, which depends on N . In order to describe in which form Σ N has a limit for N → ∞ we introduce two asymptotic spectra, which are also subsets of the complex plane. The notation follows [16].

Definition 3 (Instantaneous and strongly unstable spectrum). The set Σ
has a non-trivial solution y ∈ C 1 ([−1, 0]; C n ) is called the instantaneous spectrum.
The subset A + ⊆ Σ A of those µ with positive real part is called the strongly unstable asymptotic spectrum.
The instantaneous spectrum Σ A contains exactly n elements with imaginary part in [−π, π), counting algebraic multiplicity. We note that Σ A and A + do not depend on N but only on A. One result of our paper is that all Floquet exponents of M N with a real part that is positive uniformly in N converge to elements of the strongly unstable spectrum A + .
Yanchuk & Perlikowski [16] observed that the presence of strongly unstable spectrum is not the only possible cause of instability for large N . They observed that large numbers of Floquet exponents form bands that have a distance of order 1/N from the imaginary axis and have a spacing of order 1/N along the imaginary axis. In the limit N → ∞ these bands form curves after a rescaling of the real part by N . The limiting curves, called asymptotic continuous spectrum in [16], were defined with the help of a parametric periodic boundary value problem: has a non-trivial solution y ∈ C 1 ([−1, 0] R n ) for some ϕ ∈ R. The quantity ϕ is called the phase corresponding to γ + iω.
Again, the asymptotic continuous spectrum does not depend on N but only on A, B and τ .

2.3.
A characteristic function for Floquet exponents. We will reduce now the study of the various spectra and their relations to each other to a root-finding problem of a holomorphic function. The following lemma states the existence of a characteristic function for Floquet exponents that at the same time can be used to describe the asymptotic spectra A + and A c .
Lemma 5 (Characteristic function). There exists a function h : Ω 1 × Ω 2 ⊆ C × C → C which is holomorphic in both arguments with the following properties:

µ is in the instantaneous spectrum Σ A if and only if
and, hence, µ is in the strongly unstable asymptotic spectrum A + if h(µ, 0) = 0 and Re µ > 0, The algebraic multiplicity of µ as a Floquet exponent in the statements 1 and 2 equals its multiplicity as a root in (15) and (16).
From property 3 the motivation behind the name asymptotic continuous spectrum becomes clear: if we have a value γ 0 + iω 0 ∈ A c and the corresponding phase ϕ 0 , and ∂ 2 h(iω 0 , exp(−γ 0 − iϕ 0 )) is non-zero (which is generically the case) then a whole curve γ(ω) + iϕ(ω) satisfies h(iω, exp(−γ(ω) − iϕ(ω))) = 0 for ω ≈ ω 0 . These curves are the bands of the asymptotic continuous spectrum. Note that from the existence of the trivial Floquet exponent µ = 0 we can conclude that h(0, 1) = 0, which in turn implies that γ = ω = 0 with phase ϕ = 0 is in A c . In the generic case where the trivial exponent is contained in a single curve we call it the critical branch of A c .
The details of the construction of h, which modifies the general characteristic matrices and functions for periodic delay equations from [10,13], will be given in section 3. The general characteristic function constructed by [10,13] may have poles in the complex plane. The modification in Section 3 ensures that these poles of h(µ, z) stay in the left half-plane. Hence, the domain Ω 1 × Ω 2 of h contains all µ and z satisfying Re µ > −R and |z| < exp(R) where R > 0 is arbitrary but has to be chosen a-priori. Accordingly, statements 1-3 of Lemma 5 are valid only if both arguments of h satisfy their respective restriction. However, this is the case in the region of interest for stability and bifurcations.
The introduction of the characteristic function h clarifies how the different spectra can be calculated and reduces the analysis of the spectra to a root-finding problem. After defining h properly one could even use h to define the corresponding spectra by the properties listed in Lemma 5.

Main results.
With the help of the asymptotic spectra A + and A c we can formulate now a sharp criterion for the exponential orbital stability and instability of M N , which is our main result: Theorem 6 (Stability/Instability). The map M N (and, hence, the periodic orbit x * ) is exponentially orbitally stable for all sufficiently large N if all of the following conditions hold: S-1 (No strong instability) all elements of the instantaneous spectrum Σ A have negative real part (this implies in particular that the strongly unstable spectrum A + is empty), and S-2 (Non-degeneracy) ∂ 2 h(0, 1) = 0, that is, the Floquet exponent 0 is simple for sufficiently large N , and S-3 (Weak stability) except for the point µ = 0 with phase ϕ = 0 the asymptotic continuous spectrum A c is contained in {z ∈ C : Re z < 0}. The map M N is exponentially unstable for all sufficiently large N if one of the following conditions holds U-1 (Strong instability) the strongly unstable spectrum is non-empty, or U-2 (Weak instability) a non-empty subset of the asymptotic continuous spectrum A c has positive real part.
The weak stability condition S-3 is equivalent to stating that for all ω ∈ [−π, π) the function z → h(iω, exp(−z)) has no roots with non-negative real part (with the exception of z = 0 for ω = 0). Similarly, the weak instability condition U-2 is equivalent to stating that h(iω, exp(−z)) = 0 for some ω ∈ [−π, π) and some z ∈ C with positive real part.
Several additional corollaries follow from our analysis: Decay rate and dominant frequency. If M N satisfies S-1-S-3 and is, thus, exponentially stable for all sufficiently large N then the decay rate is at most (and generically) of order O(N −3 ) and the dominant relative frequency is of order O(N −1 ). That is, the dominant non-trivial Floquet exponents are a complex pair of the form ]. One branch of the asymptotic continuous spectrum touches the imaginary axis in a (generically quadratic) evenorder tangency at 0, and the dominant non-trivial Floquet exponent lies on this branch and next to the tangency.

Robustness.
A map M N that satisfies S-1-S-3 and is, thus, exponentially stable for sufficiently large N according to Theorem 6 remains exponentially stable for all sufficiently large N under all perturbations to A, B and τ of size less than some > 0. This does not depend on N because the quantities determining the exponential stability do not depend on N but only on A, B and τ . This means that periodic orbits that are stable for large N are uniformly robust with respect to perturbations of the system despite their weak attraction rate of order N −3 .
This permits the conclusion about the coexistence of stable periodic orbits in the example in Figure 1. The asymptotic spectra of orbits in the small neighborhood U are small perturbations of the asymptotic spectrum shown in Figure 1(f). In particular the curve touching the imaginary axis will also be curved to the left for all periodic orbits in U.
Spectral approximation. Let the instantaneous spectrum Σ A have a positive distance to the imaginary axis. Then we can guarantee that certain points z in the complex plane are in the resolvent set of M N (that is, the periodic boundary value problem (9)-(10) has only the trivial solution for µ = z): • every point z in the positive half-plane that is not in the strongly unstable spectrum A + is in the resolvent set of M N for sufficiently large N . • All points of the form γ/(N + τ ) + iω are in the resolvent set of M N for sufficiently large N if the point z = γ + iω is not in the asymptotic continuous spectrum A c . These two statements about the resolvent set of M N imply that the spectrum of M N must be close to A + or (after rescaling) close to A c . The other direction also holds if the instantaneous spectrum Σ A is not on the imaginary axis: • if µ ∈ A + has multiplicity k then k Floquet exponents of M N converge to µ for N → ∞ (counting multiplicity). • Let γ * + iω * be in the asymptotic continuous spectrum A c , that is, for some phase ϕ * ∈ [−π, π). Then we will find Floquet exponents µ N of M N that satisfy (note that this is small-o) Estimate (17) is rather weak. We need and prove a stronger and more detailed estimate under the additional condition that ∂ 2 h(iω * , exp(−γ * − iϕ * )) = 0. Then we have a local root curveγ(ω) + iφ(ω) satisfying h(iω, exp(−γ(ω) − iφ(ω))) = 0 for ω near ω * , and for sufficiently large N we find algebraically simple Floquet exponents µ k of M N satisfying where k are integers such that 2kπ/(N + τ ) is near ω * , and the O-s are smooth functions of ω. These Floquet exponents µ k form a band of discrete complex numbers approximating the curveγ(ω) + iφ(ω) of asymptotic continuous spectrum.
Degeneracies. Theorem 6 is sharp except for several degenerate cases that are excluded by the conditions S-1-S-3 and U-1-U-2. Degeneracies limiting the region of stable periodic orbits are: • An element of the instantaneous spectrum Σ A has zero real part.
• The partial derivative ∂ 2 h(0, 1) equals 0. In this case the periodic orbit cannot be exponentially stable. • (Turing and long wavelength instability) The asymptotically continuous spectrum A c touches the imaginary axis in a quadratic tangency at some value ±iω 0 . In this case the orbit is not exponentially unstable and it is still possible that the periodic orbit is exponentially stable for all N (this is different from the stationary case discussed in [5]) but it may also be weakly stable or unstable. A special case is that A c touches the imaginary axis in a quadratic tangency in the point µ = 0 with phase ϕ = π. In this case the periodic orbit is still exponentially stable for large N . • (Modulational instability) The critical branch of the asymptotic continuous spectrum A c , containing the trivial exponent µ = 0, has a tangency with the imaginary axis at ω = 0 that is of higher order than quadratic. In this case the periodic orbit is still exponentially stable according to Theorem 6.
The characteristic function h constructed in Section 3 is still valid for the degenerate cases. However, for a detailed discussion of these degeneracies, one needs not only to specify a defining equation for each degeneracy (this is straightforward) but one also has to state secondary non-degeneracy conditions. In particular the case of instantaneous spectrum with zero real part is somewhat subtle, even though it looks similar to the other listed degeneracies of co-dimension 1 at first sight. In the analogous situation for the spectrum of equilibria, it turns out that instantaneous spectrum iω 0 with zero real part generically implies a singularity of the asymptotic continuous spectrum, γ(ω 0 ) = ∞, and hence is not part of the stability boundary. Due to these reasons, we believe that a comprehensive treatment of the degeneracies is beyond the scope of the present paper. We introduce the complex variable z ∈ C and consider the periodic boundary value problem for t ∈ [−1, 0] In a first step we will construct a characteristic matrix ∆(µ, z) for (20)-(21) such that the roots of its determinant h(µ, z) = det ∆(µ, z) will be precisely those pairs of points (µ, z) in some subdomain of C × C for which (20)-(21) has a nontrivial continuously differentiable solution y ∈ C 1 ([−1, 0]; C n ). Thus, by inserting z = exp(−(N + τ )µ) we will then obtain the characteristic function h N (µ) in a subdomain of C. Consider a partition of the periodicity interval [−1, 0] into k intervals of size 1/k: Using this partition we formulate a multiple initial value problem (MIVP) for a vector of k initial (or restart) values (v 0 , . . . , v k−1 ) T ∈ C nk (similar to multiple shooting):ẏ where t ∈ [−1, 0] and z ∈ C. Notice that in (23) the solution y on the interval [t j , t j+1 ) depends on the solution y(t) in other intervals due to the term y(t − τ ) mod[−1,0] in the right-hand side of (23).
The main purpose of this partition is to reduce the length of the integration interval at the cost of increasing the dimension of the system. Indeed, this construction is very similar to the construction in [11], where in the case of rational τ a reduction to an equivalent system of ODEs could be achieved, and for the case of irrational τ corresponding rational approximations have been considered. This can be seen by introducing u j (t) := y(t + t j ), which satisfy the system of equationṡ with j = 0, . . . , k − 1. This system can now be considered as an initial value problem on the interval [0, 1/k), and for a solution on this smaller interval [0, 1/k) each component u j (t) represents the solution y(t + t j ) in the corresponding subinterval I j .
Note that instead of the delayed term u j The initial value problem for the coupled system (25)-(26) is an equivalent formulation of the multiple initial value problem (23)- (24) and, in a similar way as the boundary value problem (20)-(21), contains both delayed and advanced arguments. If τ is rational then (25)-(26) is a system of ODEs. This was the starting point for the rational approximations in [11]. We take another route by showing directly that the Picard-Lindelöf iteration for (25)-(26) converges, similar to initial-value problems of ODEs. In the sequel, we will clarify this point, using the formulation (23)- (24) which is more convenient for our purposes.
Let us denote by U (t, s, µ) ∈ R n×n the propagation matrix of the linear ODE defining the instantaneous spectrum, (11). That is, The norm of U (t, s, µ) can be estimated by In (27) we have used the notation A ∞ = max t∈[−1,0] A(t) ∞ . We will use the same notation for B. In order to clarify in which sense system (23)- (24) is an initial value problem and what it means for y to be a solution of system (23)-(24) we formulate an integral equation which is equivalent to (23)- (24): The right-hand side of the integral equation (28) is an affine map, mapping C k back to itself such that (28) is of the form where S(µ) : C nk → C k is defined by (29), and L k (µ)y is the linear part of the right-hand side in (28) (that is, the integral term). The linear map L k (µ) : C k → C k is continuously differentiable (and, thus, holomorphic) with respect to the complex variable µ. A simple estimate for the norm of L k with respect to the · ∞ -norm gives us the unique solvability of the fixed point problem (32) (which is actually the integral equation (28)-(30)):

Lemma 7 (Existence and Uniqueness of solutions for IVP).
Let R > 0 be arbitrary. If we set the number of sub-intervals, k, such that then the affine integral equation (28)-(30) has a unique solution y ∈ C k for all µ and z satisfying Re µ ≥ −R, |z| ≤ exp(R) and all tuples (v 0 , . . . , v k ) T ∈ C nk .

PERIODIC ORBITS AND LARGE DELAY 13
The details of the norm estimate for L k are in Appendix A. The solution y can be written as y = [I − zL k (µ)] −1 S(µ)v. When is this solution y (which is only in C k for general v) continuously differentiable on the whole interval? This requirement is a linear condition on the tuple v. Let us fix a constant R > 0 and choose the number of sub-intervals k > C(R).
Note that the integral equation (28) implies that v j = lim t tj y(t). In this sense the values v j are the initial (or restart) values for the differential equation (23). The construction of ∆ gives a well-defined matrix for all z and µ satisfying Re µ ≥ −R and |z| ≤ exp(R): for a given tuple v we evaluate the unique solution y = [I − zL k (µ)] −1 S(µ)v ∈ C k of the integral equation (28) and then we use this y to evaluate the right-hand side of (35). The characteristic matrix ∆ is set up such that for v in the kernel of ∆(µ, z) the solution y also satisfies the differential equation (23) on the whole interval (and does not have jumps) including the periodic boundary conditions: Proof. The rows 2 to k of the right-hand side in the definition (35) of ∆ ensure continuity of y: since v j = lim t tj y(t) for j = 1, . . . , k − 1 the condition that these rows are equal to zero reads lim t tj y(t) = lim t tj y(t), and (by inserting the right-hand side of (28)) v j+1 = U (t j+1 , t j , µ)v i + z The first row of the condition ∆(µ, z)v = 0 reads v 0 = y(0), which makes sure that y is periodic at the boundary of [−1, 0], and guarantees that the integrand in (38) is continuous at s = τ −1. Consequently, the integrand is continuous everywhere, which implies that y is continuously differentiable. Thus, we can differentiate (38) with respect to t, which implies that y satisfies the differential equation (20). This, in turn, implies that alsoẏ(−1) =ẏ(0) because the right-hand side of (20) is periodic.
Since the map L k depends analytically on µ, the matrix ∆ depends analytically on µ and z, too (as long as Re µ ≥ −R and |z| ≤ exp(R)). Calling ∆ the characteristic matrix makes sense for the following reason: The algebraic multiplicity of µ as a Floquet exponent of M N is equal to the order of µ as a root of h N .
Proof. That µ is a Floquet exponent of (6) if and only if ∆(µ, exp(−(N + τ )µ)) has a non-trivial kernel is clear because Floquet exponents of M N were defined as those complex numbers for which (20)-(21) has a non-trivial solution. The statement about the multiplicity of µ follows from arguments similar to [3], which are laid out in detail in [10] (an equivalence between the eigenvalue problem for M N , exp(µ)v − M N v = 0, and the algebraic equation ∆(µ, exp(−(N + τ )µ))v = 0 is constructed in [10]).
Both functions, h and h N , are real analytic (that is, their expansion coefficients are real numbers since A, B and τ are real).
For any given R > 0 we have chosen k > 0 and constructed a characteristic function h N (µ), defined on the half plane the roots of which are precisely the Floquet exponents (counting multiplicities) of the period map M N for the linearized delay differential equation (6). As we are interested in the location of Floquet exponents of the time-1 map M N to the right or close to the imaginary axis we have reduced the eigenvalue problem to a study of the asymptotic behavior of roots of the holomorphic function h N for N → ∞. Furthermore, the asymptotic spectra are defined as roots and root curves of h by construction of h: µ > 0 is in the strongly unstable spectrum A + if and only if h(µ, 0) = 0, and µ = γ + iω is in the asymptotic continuous spectrum A c if and only if h(iω, exp(−γ − iϕ)) = 0 for some phase ϕ ∈ R.
This proves Lemma 5 and allows us to follow an approach similar to [17,5] in the following sections. Note that in the case where the period and the delay are rationally dependent, the existence of a suitable characteristic function follows much easier (see [11]), but the structure of the asymptotic spectra discussed below will be not affected. 4. Strongly unstable spectrum. One immediate consequence of Lemma 5 is the statement asserted in [16] for the relation between eigenvalues of M N and the strongly unstable spectrum A + . Note that µ is in the strongly unstable spectrum A + if and only if h(µ, 0) = 0: if one sets z = 0 in (20) the differential equation reduces to (11), the expression defining A + . For large N unstable Floquet exponents of M N either approach the imaginary axis or A + : Lemma 11 (Convergence to strongly unstable spectrum). Let Re µ 0 > 0. If h(µ 0 , 0) = 0 then there exists a N 0 such that µ 0 lies in the resolvent set of the time-1 map M N for all N > N 0 .
If µ 0 is a root of multiplicity k of the function h(·, 0) then every sufficiently small neighborhood U of µ 0 contains exactly k Floquet exponents (counting multiplicity) of M N for all N > N 0 (N 0 depends on U ).
The resolvent set is the set of all complex numbers µ for which the map exp(µ)I − M N is an isomorphism (characterized by h N (µ) = 0 for Re µ ≥ −R/(N + τ )).
Lemma 11 shows that M N is exponentially unstable for all sufficiently large N if the strongly unstable spectrum is non-empty, which is condition U-1 of Theorem 6.
We assume for the remainder of the section that the instantaneous spectrum Σ A has a positive distance to the imaginary axis. The idea behind the construction of the asymptotic continuous spectrum is that for Floquet exponents µ close to the imaginary axis, that is, for µ of the form with a bounded factor γ in the real part, the roots of the characteristic function h N (µ) = h(µ, exp(−(N + τ )µ)) converge to a regular limit for N → ∞ after inserting the scaling (40). This can be made more specific: if the instantaneous spectrum Σ A has a positive distance from the imaginary axis then all Floquet exponents that are not converging to the strongly unstable spectrum have a real part less than R 2 /(N + τ ) where R 2 does not depend on N : Lemma 12 (Convergence to the imaginary axis). Assume that the elements of the instantaneous spectrum Σ A have a positive distance 2 > 0 from the imaginary axis. Let us denote the points in the strongly unstable spectrum, A + , by µ 1 ,. . . , µ j . There exists a constant R 2 ≥ 0 such that all Floquet exponents µ of M N satisfy either µ ∈ B (µ j ) for some j, or Proof. For all complex numbers µ that have non-negative real part and are outside of the balls B (µ j ) the matrix I − U (0, −1, µ) is invertible and the inverse has a uniform upper bound (remember that U is the propagation matrix ofẏ = [A(t) − µI]y): We also know that µ is a Floquet exponent if the integral equation (38) (which is equivalent to the differential equation (37)) has a non-trivial periodic solution y(t) for z = exp(−(N + τ )µ). This implies (y(−1) = y(0) and, hence, The factor of z in the right-hand side of this fixed-point problem for y is a linear operator K(µ) on the space of continuous functions. K(µ) is uniformly bounded for all µ that have non-negative real part and are outside of the balls B (µ j ). Let us denote an upper bound on the norm of K(µ) by C 1 . Then for all z satisfying |z| < C −1 1 the integral equation (38) cannot have a non-trivial solution. Since, for Floquet exponents of M N , z equals exp(−(N + τ )µ) this means that µ has to satisfy | exp(−(N + τ )µ)| ≥ C −1 1 , and, thus Re µ ≤ log C 1 N + τ for |z| to be larger than C −1 1 . Consequently, if we choose R 2 > log C 1 then a µ that has non-negative real part and lies outside of the balls B (µ j ) must satisfy (41) in order to be a Floquet exponent of M N .
A corollary of the construction described by equation (42) is that for each ω the intersection of the asymptotic continuous spectrum A c with the horizontal line {µ ∈ C : Im µ = ω} can be described as a set of eigenvalues of a compact operator K(iω). Thus, if Re Σ A = 0 the asymptotic continuous spectrum forms curves parametrizable by ω, which may only have poles at −∞: Corollary 13 (Asymptotic continuous spectrum consists of curves). Assume that the instantaneous spectrum Σ A has a positive distance to the imaginary axis. A point µ is in the asymptotic continuous spectrum if exp(µ) is an eigenvalue of the compact linear operator K(iω) : In particular there exists a constant R 3 ≥ 0 so that the asymptotic continuous spectrum lies to the left of the vertical line {z ∈ C : Re z = R 3 }.
The matrix U in the definition of K is the monodromy matrix of the instantaneous problemẋ(t) = [A(t) − iω]x(t). The inverse of I − U (0, −1, iω) exists because of our assumption that the instantaneous spectrum, Σ A , does not contain points on the imaginary axis. The corollary implies that the asymptotic continuous spectrum consists of curves even in points where the regularity condition on ∂ 2 h = 0 is violated. In these points K(iω) has eigenvalues of finite multiplicity larger than 1 such that finitely many curves of the asymptotic continuous spectrum cross each other (or are on top of each other).
All Floquet exponents of M N that do not converge to the strongly unstable spectrum for N → ∞ (these can be at most n, equaling the dimension of A(t)) are either stable, or they satisfy restriction (41) on the upper bound. Thus, apart from the strongly unstable spectrum, the multipliers, which are of interest from the point of view of stability and bifurcations, lie in the strip and have the form (after shifting them into the strip {z : Im z ∈ [−π, π)} by subtraction of an integer multiple of 2πi). As we focus our discussion on Floquet exponents of the scale (43) from now on it makes sense to introduce the notation of a scaled Floquet exponent. We can now formulate precisely in which sense the asymptotic continuous spectrum A c is the limit of the spectra of M N . First, we make a statement about resolvent sets.
Lemma 15 (Points distant from A c ). Let the instantaneous spectrum Σ A have a positive distance to the imaginary axis. If γ 0 > −R and γ 0 + iω 0 is not in the asymptotic continuous spectrum A c then γ 0 + iω 0 is in the scaled resolvent set of M N for sufficiently large N .
Lemma 15 implies that all spectrum of M N near the imaginary axis has to be close to A c even after blow-up of the real part. The other direction, that every point of A c in the strip −R < Re µ < R 2 is approached by scaled Floquet exponents of M N , is also true (again under the assumption that the instantaneous spectrum Σ A is not on the imaginary axis). We prove this direction first for regular root curves of h because this is the most common case, and the missing piece for the stability criterion. Our estimate is slightly sharper than mere approximation to make it useful for our proof of the stability criterion.
In short, Lemma 16 below states that any regular curve of asymptotic continuous spectrum,γ(ω) + iω (with its corresponding phaseφ(ω)), is approximated by scaled Floquet exponents of the form γ k + iω k where k are integers such that 2kπ/(N + τ ) is near ω, and the O-s are smooth functions of ω. This proves the stronger estimate (18)-(19), given in the non-technical overview in Section 2.
For all N ≥ N * there are scaled Floquet exponents µ N,k of M N near γ * + iω * . All scaled Floquet exponents µ N,k that are sufficiently close to γ * + iω * are algebraically simple and have the form where: • k is any integer satisfying • ω k is the unique solution of the fixed point problem for ω • the functionsγ N andφ N are perturbations ofγ(ω) andφ(ω) of the form where g(ω, ) is a smooth complex-valued function, which is independent of N and defined for ω ∈ (ω * − δ, ω * + δ) and ∈ [0, 1/N * ).
We notice that the scaled Floquet exponents of M N lie on bands: they are on the curve given byγ N (ω) + iω, which is a small perturbation of the curve of asymptotic continuous spectrumγ(ω) + iω. The spacing between scaled Floquet exponents along this band is given by the fixed point equation (47). The fixed point problem (47) is only weakly implicit since the right-hand side terms containing ω all have a pre-factor 1/(N + τ ), which is small. Hence, expression (45) for γ k follows immediately from (48) and expression (44) for ω k follows from (47) The root problem h (iω + γ, exp (−γ − iϕ)) = 0 is for small a small (order ) perturbation of the root problem (50). Thus, for sufficiently small , a root curve of the formγ (ω) + iφ (ω) exists for ω in some neighborhood of ω * , and it has the form where g(ω, ) is a smooth complex-valued function defined for ω in a small neighborhood of ω * and ∈ [0, max ) (with some max > 0). Note that the error term contains a factorγ(ω), making the error equal to zero on the imaginary axis. (See Appendix B for details of how to extract this factor from the error.) Inserting = (N + τ ) −1 and labeling the curvesγ (ω) asγ N (ω) andφ asφ N gives the definitions (48) and (49) in the lemma. Correspondingly, we make an initial choice for the minimal N , N * , as 1/ max . A point on the curveγ N (ω) + iω is a scaled Floquet exponent of M N if and only if its imaginary part ω satisfies exp(−iω(N + τ )) = exp(−iφ N (ω)) and, thus, for some integer k ∈ Z.
After dividing by N + τ this becomes the fixed point equation (47)  We choose a neighborhood U of ω * of the form (ω * − 2δ, ω * + 2δ) such thatφ N is well-defined for all ω ∈ U and all N ≥ N * , and satisfies |φ N (ω)| < L for some constant L (which is independent of N ). Next, we increase N * such that L N + τ < 1 for all N > N * , and ϕ N (ω) N + τ < δ for all N ≥ N * and ω ∈ U.
Then the right-hand side of the fixed point problem (47) is contracting with a Lipschitz constant L/(N + τ ) for all ω ∈ U and it is mapping U = (ω * − 2δ, ω * + 2δ) back into itself: Thus, the Banach Contraction Mapping Principle guarantees that (47) has a unique solution ω k for all integers k satisfying (52). Finally, we confirm the algebraic simplicity of the scaled Floquet exponents µ N,k =γ N (ω k ) + iω k by checking the multiplicity of the corresponding root of h N : the derivative of h N , divided by N + τ , is Insertingγ N (ω k )/(N +τ )+iω k for z on the right, and the relation exp(−i(N +τ )ω k ) = exp(−iφ N (ω k )) we get Since ω k is in the neighborhood U of ω * in which ∂ 2 h is non-zero the overall derivative is non-zero for sufficiently large N .
Lemma 16 is based on a perturbation argument assuming that the complex function h(iω, exp(−γ − iϕ)) has a regular root curve ω →γ(ω) + iφ(ω). Thus, it is also valid if the instantaneous spectrum, Σ A , does not have a positive distance to the imaginary axis as long as one restricts consideration to Floquet exponents of the form γ/(N + τ ) + iω (with bounded γ). Positive distance of Σ A to the imaginary axis merely ensures that all Floquet exponents µ with Re µ > −R/(N + τ ) are of this form (except for those approximating the strongly unstable spectrum A + ).
Lemma 15 and Lemma 16 about the approximation of the asymptotic continuous spectrum, together with Lemma 11 about the approximation of the strongly unstable spectrum, are the tools that we need to prove the criterion for asymptotic stability from Theorem 6. Before turning to asymptotic stability let us prove the remaining statement about spectral approximation of the asymptotic continuous spectrum A c . Let γ * + iω * be an element of A c with phase ϕ * (that is, h(iω * , exp(−γ * − iϕ * )) = 0).
Then we find scaled Floquet exponents of M N that approximate γ * + iω * even if the non-degeneracy condition ∂ 2 h(iω * , exp(−γ * − iϕ * )) = 0 is not satisfied: Lemma 17 (Approximation of A c ). Assume that the instantaneous spectrum Σ A is not on the imaginary axis. Let γ * ∈ [−R, R 2 ] and let γ * + iω * be in the asymptotic continuous spectrum A c . Then there exists a sequence γ N + iω N of scaled Floquet exponents of M N such that Lemma 17 covers the claim of Theorem 6 about exponential instability of M N for sufficiently large N under the condition of weak instability U-2.
Proof. Since the instantaneous spectrum Σ A does not contain points on the imaginary axis, exp(γ * + iϕ * ) is a non-zero eigenvalue of the compact operator K(iω * ) as introduced in Corollary 13. Consequently, exp(γ * + iϕ * ) is isolated and has finite multiplicity, which implies that γ * + iϕ * has finite multiplicity as a root of z → h(iω * , exp(−z)). Let us define k(N ) for large N as By construction of k(N ) we have that Since z → h(iω * , exp(−z)) has an isolated root at γ * + iϕ * , the functions also have roots z N = γ N + iϕ N for sufficiently large N which converge to γ * + iϕ * for N → ∞. Let us define Then, by construction, γ N + iω N is a scaled Floquet exponent of M N since Moreover, γ N → γ * and ω N → ω * , which proves the claim of the lemma.
6. Asymptotic stability for large delay. The convergence results for the strongly unstable spectrum in Lemma 11 and the asymptotic continuous spectrum in Lemma 15 and Lemma 16 can be combined to give a criterion for the stability of the periodic orbit x * of the original nonlinear system (5) depending on the triplet (A, B, τ ) that ensures stability of x * for all sufficiently large N . If the instantaneous spectrum Σ A has a positive distance to the imaginary axis then the point 0 is part of the asymptotic continuous spectrum because 0 is a Floquet exponent for all N . If ∂ 2 h(0, 1) = 0 then a regular curveγ(ω)+iω of the asymptotic continuous spectrum is passing through 0 (that is,γ(0) = 0). This curve at least touches the imaginary axis because γ(0) = 0 andγ (0) = 0 (γ(ω) is an even function, thus, all odd derivatives of γ are zero).
Proof. (Lemma 18) Since ∂ 2 h(0, 1) = 0 we know that the Floquet exponent 0 is simple for M N if N is sufficiently large. Also, since the instantaneous spectrum is in the negative half-plane the strongly unstable spectrum is empty. Hence, M N is the return map of a stable periodic orbit x * if it has no non-zero scaled Floquet exponent γ + iω for which γ ∈ [0, R 2 ] (and ω ∈ [−π, π)).
Proving the statement by contradiction, we assume that, for a sequence of increasing N , M N has a scaled Floquet exponent γ N + iω N = 0 where γ N ∈ [0, R 2 ].
Consequently, the accumulation point must be an element of the asymptotic continuous spectrum A c,ϕ . Since γ N ≥ 0 for all N of the sequence, γ * must be greater or equal 0, too. By assumption, the only element of A c,ϕ with non-negative γ is γ * + iω * = 0. Thus, γ * = ϕ * = ω * = 0. We choose > 0 sufficiently small such that we can apply Lemma 16 to (γ * , ϕ * , ω * ). We have that ω N ∈ (− , ) and γ N + iϕ N ∈ B (0) for sufficiently large N of the sequence (γ N , ω N , ϕ N ) (how large N has to be, depends on ). Since ∂ 2 h(0, 1) = 0 this guarantees that γ N lies on the curveγ N (ω) given in (48) in Lemma 16: Since for all sufficiently large N the factor in front ofγ(ω N ) is positive,γ N (ω N ) must have the same sign asγ(ω N ), which is negative if ω N = 0 due to the assumptions of the lemma. Since γ N is assumed to be non-negative, this implies that ω N = γ N = 0. Thus, the scaled Floquet exponents γ N + iω N of M N are zero for the converging sub-sequence, which is in contradiction to our assumption γ N + iω N = 0.
Lemma 18 proves the exponential stability claim of Theorem 6.

7.
Conclusions. We have shown that Floquet exponents of periodic solutions of delay differential equations (1) with large delay can be approximated by a set of continuous curves (asymptotic continuous spectrum) that are independent of the delay, and a finite set Floquet exponents (strongly unstable spectrum). Although the structure of the spectrum is shown to be similar to the case of equilibria [5], there are some unique features, which occur specifically for periodic orbits. Our results are based on the construction of the characteristic function, the roots of which give Floquet multipliers of the periodic orbit.
Using the asymptotic spectra we have been able to provide necessary and sufficient conditions for the exponential stability of periodic solutions for all sufficiently large delays. Our results are applicable to the case when the delay τ is large compared to the period T of the solution. In this case, the large parameter N , which controls precision of the asymptotic approximation is proportional to τ /T .
Let us mention some of the specific features of the spectrum. In contrast to the equilibrium case, the asymptotic continuous spectrum of Floquet exponents for periodic solutions contains generically a curve with a tangency to the imaginary axis (see Figure 1(e,f)). We have proved that even in the presence of this tangency, the stability (or instability) of the asymptotic continuous spectrum implies the exponential stability (resp. instability) of the corresponding periodic orbit. We have shown that the generic decay rate of perturbations of the exponentially stable periodic orbit of system (1) is of the order N −3 .
From the practical point of view, our results can be useful for studying periodic regimes in applications that involve feedback with large delays, for example, semiconductor lasers with optical feedback [4,17], or systems with feedback control [9].
From a mathematical point of view our result may provide a rigorous approach to proving the existence of a large number of stable rapidly oscillating periodic solutions for some special cases in which the asymptotic spectra can be computed explicitly. This would provide a contrast to the results for scalar feedback equations [14].