Multiple time-scale delay systems in mathematical biology and laser dynamics

(1)

Multiple Time-Scale Delay Systems

in Mathematical Biology and Laser Dynamics

vorgelegt von

M. Sc.

Stefan Ruschel

von der Fakultät II - Mathematik und Naturwissenschaften

der Technischen Universität Berlin

zur Erlangung des akademischen Grades

Doktor der Naturwissenschaften

Dr. rer. nat.

genehmigte Dissertation

Promotionsausschuss:

Vorsitzender:

Gutachterin:

Gutachter:

Prof. Dr. Wilhelm Stannat

Prof. Dr. Kathy Lüdge

Prof. Dr. Dmitry Turaev

PD Dr. Serhiy Yanchuk

Tag der wissenschaftlichen Aussprache: 29.07.2019

(2)

(3)

Abstract

This thesis addresses the dynamical behavior of Delay Differential Equations (DDEs) with a mul-tiple time-scale structure as a consequence of large delays, or additional small parameters. More specifically, it gives attention to the effects of delay-induced instabilities of equilibrium solutions as well as families of equilibrium solutions and discusses the corresponding nonlinear dynamical phenomena that are not direct consequences of finite-dimensional geometric theory for Ordinary Differential Equations. When studying such systems of DDEs close to equilibrium, one is first concerned with the spectral properties of a corresponding linearized system. Using an asymptotic approach, a rigorous description of the spectrum of linear DDEs with multiple hierarchical large delays is provided. It is shown that the spectrum splits into two distinct parts: the strong spectrum and the pseudo-continuous (or weak) spectrum reflecting the hierarchical structure of the delays. It is shown that a generic destabilization, a so-called weak instability, is mediated by a subset of the pseudo-continuous spectrum crossing the imaginary axis, that is associated with the largest delay in the system.

On the basis of three specific examples motivated from Mathematical Biology as well as Laser Dynamics, and applying the available invariant manifolds theory for semiflows, these results are then used to illuminate the nonlinear dynamical behavior of DDEs close to families of weakly un-stable equilibrium solutions. In order to demonstrate how this local information can be lifted to understand the global dynamics, a specific 2-delay epidemiological model is analyzed. Here, the dynamics away from the manifold can be studied independently of the center direction, and the case of weak-instability is studied in detail. The obtained results are then interpreted in the biological context and allow for valuable insight into consequences of delays and imperfect implementation of isolation in infectious disease management. Thereafter, the focus changes towards DDEs with a single large delay and an additional small parameter. Here, the obtained results are used to study reduced systems, where the slow time-scale has been eliminated and the problem is reduced to the local study of families of equilibrium solutions. In particular, it is shown how the weak instability of the “laser off” state of the Lang-Kobayashi laser model translates into weakly chaotic solutions bearing some analogy to the weak-strong splitting of the spectrum of the autonomous linear system. Finally, it is shown that geometric singular perturbation theory of DDEs has aspects that are funda-mentally different from singular geometric perturbation theory for Ordinary Differential Equations. A minimal model is studied, where delay induces switched states. Necessary conditions for their existence are derived, and the case of weak-instability of the corresponding family of equilibrium solutions is studied in detail.

(4)

(5)

Zusammenfassung

Gegenstand dieser Dissertation ist die Dynamik von Delay Differentialgleichungssystemen (DDEs) mit Zeitskalentrennung aufgrund großer Zeitverzögerungen, oder zusätzlicher kleiner Systempara-meter. Insbesondere werden die Auswirkungen verzögerungsinduzierter Instabilitäten auf Gleich-gewichtslösungen, oder Familien von Gleichgewichtslösungen aufgezeigt, und nichtlineare dynamische Phänomene untersucht, die nicht unmittelbar mit Mitteln der geometrischen Theo-rie endlich-dimensionaler gewöhnlicher Differentialgleichungen beschreibbar sind. Werden Sys-teme dieser Art auf ihre Gleichgewichtslösungen hin untersucht, so steht zunächst das Spektrum des entsprechenden linearisierten Systems im Vordergrund. Als eines der Hauptergebnisse dieser Arbeit wurde gezeigt, dass das Spektrum von linearen Delay Differentialgleichungen mit mehreren hierarchischen großen Zeitverzögerungen in zwei Teile zerfällt: das starke Spectrum und das pseudo-continuierliche (oder schwache) Spektrum. Eine generische Destabilisierung, eine in dem Fall so-genannte schwache Instabilität, wird durch eine Teilmenge des pseudo-kontinuierlichen Spektrums vermittelt, welche mit der größten Zeitverzögerung in Verbindung gebracht werden kann.

Anhand dreier spezifischer Beispiele aus der Mathematischen Biologie und der Laser Dynamik, und mithilfe der verfügbaren Theorie invarianter Mannigfaltikeiten von Halbflüssen, wurden diese Ergebnisse auf das Verhalten von Lösungen in der Umgebung von Familien von schwach insta-bilen Gleichgewichtslösungen angewendet. Um zu demonstrieren wie das Wissen um das lokale Verhalten auf die globale Dynamik übertragen werden kann, wird ein spezifisches epidemiolo-gisches Model mit zwei Delays untersucht. In diesem speziellen Fall, kann die Dynamik unab-hängig von der Richtung entlang der Mannigfaltigkeit betrachtet und der Fall schwacher Insta-bilität im Detail untersucht werden. Die Ergebnisse werden anschließend im biologischen Kontext interpretiert und erlauben mathematische Einsicht in die Auswirkungen von Zeitverzögerungen und Unzulänglichkeiten in der Isolation im Hinblick auf die Eindämmung ansteckender Krankheiten. Im Anschluss, richtet sich der Fokus der Arbeit auf DDEs mit großer Zeitverzögerung und einem kleinen Parameter. Die erbrachten Resultate werden hier angewendet um formale, reduzierte Sys-teme und die entsprechenden Familien von Gleichgewichtslösungen zu untersuchen. Insbeson-dere wird gezeigt, wie die schwache Instabilität der trivialen Lösung des Lang-Kobayashi Laser Models sich, in Analogie zum autonomen Fall, in den schwach-chaotischen Lösungen des Systems wiederfindet. Abschließend wird gezeigt, dass die singuläre Störungstheorie von DDEs Lösun-gen zulässt, die nicht mit Mitteln Gewöhnlicher DifferentialgleichunLösun-gen erklärt werden können, etwa sogenannte delay-induced switched states. Notwendige Bedingungen für die Existenz solcher Lösungen und der Fall schwacher Instabiltät von Gleichgewichtslösungen des entsprechenden re-duzierten Systems, werden im Detail untersucht.

(6)

(7)

Relevant publications

This thesis is based on research published in the following research articles:

1. Chap. 3: S. Ruschel and S. Yanchuk, “The Spectrum of Delay Differential Equations with

Multiple Hierarchical Large Delays”, arXiv:1902.00404 (to appear Discr. Cont. Dyn-S)

SR wrote, proofread and edited the text. SR performed rigorous analysis, provided the proof of Theorems 1 and 3, and generated all figures. SY advised on the investigation, proof-read and edited the text. SY drafted the proof of Theorem 2, and provided the algorithm to compute the spectrum.

2. Chap. 4: S. Ruschel, T. Pereira, L.-S. Young and S. Yanchuk, “An SIQ delay differential

equations model for disease control via isolation”, J. Math. Biol. 79, 249-279 (2019)

SR performed the rigorous analysis and investigation including the implementation of all computations; generated all figures; wrote, proofread and edited the text. TP and LSY de-signed the study, advised on the investigation, wrote, proofread and edited the text. SY contributed to the theoretical investigation, in particular Theorem 9, proofread and edited the text.

3. Chap. 4: L.-S. Young, S. Ruschel, S. Yanchuk and T. Pereira, “Consequences of delays and

imperfect implementation of isolation in epidemic control”, Sci Rep 9, 3505 (2019)

All authors carried out the research and edited the manuscript. L.-S.Y. and T.P. wrote the main body of the manuscript. T.P. prepared Figures 1 and 2 and S.R. prepared Figure 3. 4. Chap. 5: S. Ruschel and S. Yanchuk, “Chaotic bursting in semiconductor lasers”, Chaos 27,

114313 (2017)

SR performed the rigorous analysis and investigation including the implementation of all computations; generated all figures; wrote, proofread and edited the text. SY designed the study, advised on the investigation, implementation and figures, proofread and edited the text. SY contributed to the numerical computation and provided the program to compute Lyapunov Exponents.

5. Chap. 6: S. Ruschel and S. Yanchuk, “Delay-Induced Switched States in a slow-fast system”,

Phil. Trans. R. Soc. A. 377, 20180118 (2019)

SR performed the rigorous analysis and investigation including the implementation of all computations; generated all figures; wrote, proofread and edited the text. SY designed the study, advised on the investigation, proofread and edited the text.

(8)

(9)

Acknowledgment

I would like to express my sincere appreciation toward my supervisor Serhiy Yanchuk for taking me on as a doctoral student and his guidance during my doctoral studies in Delay Differential Equa-tions. He never failed to be encouraging, yet critical at the same time. I am thankful for the collabo-ration with Tiago Pereira at the University of São Paulo, and Lai-Sang Young at New York Univer-sity. I am indebted to Serhiy, Tiago and Lai-Sang for their continuous support and patience with me. I also thank Matthias Wolfrum, Jan Sieber and the Applied Dynamical Systems group at TU Berlin for many inspiring discussions. I profited greatly from their expertise and individual viewpoints on dynamical systems and nonlinear dynamics. I am thankful to the German Research Foundation for the financial support within the International Research Training Group “Dynamical Phenomena in Complex Networks: Fundamentals and Applications” IRTG 1740/ TRP 2015/50122-0 at Hum-boldt University Berlin and the University of São Paulo, co-funded by FAPESP; as well as CRC 910 “Control of self-organizing nonlinear systems: Theoretical methods and concepts of applica-tion” at Technical University Berlin. I enjoyed the hospitality of the University of São Paulo in São Carlos for a six months stay as well as of the New York University, Imperial College London and University of Exeter for multiple short visits.

I want to express my deep gratitude toward my family and my wife, for their unfailing support studying Mathematics and in life.

(10)

(11)

1. Introduction

It is the aim of this thesis to provide a nonlinear dynamics perspective for the multiple-time scale analysis of Delay Differential Equations (DDEs) and give a motivated introduction to what multiple

time-scales refers to in this context. Many nonlinear dynamical systems exhibit oscillations on

different time-scales that are characterized by interchanged phases of slow and fast motion. Their analysis relies on the study of reduced systems, the properties of which are then recombined using geometric singular perturbation theory [Fen79, Jon06, Kue15]. Analogous tools are increasingly available for Delay and Functional Differential Equations [DvGVW95, BLZ98, BLZ00, HMO02], and despite several advances [CMP83, IA89, AS01, FM02, Chi03, CCLG05, CSE09, GRMP10, WEH+13, KT16, TTC+15], the specific effects of time-delay on multiple time-scale systems are

largely uncharted territory.

This is a compelling subject in nonlinear dynamics, as systems with time delay are on the fore-front of modern methods and advances in dynamical systems, such as chaos [Far82, adHW83, HDY+_{11, LPM15] and its control [Pyr92, Pyr96, Pyr98, SS08]. At the same time, DDEs are highly}

relevant in various fields of modern applications including secure communication [ASL+_05],

in-formation processing [ASV+_{11], machine learning [LSB}+_{12], and across many fields such as}

mathematics, physics, biology, engineering, economics, and many others [Ern09, Ata10, Smi11, EJWY17]. Usually the delays are caused by finite-time communication speeds between constituent units of these systems, such that the corresponding model equations depend on past states. Consider for example the DDE

x1(t) = f (x(t), x(t_{´ τ)),} (1.1) where x(t) P Rd_{, t, τ} _{P R, f P C}1_(Rd_ˆRd_{, R}d). Equation (1.1) is similar to an Ordinary

Differen-tial Equation (ODE) only that the right hand side of Eq. (1.1) that determines the value of the deriva-tive at time t, depends on the value of x at moments τ time ago, as well as on the present value. It can be considered as the simplest case with only one discrete delay, or rather depending on only one de-layed variable. The τ-notation for the delay is commonly used in the field and is adopted throughout the manuscript. The text requires the reader to have a basic understanding of the geometric theory and nonlinear dynamics of Ordinary Differential Equation (ODE) [Chi99, SSTC01, GH04, Kuz13], yet provides all necessary concepts for Delay and Functional Differential Equations (such as notion of solution, semi-flow, etc.) in a brief preliminary section 2.1, that is based on classic monographs in the field [BC63, Hal66, Dri77, HV93, DvGVW95, HMO02]. Additionally, the relevant con-cepts of the available invariant manifolds theory for semiflows [HV93, BLZ98, BLZ00, HMO02] are discussed in Sec. 2.2.

(14)

1.1. Delays introduce additional time-scales

This section introduces the basic mechanisms, how delay introduces additional time-scales. To start off, consider the case of a single small delay.

1.1.1. Small delay

Assume for now that the delay is small τ ! 1. In this case, the dynamics of Eq. (1.1) is es-sentially low dimensional and can be related to a system of ODEs with a multiple time-scale structure [Chi03, Chi04]. The idea is to formally Taylor-expand the functions τ ÞÑ x(. ´ τ) and τ _{ÞÑ f(x(.), x(. ´ τ)) as series in powers of τ. For each fixed t, one obtains}

x(t´ τ) = x(t) ´ τx1_{(t) +} τ2

2x

2_{(t) + O(τ}3_),

and intuitively, if the derivatives are bounded uniformly as τ Ñ 0, f is sufficiently smooth and τ sufficiently small, one can approximate the right hand side of Eq. (1.1) using

f(x(t), x(t_{´ τ)) =f(x(t), x(t)) ´ τB}2f(x(t), x(t))x1(t) + . . . (1.2) +τ 2 2 ( B2f(x(t), x(t))x2(t) +B22f(x(t), x(t))(x1(t) )2) + O(τ3),

where B2f(x(t), x(t)) is the Jacobian of f with respect to its second argument evaluated at x(t).

This asymptotic expansion is referred to as post-Newtonian approximation [Chi03, Chi04]. The term is physically motivated; finite propagation speeds, e.g. the speed of light, between constituent units of a physical system may lead to delays, which are not considered in Newtonian/Classical Physics.

Instead of investigating Eq. (1.1) directly, one can now study the resulting system of ODEs when terms of higher order are neglected. Well known examples of this technique include the study of high-speed machining [CSE09] and gravitational radiation damping [CKMR01]. Indeed, introduc-ing w(t) = x1_(t)_{, ε = τ}2_/2_{, and assuming B}

2f(x(t), x(t))is invertible, Eq. (1.1) can be formally

rewritten as εw1 =_B2f(x, x)´1 ( w_{´ f(x, x) +}?2ε_B2f(x, x)w´ εB22f(x, x)w2 ) , (1.3) x1 = w, (1.4)

where the dependence on time has been omitted for brevity. Note that if the matrix B2f(x, x)is not

invertible, Eq. (1.1) can be viewed as a mixed system of Ordinary and Delay Differential Equa-tions and can be reformulated respectively. One should now pay special attention to the multiple

(15)

1.1. DELAYS INTRODUCE ADDITIONAL TIME-SCALES time-scale structure of Eqs. (1.3)–(1.4); while x1_(t)stays bounded as ε Ñ 0, outside an ε-small

neighborhood of the critical manifold C0 =

!

(w, x)_{P R}2d_{| w = f(x, x)}),

w1(t)grows unbounded as ε Ñ 0. Generally, x is referred to as slow variable and w is called fast

variable respectively. Formally setting ε = 0, yields the Differential-Algebraic Equation

0 = w_{´ f(x, x),} (1.5)

x1 = w, (1.6)

and the flow defined by (1.5)–(1.6), restricted to the critical manifold coincides with the flow of Eq. (1.1) for τ = 0. Complementarity, rescaling time in Eqs. (1.3)–(1.4) as t ÞÑ ˜t= t/ε and letting

˜

w(˜t) = w(t), ˜x(˜t) = x(t), one obtains the “desingularized” system ˙˜ w=_B2f(˜x,x)˜ ´1 ( ˜ w_{´ f(˜x, ˜x) +}?2ε_B2f(˜x,x) ˜˜ w´ εB22f(˜x,x) ˜˜ w2 ) , (1.7) ˙˜x = ε ˜w, (1.8)

where ˙ = d/d˜t is the derivative with respect to the fast time-scale ˜t (in this case). Setting ε = 0 here, one recovers C0as the set of equilibria of

˙˜

w=_B2f(˜x,x)˜ ´1( ˜w´ f(˜x, ˜x)) , (1.9)

with respect to the fixed parameter ˜x. The behavior of solutions in the neighborhood of the in-dividual equilibrium ˜w(˜x) = f (˜x,x)˜ is determined by the eigenvalues of the matrix B2f(˜x,x)˜ ´1.

One may now ask, whether C0 persists for ε ą 0; as well as what is the flow on this object,

and whether it attracts nearby solutions. These questions can be addressed using geometric sin-gular perturbation theory [Fen79, O’M91, Jon06, Kue15]. The following vocabulary is needed:

˜

w(˜x)is called hyperbolic, if B2f(˜x,x)˜ ´1 does not have eigenvalues with real part zero, and C0 is

called normally hyperbolic if B2f(˜x,x)˜ ´1 is hyperbolic for all ˜x with ( ˜w(˜x), ˜x)P M0. Indeed, for

sufficiently small ε and normally hyperbolic C0, there exists a locally invariant, normally

hyper-bolic manifold Cε (ε-close and diffeomorphic to C0) for ε ą 0; and the dynamics of Eqs. (1.3)–

(1.4) can be projected onto Cε in a similar way as in the limiting case (1.5)–(1.6). The rigorous

result is given in Sec. 2.2, after the needed mathematical concepts have been established. The most important issue should now be, whether and in which sense solutions of Eqs. (1.3)–(1.4) relate to solutions of the original DDE (1.1). In fact, one can identify a low-dimensional equivalent of Cε in Eq. (1.1), the so-called inertial manifold, such that the vectorfield on Cε coincides with the

(16)

order of approximation (1.2), see [Chi03, Chi04] for details. In this sense, a solution to (1.3)–(1.4) corresponds to a solution of Eq. (1.1). This is certainly only a glimpse into geometric singular perturbation theory and the mathematical framework for DDEs, yet shows that studying DDEs de-mands geometric theory for Ordinary as well as Delay Differential Equations, and that “adding” a small delay introduces an additional fast time-scale to the original system.

1.1.2. Large delay

On the other hand, motivated by the discovery of chaos in iterative maps x ÞÑ h(x), x P Rd_{, h}

P C1(Rd), a wide range of real-world dynamical phenomena have been modeled as linear filters

(convolutions) of such mappings

εx1(t) + x(t) = h(x(t´ 1)); (1.10) and rescaling time as t ÞÑ ˜t= t/ε, as well as letting ˜x(˜t) = x(t), reveals an equivalent DDE

˙˜x(t) + ˜x(t) = h(˜x(˜t´ ε´1_)), _(1.11)

with large delay τ = ε´1, where˙ = d/d˜tis again the derivative with respect to ˜t. Classical

exam-ples of this approach cover oscillatory, excitable as well as chaotic behavior [adHW83] in phys-iological control systems [MG77], nonlinear optics [Ike79], population dynamics [GBN80] and neuroscience [MW89]. Particularly in nonlinear optics, the finite-time communication delays are typically much larger than the device’s internal timescales [SGOMF13, Lüd11], and therefore, give rise to rich dynamical phenomena [GHKS81, FvTL+_{96, CCLG05, LPM13]. Specific examples of}

delay systems with (multiple) large delays include semiconductor lasers with two optical feedback loops of different length [YG14, YPW+_{15, MJB}+_{15], and ring-cavity lasers with optical feedback}

[FRSS08, OLV+_{12, JKL17]. Additional examples can be found in applications to biological}

sys-tems, when a corresponding separation of time-scales is justified [CV96, SC00, SF17, RPYY18]. The motivation is clear; for small ε ą 0, one might expect the solution to Eq. (1.10) to resemble the solution to the Difference Equation

x(t) = h(x(t_{´ 1))} (1.12) with continuous argument t P [t0´ 1, 8) and same initial condition, at least for some time [IA89].

See [SMR93] for details on Difference Equations and [ASY96, You13b] for an introduction to chaos. The solution to Eq. (1.12) is determined via the iterative mapping ξ ÞÑ h(ξ) P Rd_,which is

called reduced or singular map in this context. Clearly, the fixed points ␣ξ _{P R}d

(17)

equilib-1.1. DELAYS INTRODUCE ADDITIONAL TIME-SCALES rium solutions of (1.10); the local behavior of solutions is in general different however. Specific cases have been studied in much detail, see [Hay50, CvdD86, BC94, AH80, RW03] and references therein. However recently [YW10, LWY11], the asymptotic behavior as ε Ñ 0 of spectrum of equi-librium solutions to (1.11) has been studied rigorously, also for the general equation (1.1) with large delay. Many of the underlying ideas can be generalized to the case of multiple large delays, which is rigorously studied in Chap. 3, and the results of which are used throughout the remaining Chaps. 4– 6. The spectrum in this case resembles the spectrum of spatially extended systems, and the analo-gies have been explored to a large degree [AGLM92, GP96, YG14, YLWM15, YG15, YG17].

Assume for now that d = 1. Then, one can easily show that |h1_(ξ)_{| ă 1 implies x(t) Ñ ξ in both}

Eqs. (1.10) and (1.12) for all x(t), t P [t0´ 1, t0)in a small neighborhood of ξ (See Corollary 3.11

and Example 3.14); and in this case the solution to Eq. (1.10) to resembles the solution to Eq. (1.12). Suppose otherwise, and assume that singular map posses a stable period-2 orbit ξ. It should be noted that, similar to the case of small delay, Eq. (1.10) holds a multiple time-scale structure. If for some fixed t, (x(t), x(t ´ ε´1₎₎_{P R}2_{lies outside an ε-small neighborhood of the singular manifold (the}

solution manifold of the singular map) S0 =

!

(x, y)_{P R}2d_{| x = h(y)}),

|x1_(t)_{| grows unbounded as ε Ñ 0. Compared to Subsec. 1.1.1, x(t) plays the role of a fast variable.}

In order to build some intuition, rewrite Eq. (1.10) as

x(t) =h(x(t_{´ 1)) +}[x(t0)´ h(x(t0´ 1))] e´(t´t0)/ε (1.13)

+ żt

t0

e´(t´s)/εh1(x(s´ 1))x1_(s_{´ 1)ds}

using the Variation of Constants formula. As long as the derivative of x is bounded independent of ε, and ε sufficiently small, this is a small perturbation of (1.12). Only when the derivative is large with respect to 1/ε this viewpoint becomes inadequate. The integral representation (1.13) allows for the following insight: As long as |x1_(t)_{| is bounded and independent of ε, i.e. |x}1_(t)_{| ! C/ε}

for some C ă 8, it immediately follows that x(t) = h(x(t ´ 1)) + O(ε) for all t. On the other hand, if the derivative is large, yet only on a short interval (a, b) Ă [t0´ 1, t0],with |b ´ a| Ñ 0 as

ε_{Ñ 0, and c/ε ă |x}1(θ)_{| ă C/ε for some c ă C; then x(t) = h(x(t ´ 1)) + O (1) for θ P (a, b).} So the derivative may act as a large perturbation within (a, b). However then, for θ P (b, 1] one can easily show that this perturbation decays fast as x(t) = h(x(t ´ 1)) + O(e´(t´b)/ε)+ O (ε). This

results in the solution segment on [t0, t0 + 1]to appear as closely related to the solution segment

on [t0´ 1, t0], but possibly shifted with respect to that segment. Indeed, under certain conditions

on h, the period-2 orbit of singular map persist as periodic solution for sufficiently small ε, where the “plateaus” given by the period-2 solution of the singular map are connected by timely fast

(18)

transition layers on time-scale ε. These so called square wave solutions have been extensively studied using matched asymptotic expansions [FM02, Fow05], geometric methods for Delay and Functional Differential Equations [HL86, CMP83, CLMP89, Niz04]; and their stability properties has received a lot of attention [HH94, MPS96a]. Under certain conditions on h, Eq. (1.10) displays periodic solutions with period slightly larger than one [FMP89, Kri08, GRMP10], that give rise to interesting spatiotemporal phenomena such as coarsening, nucleation and spatially localized solutions in DDEs [WED+12, GMZY12, GMZY13, YG14, YG15, EJWY17, YG17, YRSW19].

Chap. 6 studies delay induced switched states, as a particular class of solutions to Eq. (1.10), which are not necessarily periodic and only characterized fast transitions layers induced by the delay.

1.2. Overview of main results

Chapter 3 provides a rigorous description of the spectrum of linear DDEs with multiple hierarchical large delays. This is an extension of the results obtained by Lichtner et al. [LWY11] for a single large delay to multiple large hierarchical delays, and under more general non-genericity conditions. It is shown that the spectrum of the linear delay differential equation x1_{(t) = A}

0x(t) + A1x(t´

τ1) + . . . + Anx(t´ τn)with multiple hierarchical large delays 1 ! τ1 ! τ2 ! . . . ! τn splits

into two distinct parts, namely the strong spectrum and the pseudo-continuous spectrum. As the delays tend to infinity, the strong spectrum converges to specific eigenvalues of A0, the asymptotic

strong spectrum. Eigenvalues in the pseudo-continuous spectrum however, converge to the imag-inary axis. It is shown that after rescaling, the pseudo-continuous spectrum exhibits a hierarchical structure corresponding to the time-scales τ1, τ2, . . . , τn. Each level of this hierarchy is

approxi-mated by spectral manifolds that can be easily computed. The set of spectral manifolds comprises the so-called asymptotic continuous spectrum. The position of the asymptotic strong spectrum and asymptotic continuous spectrum with respect to the imaginary axis completely determines stabil-ity. In particular, a generic destabilization is mediated by the crossing of an n-dimensional spectral manifold corresponding to the timescale τn.

Chapter 4 discusses the nonlinear dynamical behavior of solutions to a susceptible-infected-isolated model (SIQ) with two delays. The derivation of the model is provided in Appendix A. For the parameter ranges considered, this system is influenced by one large and one small delay, as well as two 1-parameter families of equilibrium solutions. After applying the results of Chap. 3, it is shown that the state space is foliated by codimension-1 manifolds that can be computed semi-explicitely. This allows to lift the local stability information along the 1-parameter families of equilibrium solutions to the global behavior of solutions. This kind of information is usually not available, and allows for valuable insight into consequences of delays and imperfect implementa-tion of isolaimplementa-tion and the predicimplementa-tion the outcome of an incipient disease outbreak.

(19)

1.2. OVERVIEW OF MAIN RESULTS chapter focuses on the dynamic mechanisms for low frequency fluctuations in semiconductor lasers subjected to delayed optical feedback, using the Lang–Kobayashi model. This system of DDEs displays pronounced envelope dynamics, ranging from erratic, so called low frequency fluctuations to regular pulse packages, if the time scales of fast oscillations and envelope dynamics are well separated. The additional timescale introduces envelope-dynamics along the solution of the linear fast system that, using Chap. 3, exhibits a weak-strong spectral splitting. Using invariant manifolds theory, one can establish low frequency fluctuations as a weakly chaotic slow motion along the the center stable and center unstable manifold of the laser off-state. As additional evidence, the Lyapunov exponents of the system is computed and one shows the corresponding weak-strong splitting of the Lyapunov spectrum. It can be shown that the onset of low frequency fluctuations coincides with events of finite time strong chaos in the system.

Chapter 5 discusses two-component delay system as a minimal model for delay-induced switched states. Previously, such systems have been reported to model switching in optoelectronic experi-ments [WED+_{12, WEH}+_{13, LPM13, LPM15, EWBH16], where each switching induces another}

one after approximately one delay time. It is shown that these solutions are fundamentally different from the typical relaxation oscillations expected from singular geometric perturbation theory for Ordinary Differential Equations [Fen79, Jon06, DGK+12, Kue15]. Necessary conditions for their

existence are derived, and the case of weak-instability of the corresponding family of equilibrium solutions is studied in detail.

(20)

(21)

2. Preliminaries

For completeness, this chapter discusses results on existence, uniqueness, and smooth dependence on initial conditions of the initial value problem

x1(t) = f (x(t), x(t_{´ τ}1), . . . , x(t´ τn)), x(t) = φ(t) for all t0´ τnď t ď t0, (2.1)

where x(t) P Fd, with Fd _{= R}d or Fd _{= C}d, n, d P N, f : Fdˆn _{Ñ F}d, φ : Fd

Ñ Fd, t 0 P R

and without loss of generality 0 ă τ1 ď τ2 ď . . . ď τn. Notice that in order to uniquely define

the future of the variable x(t), one must prescribe a function ψ on the interval [t ´ τ, t], which therefore can be thought as the state of the system. The function ψ is sometimes called the history or history function. Thereafter, a brief overview of the available geometric theory for DDEs is provided. These known results will be used throughout the manuscript.

2.1. Existence, uniqueness, and linearized stability

This section follows the classic textbook of Hale and Verduyn Lunel [HV93] in spirit and in nota-tion. Denote Ck_{:= C}k_([´τ

n,0] , Fd

)

the space k-times continuously differentiable functions and }.}kthe corresponding norm, such that for all φ P C

k, }φ}k= sup ´rďθă0 › ›φ (θ)›› +ÿk j=1 sup ´rďθă0 › › ›φ(j)(θ)››_{› ă 8,} where }.} is a norm on Fd. (Ck_,_}.} k )

is a Banach space. Let α ě 0. If x P C0_([t

0´ τn, t0+ α] , Fd

) for any t P [t0, t0+ α], let xtP C0be defined by xt(θ) = x (t + θ), for all ´τn ď θ ď 0. Hence, for

a fixed t P [t0, t0+ α] ,the segment xtof x is an element of the space (Ck,}.}_k

)

of Ck-continuous

functions from [´τn,0]into Rn. If Ω Ă R ˆ C0, then a function F : Ω Ñ R defines an RFDE(F),

referring to Retarded Functional Differential Equation, by

˙x (t) = F (t, xt) (2.2)

with initial conditions (t0, φ), such that xt0 = φ.

Theorem 2.1 (Hale and Verduyn Lunel, [HV93] p. 49,51). If F _{P C}p_{(Ω, F}d_{), p ě 1, and}

(t0, φ) P Ω, then the solution x (t0, φ, F) (t) of the RFDE(F ) through (t0, φ) is unique and

p-times continuously differentiable with respect to (φ, F ) for t in any compact set in the domain of definition of x (t0, φ, F) (t). Furthermore, for each tě t0, the derivative of x with

(22)

respect to φ, Dφx(t0, φ, F) (t) is a linear operator from C0 to Fd, Dφx(t0, φ, F) (t0) = I, the

identity, and Dφx(t0, φ, F) ψ (t) for each ψ P C0 satisfies the linear variational equation

˙y (t) = DφF(t, xt(t0, φ, F)) yt. (2.3)

The proof is based on the Banach fixed-point theorem, similar to the case of an Ordinary Dif-ferential Equation. To apply Theorem 2.1 to the specific setting of DDEs, notice that one recovers Eq. (2.1), if F can be written as F (t, xt) = f(t, ev (τ1, xt) , . . . , ev (τn, xt)

)

, where ev (θ, φ) = φ(θ)is the evaluation (operator). Note that the evaluation operator is Cp, if the function to be

eval-uated is. It is important to note, that the solution starts (only) continuous and its regularity improves with every delay interval. In general, let f be Cp, p ě 1. If φ P Cp, then the solution x (t

0, φ, f) (t)

of (2.4) is p-times continuously differentiable in (φ, f). Moreover, if φ P C0_{, then the solution}

x(t0, φ, f) (t)of (2.4) is p-times continuously differentiable in (φ, f) for all t ą t0+(p´1)τn. This

is referred to as the smoothing property. Thus for the system at hand, choose C = C0_([´τ

n,0] , Fd

) as space of initial conditions. The solution defines a C1 semi-flow Tt _{: φ}

ÞÑ xt(φ)on C [HV93].

Clearly, T0 ₌ id and Tt+s _{= T}t

˝ Ts. For a given initial state φ the orbit through φ is the subset

␣ Tt_φ(

tě0of the state space, and one refers to the dynamics of Eq. (2.1) as the collection of orbits

filling the state space.

Recall that the ω-limit set ω(φ) of φ P C under the semi-flow Ttis defined to be

ω(φ) =tψ P C | Tt_φ

Ñ ψ for some sequence tnÑ 8u.

If xt(φ)exists for all t, xt(φ)is C1 with a uniform bound on its derivatives for all t ě t0+ τn. By

the Arzela-Ascoli Theorem then, ω(φ) is nonempty and compact in C (with its C0 norm). On the

basis of Theorem 2.1, one relies on the principle of linearized stability to determine the asymptotic behavior of small perturbations to ω(φ).

In the simplest case ω(φ) may consist of an equilibrium solution, that is a constant function ψ _{P C, such that ψ(θ) = 0 for all θ P [´τ}n,0], without loss of generality. Throughout the paper,

the shorthand ˆ. is used, to refer to constant equilibrium solutions ˆ0 P C, such that ˆ0(θ) = 0 for all θ _{P [´τ}n,0]and f(0, 0, . . . , 0) = 0. For Eq. (2.1), the corresponding linearized system (2.3) takes

the form

ξ1(t) = A0ξ(t) + A1ξ(t´ τ1) +¨ ¨ ¨ + Anξ(t´ τn), (2.4)

where Ai is the Jacobian with respect to the i-th argument of f, i.e. Ai = Bif(0, 0, . . . , 0),

0 ď i ď n. This can be easily seen inserting the ansatz x(t) = 0 + ξ(t) into Eq. (2.1), and neglecting terms of higher order after Taylor-expansion with respect to ξ. Equation (2.4) can be thought of as similar in spirit to an Ordinary Differential Equation (ODE) except that it may exhibit so-called small solutions; those are solutions that “collide” with the trivial solution in finite time,

(23)

2.1. EXISTENCE, UNIQUENESS, AND LINEARIZED STABILITY say t1, and equal zero for all t ě t1. Apart from this peculiarity, that is up to small solutions,

any solution of (2.4) can be written as a superposition of exponential functions as in the case of ODEs [HV93]. In particular, the long term behavior of the solution as t Ñ 8 is governed by the characteristic exponents.

In this sense, solving (2.4) is equivalent to finding nontrivial solutions to the matrix-valued quasi-polynomial equation ∆(λ)v = 0, where ∆ : C Ñ Cdˆd_,

∆(λ) :=_{´λI + A}0+ n

ÿ

k=1

Ake´λτk (2.5)

is the characteristic matrix. A nontrivial solution v exists, if and only if there is λ P C such that ker ∆(λ) ‰ H, or equivalently det ∆(λ) = 0. For simplicity, assume λ is a simple root of det ∆(λ). Together with a corresponding 0 ‰ v P ker ∆(λ) Ď Cd, it gives rise to a solution t ÞÑ v exp(λt) of

Eq. (2.4). See Ref. [HV93] for further details. The pair (λ, v) P C ˆ Cd is called an

eigenvalue-eigenvector pair and the entirety of eigenvalues λ is called the spectrum Σ := ␣λ_{P C | det ∆(λ) = 0}(

of (2.4). Hence, the problem consists of describing the location of complex-valued solutions to the characteristic equation

χ(λ) := det ∆(λ) = 0. (2.6) A complete description of the spectrum can be formidable task even for a single delay, and is gener-ally unfeasible for two or more. From the point of view of applications, one is genergener-ally interested in the case when ℜ(λ) ă 0 for all λ P Σ and the equilibrium is exponentially stable. Specific cases therefore have been studied in much detail, see [Sch11, Hay50, CvdD86, BC94, AH80, RW03] and references therein. For each fixed ε ą 0, much is known about the solutions of (2.6). Firstly, there are countably many solutions that continuously depend on parameters. Secondly, solutions accu-mulate at minus infinity and within each vertical stripe [α, β]ˆiR Ď C there are only finitely many solutions [BC63, HV93]. In particular, β can be chosen finite [HV93]. The following theorem was proven for a single delay, i.e. n = 1.

Theorem 2.2 (Small Delays are harmless, Smith [Smi11, p. 48]). Let λ1, λ2, ..., λd be the

distinct eigenvalues of A0 + A1, let δ ą 0, and let s P R satisfy s ă miniℜ(λi), 1 ď i ď d.

Then, there exists r0 ą 0 such that if 0 ă r ă r0 and χ(λ) = 0 for some λ then either

ℜ_(λ)_{ă s or |λ ´ λ}_i_{| for some i.}

This is a well known and sometimes misunderstood concept in DDE nonlinear dynamics. For sufficiently small delay, Σ can be considered as a small perturbation of the set tλ1, λ2, ..., λdu and

(24)

de-lays are however not harmless in the sense, that passing to negative delay produces an Advanced Differential Equation, which is not even generate a semi-flow (in forward time). In other words, small delays are not a generic small perturbations to the corresponding equation with delay zero. As one of the main contributions, Chap. 3 provides a rigorous description of the spectrum Σ, when the involved time delays are hierarchically large, that is 1 ! τ1 ! τ2 ! . . . ! τn.As in this case, a

complete description is possible as well, one might say that large delays are equally harmless. The case when they are not, corresponds to the so-called weak instability, see Chap. 3 for the rigorous analysis.

2.2. Geometric theory and invariant manifolds

The eigenvectors v = vλ corresponding to λ P Σ with negative, zero, and positive real part define

the (strong) stable eigenspace Es

Ă C, the center eigenspace Ec

Ă C, and (strong) unstable eigenspace Eu

Ă C of the trivial equilibrium of Eq. (2.4), respectively. They give rise to the corresponding splitting of state space C = Es

‘ Ec

‘ Eu_.One should remark here that, for simple

geometric reasons, the equilibrium is called hyperbolic if Ec ₌

H. The following theorem states that this splitting persists in the nonlinear system Eq. (2.1). Recall that ˆ0 P C is the function such that ˆ0(θ) = 0 for all θ P [´τn,0]. Fix a given neighborhood V Ă C of ˆ0 in Eq. (2.1), and define

the (local) strongly stable manifold Wssof ˆ0 of Eq. (2.1), as the nonlinear equivalent of Es, i.e. the

set of points φ P C (φ P V ) such that the solution xt(φ)P C approaches ˆ0 exponentially as t Ñ 8.

The (local) strongly unstable manifold Wsu of ˆ0 can be defined analogously. Similarly, define

Wc_{X V =} ␣ψ _{P C|ψ = φ + h(φ), φ P E}c_{X V}(, where h : Ec _{Ñ E}s_{‘ E}u is a C1_{-mapping with h(ˆ0) = ˆ0 and D}

φh(ˆ0) = ˆ0,as the (local) center

manifold of ˆ0; the nonlinear equivalent of Ecin Eq. (2.1).

Theorem 2.3 (Hale and Verduyn Lunel, [HV93] p. 314.). If F is a Ck_{-function, then there}

is a neighborhood V of ˆ0 _{P C, such that each of the sets W}ss_{, W}su _{and W}c _{exists and is a}

Ck-manifold of C. The manifolds Wss_{, W}su_, _{are uniquely defined, whereas the manifold W}c

is not. Furthermore, every invariant set of Eq. (2.1) that remains in V must belong to Wc_.

In particular, one has a (local) splitting V = Wss

‘ Wc

‘ Wsu of the neighborhood of the

trivial equilibrium in Eq. (2.1). Wss_{, W}c_{, W}suare positively invariant under the semi-flow Tt_,i.e.

Tt_Wss_{Ă W}ss, Tt_Wc _{Ă W}c_,and Tt_Wsu _{Ă W}sulocally for all t ě 0. By slight abuse of notation,

one says Wss_{, W}c_{, W}su are invariant for short. If Ec _{= E}u ₌_{H it follows that W}su _{= W}c ₌_H,

and the equilibrium is asymptotically stable in Eq. (2.1). A rigorous proof of Theorem 2.3 involves the extension of the concept of solution to the space Fd_{ˆ L}8_([_´τ

(25)

2.2. GEOMETRIC THEORY AND INVARIANT MANIFOLDS general history functions, contained only in larger spaces such as measurable, or bounded functions on the interval [´τn,0]would uniquely define solutions Eq. (2.1) to provided x(0) is defined as

well; it turns out however, that the space of functions above is the correct space to define an inner product with respect to functions of normed bounded variation that is needed to project to the corresponding eigenspaces [DvGVW95]. Suppose that that Eq. (2.1) possesses a continuous 1-parameter family M of equilibrium solutions. Then, M is positively invariant with respect to the semiflow Tt(invariant for short), and by the Arzela-Ascoli Theorem, M is a Ckcompact connected

manifold. As an immediate consequence of Theorem 2.3, the splitting Wss

m ‘ Wmc ‘ Wmsu varies

continuously along M, and M itself belongs to the center manifold Wc of at least one m P M. As

a result, one can study the local neighborhood of M as foliated by the (strong) stable and (strong) unstable manifolds Wss

m, Wmsu, see Chap. 4 for an explicit example. One may now ask whether

under small perturbations, this splitting persist, as well as whether M itself persist as an invariant manifold [HMO02]. The rigorous answers to these questions for a general C1_{-semiflow are due}

to [BLZ98, BLZ00]. The most important results are collected below. Let ε ą 0, and M Ă C be a C2 _{compact connected positively invariant manifold (it is argued in [BLZ00] that C}1 _{is also}

sufficient).

Definition 2.4. M is said to be normally hyperbolic if the local neighborhood along M (the so-called vector bundle) splits continuously into

Xu_{‘ X}c_{‘ X}s where Xc _{is tangent to M , such that}

1. This splitting is invariant under the linearized semiflow DTt_{, such that for m}

P M, t_{ě 0, and m}1 = Tt(m), DTt(m)ˇˇ_Xα m : X α m Ñ X α m1

holds for α = u, c, s and DTt_{(m) is an isomorphism from X}α

m onto Xmα1.

2. DTtˇˇ

Xu expands and does so to a greater degree, than does DT

tˇˇ

Xc while DT

tˇˇ Xs

contracts and does so to a greater degree, than does DTtˇˇ

Xc. More specifically, there

exists t1 ě 0, and λ ă 1 such that for all t ě t1 and m P M,

λ_{¨ inf}!››DTt(m) xu›_{› : x}u P Xu m,}x u } = 1)ą max " 1,››_›DTt_(m)ˇ_ˇ Xc m › › › * , λ_{¨ min} " 1,inf!››DTt_{(m) x}c›› : xc P Xc m,}x c } = 1)* ą›››DTt(m)ˇˇ_Xs m › › › .

(26)

as long as it is normally hyperbolic, that is as long as dim Ec _{= 1.}Therefore, the following results

apply immediately to Chaps. 5–6. Denote

Xα(ε) =␣(m, xα) : m_{P M, x}α _{P X}mα,}x α

} ă ε( for α = u, c, s, such that Θ (Xu_(ε)

‘ Xs_(ε))

is an small tubular neighborhood of M, where Θ is defined as Θ : Xu

‘ Xs

Ñ X, Θ (m, xu_{+ x}s_{) = m + x}u _{+ x}s. There is a more technical

smoothness assumption on Xα, which the authors argue can be removed. For more details refer

to [BLZ00, BLZ98]. Then following Theorem is analogous to the Fenichel’s Theorem for ODEs [Fen79].

Theorem 2.5 (Bates, Lu, Zeng [BLZ98]). Suppose that Tt_{is a semiflow on C and M is a C}2

-compact connected invariant manifold on which Tt _{is normally hyperbolic for t sufficiently}

large. Suppose also that for each m_{P M, DT}t

|Xu

m is an isomorphism. Let t1 be large enough

and be fixed and N be a fixed neighborhood of M. For any ε_{ą 0 there exits σ ą 0 such that} if ˜T is a C1_{-map which satisfies}_{|| ˜}_T _{´ T}t1_||

C1_{(N )} ă σ, then

1. Persistence: ˜T has a unique C1 _{compact connected invariant normally hyperbolic}

invariant manifold ˜M near M .

2. Convergence: M converges to M in the C1_{-topology as} _{|| ˜}_T _{´ T}t1_{|| tends to 0;}

3. Existence: ˜T has unique C1_{-invariant manifolds W}cs_{(ε) and W}cu_{(ε) in a tubular}

neigh-borhood N (ε) of M, which at M are tangent to the center-stable vector bundle ˜Xc_{‘ ˜}_Xs

and the center-unstable vector bundle ˜Xc_{‘ ˜}_Xu_,_{respectively.}

4. Characterization: ˜ Wcs(ε) = !x0 P N(ε)| ˜Tk(x0)P N(ε), for k ě 1, ˜Tk(x0)Ñ ˜M as kÑ 8 ) , ˜ Wcu(ε) =␣x0 P N(ε)|Dtxkuką0 Ă N(ε) , satisfying ˜Tk(xk) = xk´1 for k ě 1, xk Ñ ˜M as kÑ 8 ) , and ˜M = ˜Wcs_(ε)_{X ˜}_Wcu_(ε)

Furthermore, if ˜Tt is a C1-semiflow on X which satisfies || ˜Tt1 _{´ T}t1_||

C1_{(N )} ă σ and || ˜Tt ´

Tt

||C1_{(N )} ă σ for 0 ď t ď t₁,then ˜M is a normally hyperbolic invariant manifold for ˜Tt with

(27)

2.2. GEOMETRIC THEORY AND INVARIANT MANIFOLDS Theorem 2.6 (Bates, Lu, Zeng [BLZ98]). Let t1 ą t0. For each small εą 0, Tt has a unique

C1 _{center stable manifold W}cs_{(ε) and a center unstable manifold W}cu_{(ε) in the tubular}

neighborhood Θ(Xu_(ε)_{‘ X}s_{(ε)). These have the following properties:}

1. M = Wcu_(ε)

X Wcs_{(ε). For each m}

P M, the tangent spaces at m are 2. Tt_Wcs_(ε) XΘ(Xu_(ε) ‘ Xs_(ε)_{) Ă W}cs_{(ε). T}t_Wcs_{(ε) converges to M as t} Ñ +8, and Wcs(ε) = !x_{P Θ}(Xu(ε)_{‘ X}s_{(ε)) : T}kt1_x_{P Θ}(Xu_(ε)_{‘ X}s(ε)) , for all k ą 0) 3. Tt1(Wcs_(ε)) Ă Wcs_(ε) 4. Tt1 _{: W}cu_(ε)_X(Tt1)´1(Wcu_(ε)) Ñ Wcu_{(ε) is a diffeomorphism. If we define T}´t _on

Wcu_{(ε) in this way, then} _(Tt1)´1(Wcu_(ε)) converges to M as t Ñ 8 and

Wcu(ε) =!x_{P Θ}(Xu (ε)_{‘ X}s(ε)) : for all k ą 0, Dyk P Θ(Xu(ε)‘ Xs(ε) ) s.th. Tkt1_y k= x ) . Moreover, Wcs

ε and Wεcuare invariantly fibered as states in the following two theorems.

Theorem 2.7 (Bates, Lu, Zeng. Thm. 2.2. [BLZ00]). For small ε there exists a unique family of C1_{-submanifolds}␣_Wss

m (ε) : mP M

(

of Wcs_{(ε) satisfying:}

1. For each m_{P M, M X W}ss

m (ε) =tmu, the tangent space TmWmss(ε) = Xms and Wmss(ε)

varies continuously with respect to m on M . 2. If m1, m2 P M, m1 ‰ m2, then Wmss1(ε) = W ss m2(ε) = H and W cs_{(ε) =} Ť mPMWmss(ε). 3. For m_{P M, T}t1(Wss m (ε)) Ă WTsst1_(m)(ε) .

4. For all m_{P M and t ą 0, T}t_(Wss

m (ε)) X Θ (Xu(ε)‘ Xs(ε)) Ă WTsst_(m)(ε). 5. For x_{P W}ss m (ε) and m‰ m1 P M, we have | Tt_(x)´Tt_(m)_| |Tt_(x)´Tt_(m 1)| Ñ 0 exponentially as t Ñ +8. 6. For x, y _{P W}ss m (ε) , ˇ ˇTt_(x)_{´ T}t_(y)ˇˇ Ñ 0 exponentially as t Ñ +8.

For small ε there exists a unique family of C1_{-submanifolds} ␣_Wuu

m (ε) : mP M

(

of Wcu_(ε)

satisfying:

1. For each m_{P M, M X W}uu

m (ε) =tmu, TmWmuu(ε) = Xmu and Wmuu(ε) varies

(28)

2. If m1, m2 P M, m1 ‰ m2, then Wmuu1(ε) = W uu m2(ε) = H and W cu_{(ε) =} Ť mPMWmuu(ε). 3. For all m P M, Tt1 _{: W}uu m (ε)X(Tt1 )´1( Wuu Tt1(m)(ε) ) Ñ Wuu Tt1(m)(ε) is a diffeomor-phism. 4. For x _{P W}uu

m (ε) , if Tt(x) P Θ(Xu(ε)‘ Xs(ε)) for all 0 ă t ă t2, for some t2, then

Tt_(x)_{P W}uu Tt_(m)(ε) , for all 0ă t ă t2. 5. For x _{P W}uu m (ε) and m ‰ m1 P M, we have | T´t_(x)´T´t_(m)| |T´t_(x)´T´t_(m₁₎| Ñ 0 exponentially as t _{Ñ +8.} 6. For x, y _{P W}uu m (ε) , we have ˇ ˇT´t_(x)_{´ T}´t_(y)ˇˇ Ñ 0 exponentially as t Ñ +8. Each submanifold Wss m (ε)(Wmss(ε) )

is called a stable (unstable) fiber. The stable (unstable) foliation is the decomposition of Wsu _(Wcu₎into the disjoint union of stable (unstable) fibers.

(29)

(30)

(31)

3. The spectrum of linear DDEs with

multiple hierarchical large delays

This chapter provides a rigorous description of the spectrum of linear DDEs x1(t) = A0x(t) +

n

ÿ

k=1

Akx(t´ τk), (3.1)

with hierarchical large delays τk = σkε´k, σk ą 0, 1 ď k ď n, by studying the asymptotic

behavior of solutions to the corresponding characteristic equation χε(λ) :=det ( ´λI + A0+ n ÿ k=1 Ake´λσkε ´k ) = 0 (3.2)

as ε Ñ 0. Throughout, it is assumed that x(t) P Cd is a complex-valued, Euclidean vector of

size d and Ak P Cdˆd, Ak ‰ 0 , 0 ď k ď n are given matrices independent of time and ε ą 0.

From Sec. 2.1, one recalls the definition of the characteristic equation. As a reminder, a solu-tion to Eq. (3.2) is called an eigenvalue, and the entirety of eigenvalues is called the spectrum Σε ₌␣_λ_{P C|χ}ε_{(λ) = 0}(of Eq. (3.1). Section 3.2 shows that the spectrum splits into two distinct

parts with different scaling behavior: the strong spectrum and the pseudo-continuous spectrum. As ε Ñ 0, the strong spectrum converges to specific eigenvalues of A0, the so-called

asymp-totic strong spectrum. Eigenvalues in the pseudo-continuous spectrum converge to the imaginary axis as ε Ñ 0. It is shown that after rescaling the pseudo-continuous spectrum exhibits a hierarchi-cal structure corresponding to the time-shierarchi-cales τ1, τ2, . . . , τn.In particular, the pseudo-continuous

spectrum can be represented as a union of subsets corresponding to different timescales τk. For

1 _{ď k ă n, each of these subsets can be approximated with a k-dimensional spectral manifold,} see Sec. 3.3. The manifolds corresponding to k ă n, generically lie in the positive half plane; they extend to the negative half-plane only under certain degeneracy conditions related to the rank of the matrices Ak+1, . . . , An, which are addressed in Sec. 3.1.

The presented asymptotic approach has many advantages with respect to the direct analysis of the characteristic equation. Firstly, the spectral manifolds can often be computed explicitly, or even analytically, and the eigenvalues can be found by projecting the manifolds to the complex plane. Section 3.4 contains a corresponding numerical method to compute the exact spectrum of Eq. (3.1) for sufficiently small ε ą 0 that allows to compare the exact spectrum and its approximation given by the asymptotic spectrum. Secondly, the asymptotic spectra are by construction exact at the

(32)

imaginary axis, and therefore, the stability boundaries are completely determined by the position of the asymptotic strong spectrum and asymptotic continuous spectrum with respect to the imagi-nary axis. It is shown that a generic destabilization of Eq. (3.1) takes place via the crossing of an n-dimensional spectral manifold corresponding to the timescale τn. This destabilization scenario

is reminiscent of spatially extended systems [YLWM15, YG17]. The corresponding bifurcation is low-dimensional, i.e. at the bifurcation point a finite number of eigenvalues crosses the imaginary axis. However, the primary bifurcation is followed by further eigenvalues crossing the imaginary axis as the bifurcation parameter is increased, and the distance between these bifurcations in param-eter space scales with ε as ε Ñ 0. See also Sec. 3.4 for an example of this type of destabilization as parameters are varied.

3.1. Degeneracy spectrum

This section discusses certain genericity conditions of Eq. (3.2). At first, one might assume that the matrices Ak,0 ď k ď n have full rank, i.e. det(Ak) ‰ 0, 0 ď k ď n, as zero generically is

not an eigenvalue of Ak, or in other words, µ = 0 is generically not a solution of the polynomial

det(µI ´Ak) = 0. From the point of view of applications, however, it certainly cannot be expected

that every single component is influenced by delay. This leads to matrices Akwith low rank and is

reflected in the degeneracy spectrum of Eq. (3.2) (See Theorem 3.3).

Consider for example the case when Anis not invertible, i.e. dn:=rankAnă d. Then there are

unitary matrices Un, Vnsuch that

An= Un ( 0 0 0 A(n)_n,4 ) Vn˚, (3.3) where A(n)

n,4 P Cdnˆdn is a diagonal matrix of full rank, and Vn˚ is the conjugate transpose of Vn.

Equation (3.3) is the singular value decomposition of An; it characterizes the directions along which

An acts as a low-rank perturbation. The columns u1n, u2n, . . . , udn of Un = [u1n, u2n, . . . , udn] and

vn1, vn2, . . . , vnd of Vn = [v1n, v2n, . . . , vdn] are the left and right singular vectors of An, respectively.

Without loss of generality, the singular vectors can be arranged such that the matrices Un,1 :=

[u1

n, u2n, . . . , ud´dn k] and Vn,1 := [vn1, vn2, . . . , vnd´dk] consist of the left and right singular vectors

corresponding to the singular value zero; and U˚

n,1AnVn,1 = 0corresponds to the projection onto

the cokernel (the complement of the image of An in Cd) of An. This projection allows to define

the following spectral sets.

Definition 3.1 (non-generic spectral subsets). Let dn:= rankAn ă d.

(33)

3.1. DEGENERACY SPECTRUM vectors of An corresponding to the singular value zero. Denote

J₁(n) := U_n,1˚ Vn,1,

A(n)_j,1 := U_n,1˚ AjVn,1, j = 0, . . . , n´ 1,

and the corresponding projected characteristic equation

˜ χε_n´1(λ) :=det ( ´λJ1(n)+ A (n) 0,1 + n´1_ÿ k=1 A(n)_k,1e´λσkε´k ) . The set ˜ Σε_n´1 :=␣λ_{P C | ˜}χε_n´1(λ) = 0, ℜ(λ)_{ă 0}( is called the truncated stable τn´1-spectrum.

If A(n)_n´1,1 is again not invertible, this procedure is applied iteratively.

(ii) Recursively for all 1 ď k ď n ´ 1 (starting from n ´ 1), if det A(k+1)k,1 = 0, define

˜

Uk,1, ˜Vk,1 (notice the tilde notation) containing left and right singular vectors of A(k+1)k,1

corresponding to the singular value zero. Denote J₁(k):= ˜U_k,1˚ J₁(k+1)V˜k,1,

A(k)_j,1 := ˜U_k,1˚ A(k+1)_j,1 V˜k,1, j = 0, . . . , k´ 1

and the corresponding truncated characteristic equation

˜ χε_k´1(λ) :=det  ´λJ (k) 1 + A (k) 0,1 + k´1_ÿ j=1 A(k)_j,1e´λσjε´j  , ˜ χ0(λ) := ´λJ1(1)+ A (1) 0,1.

(iii) Define k, 1 _{ď k ď n ´ 1 as the smallest index such det A}(k+1)_k,1 = 0 for all k _{ď k ď n ´ 1} and det An= 0.

(iv) For 0_{ď k ď n ´ 1 and ˜}χε

k(λ) nontrivial, define the set

˜ Σεk :=

␣

λ_{P C | ˜}χε_k(λ) = 0, ℜ(λ)ă 0( for k _{ě k, and ˜}Σε

(34)

-spectrum. If k = 1, set ˜

S₀´ :=␣λ_{P C | ˜}χ0(λ) = 0, ℜ(λ)ă 0

( ,

and ˜S₀´ :=_{H otherwise. ˜}S₀´ is called asymptotic strong stable spectrum.

(v) If k = 1 and det J₁(1) = 0, define the matrices U , V containing left and right singular vectors of J₁(1) corresponding to the singular value zero.

These sets correspond to spectral directions along which Eq. (3.1) acts as a DDE with fewer delays or an ODE. In order to guarantee that Eq. (3.1) is indeed a DDE and cannot be transformed into a system of ODEs through variable transformations, one has to demand the following non-degeneracy condition.

Condition (ND). If detAn = 0, k = 1 and det J (1)

1 = 0,then det

(

U˚_A(1)_0,1V)_{‰ 0.}

This is a rather abstract condition. In order to build some intuition, consider the following example. Example 3.2. Let d = 2, n = 1, and the matrices A0 and A1 are given by

A0 = ( a0,1 a0,2 a0,3 a0,4 ) , A1 = ( 0 1 0 0 ) .

Clearly, rankA1 = 1ă 2, k = 1 and one readily computes

U1,1 = ( 0 1 ) , V1,1 = ( 0 1 ) , J₁(1) = U_1,1˚ V1,1 = 0, as well as U ₌ ( 0 1 ) , V ₌ ( 0 1 ) , U˚_A(1)_0,1V _{= a}_0,3_.

Clearly, a0,3 = 0 violates (ND). In this case, the system is degenerate; it corresponds to an

ODE. Indeed, straightforward computation shows that the characteristic equation χε(λ) = (a0,1´ λ)(a0,1´ λ) ´ a0,3(a0,2+ e´λσ1ε

´1

) = 0 (3.4) does not depend on e´λσ1ε´1 _{if a}

0,3 = 0, and the spectrum consists of exactly two eigenvalues

Σε ₌ _ta

0,1, a0,4u for all ε ą 0. In the following sections, it will become clear that all

eigenvalues λ _{P Σ}ε_zB

r(ta0,1, a0,4u) of Eq. (3.4), for some sufficiently small r ą 0, approach

minus infinity as a0,3 Ñ 0. Here, Br(X) = Ť_xPX

␣

z _{P C| |z ´ x| ă r}( denotes the set of balls of radius r around a set X _{Ă C.}

(35)

3.1. DEGENERACY SPECTRUM On the basis of Def. 3.1, the following Theorem 3.3 provides a hierarchical approximation of

Σε_{X tλ P C|ℜ(λ) ă 0u} by spectral subsets of truncated characteristic equations ˜χε

k, when some of the matrices Ak do not

have full rank.

Theorem 3.3. Let det An = 0, k be such that k ´ 1 ď k ď n ´ 1, and (ND) be satisfied.

Further, let ε ą 0 be sufficiently small and µε P ˜Σεk (µε P ˜S0´ for k=0). Then there exits a

small neighborhood Uε_(µ

ε)Ă C of µε such that the number of eigenvalues λεP ΣεX Uε(µε)

equals the multiplicity of µε as a zero of ˜χεk.

Proof. We proof by induction starting from the highest order k = n. We assume that det An = 0, An ‰ 0 and consider a µε such that ˜χn´1ε (µε) = 0, i.e. µε P ˜Σεn´1 as specified in

Definition 3.1. We show that for a sufficiently small ε and neighborhood Uε(µε) the number

of zeros µε of ˜χε_n´1 counting multiplicity equals the number of eigenvalues λεP ΣεX Uε(µε).

Again let the matrices Un = [Un,1, Un,2] and Vn = [Vn,1, Vn,2] contain the left and right

singular vectors corresponding to the cokernel (Un,1 and Vn,1) and image (Un,2 and Vn,2) of

An, see Def. 3.1 for details. Consider z P C, |z| sufficiently small and define

fε(z) := χε(z + µε) =det ( Cε 1(z) C2ε(z) Cε 3(z) C4ε(z) ) , where C₁ε(z) =_{´(z + µ}ε)J1(n)+ A (n) 0,1 + n´1_ÿ k=1 A(n)_k,1e´(z+µε)σkε´k_, C₂ε(z) =_{´(z + µ}ε)J2(n)+ A (n) 0,2 + n´1_ÿ k=1 A(n)_k,2e´(z+µε)σkε´k_, C₃ε(z) =´(z + µε)J3(n)+ A (n) 0,3 + n´1_ÿ k=1 A(n)_k,3e´(z+µε)σkε´k_, C₄ε(z) =´(z + µε)J4(n)+ A (n) 0,4 + n´1_ÿ k=1 A(n)_k,4e´(z+µε)σkε´k _{+ A}(n) n,4e´(z+µε)σnε ´n ,

is the block structure obtained from multiplying χε_{(λ) by} det U˚

n and det Vn from left and

right with the corresponding projected matrices

A(n)_k,1 = U_n,1˚ AkVn,1, J1(n)= Un,1˚ Vn,1 P C(d´dn)ˆ(d´dn),

(36)

A(n)_k,3 = U_n,2˚ AkVn,1, J3(n)= Un,2˚ Vn,1 P Cdnˆ(d´dn),

A(n)_k,4 = U_n,2˚ AkVn,2, J4(n)= Un,2˚ Vn,2 P Cdnˆdn,

for all 0_{ď k ď n ´ 1. Using the Schur complement formula, we obtain} fε(z) =

(

e´(z+µε)σnε´n

)dn

det( ˜C₁ε(z)) det ( ˜C₄ε(z)), where the matrices ˜Cε

1(z) and C4ε(z) are given by

˜ C₁ε(z) =_{´ (z + µ}ε)J (n) 1 + A (n) 0,1 + n´1_ÿ k=1 A(n)_k,1e´(z+µε)σkε´k ´ e(z+µε)σnε´n_Cε 2(z)( ˜C4ε(z) )´1 C₃ε(z), ˜ C₄ε(z) =A(n)n,4+ e(z+µε)σnε ´n [ ´(z + λ)J4(n)+ A (n) 0,4 + n´1_ÿ k=1 A(n)_k,4e´(z+µε)σkε´k ] . Note that det ˜C₄ε(z) = det A(n)_n,4+ O (ˇ ˇ ˇe(z+µε)ε´1 ˇ ˇ ˇ ) , Choose Uε_(µ

ε) such that ˜χεn´1(z+µε)‰ 0 and ℜ(z+µε)ă 0 for all z such that z+µε P Uε(µε).

As a result, we have ( ˜_Cε 4(z) )´1 =(A(n)_n,4)´1+ O (ˇ ˇ ˇe(z+µε)ε´1 ˇ ˇ ˇ ) , det ˜C₁ε(z) = ˜χε_n´1(z + µε) + O (ˇ ˇ ˇe(z+µε)ε´n ˇ ˇ ˇ ) , and fε(z) ( e(z+µε)ε´n )dn =det(A(n)_k,4) ˜χε_n´1(z + µε) + O (ˇ ˇ ˇe(z+µε)ε´1 ˇ ˇ ˇ ) . where ˜χε

n´1(z + µε) is as in Definition 3.1(i), and by assumption ˜χε_n´1(µε) = 0. The factor

(e(z+µε)ε´n₎dn _{remains bounded as ε}_{Ñ 0. Therefore, Rouch´e’s Theorem implies that f}

ε has

the same number of zeros as ˜χε

n´1(¨ + µε) counting multiplicity. This proves the theorem

for k = n_{´ 1. If k ă n ´ 1 this procedure has to be applied again to show that elements} of ˜Σε

k can be approximated by elements of ˜Σεk´1. The induction step k ÞÑ k ´ 1 obtaining

˜ χε

k´1 : C Ñ R and ˜Σεk´1 is completely analogous for all kě k. If k = 1, we have to grantee

that after the induction step k = 1_{ÞÑ k = 0, the obtained truncated characteristic equation} ˜ χ0(z) =det ( ´zJ1(1)+ A (1) 0,1 ) , (3.5)

(37)

3.2. HIERARCHICAL SPLITTING AND ASYMPTOTIC SPECTRUM is nontrivial, i.e. there exits µ such that ˜χ0(z)‰ 0. If det J1(1) ‰ 0, then ˜χ0(z) = 0 if and only

if z is an eigenvalue of the matrix (J₁(1))´1_A(1)

0,1, and hence ˜χ0(z) is nontrivial. Otherwise,

let the matrices U = [U1, U2] and V = [V1, V2] contain the left and right singular vectors

corresponding to the cokernel (U1 and V1) and image (U2 and V2) of J1(1). Condition (ND)

implies det(U˚ 1A

(1)

0,1V1) ‰ 0. Thus, using the Schur complement formula, Eq. (3.5) can be

recast as ˜ χ0(z) =det(U1˚A (1) 0,1V1)g(z), where g(z) =det ( ´zU2˚J (1) 1 V2+ U2˚A (1) 0,1V2+ U1˚A (1) 0,1V2 ( U₁˚_A(1)_0,1V₁)´1U₂˚_A(1)_0,1V₁ ) . Thus, ˜χ0(z) = 0 if and only if z is an eigenvalue of the matrix

( U˚ 2J (1) 1 V2 )´1( U˚ 2A (1) 0,1V2+ U1˚A (1) 0,1V2 ( U˚ 1A (1) 0,1V1 )´1 U˚ 2A (1) 0,1V1 )

This proves the Theorem.

3.2. Hierarchical splitting and asymptotic spectrum

This section discusses the behavior of eigenvalues λ = λε _{P Σ}ε as ε Ñ 0. Fix λ = λε and

assume that there exists a series of eigenvalues (λε₎εą0as ε decreases. Further assume that λεhave

positive real part and remain bounded away from the imaginary axis as ε Ñ 0, i.e. let λε

P Σεwith

ℜ_(λε₎

ą δ for an arbitrary δ ą 0 and all ε ą 0. Then it is easy to see that }´λε_I + A0} ď n ÿ k=1 }Ak} e´ℜ(λ ε_)σ kε´k_, _(3.6)

where }.} is an induced matrix norm, and the delay terms correspond to a small perturbation řn

k=1}Ak} exp(´ℜ(λ)ε´k) Ñ 0 as ε Ñ 0. The limiting solution λ0 =limεÑ0λεis an eigenvalue

of A0with positive real part (if it exists). This suggests that part of the spectrum (the so-called strong

unstable spectrum, see Definition 3.4, with this specific scaling property can be approximated by eigenvalues of A0with positive real part, and one expects an error, which is exponentially small in

εas ε Ñ 0 (Theorem 3.6). This motivates the following definition. Definition 3.4. Let S0 := ␣ λ_{P C | det [´λI + A}0] = 0 ( .

(38)

The set

S₀+ := S0X

␣

λ_{P C | ℜ(λ) ą 0}(, (3.7) is called the asymptotic strong unstable spectrum and the set

A₀ _{:= S}₀+_{X ˜}_S₀´ is called the asymptotic strong spectrum. Let r0 :=min

␣ |λ ´ µ| λ, µ P S0, λ‰ µ ( and r := 1 3 min ␣ r0,dist(S0, iR) (

, then the sets Σε su := Σ ε X Br(S0+), Σεss:= Σ ε X Br( ˜S0´), Σεs := Σ ε X Br(A0) (3.8)

are called strong unstable spectrum, strong stable spectrum and strong spectrum, respec-tively. The set Σε

c := ΣεzΣεs is called the pseudo-continuous spectrum.

This definition uses the fact that under certain conditions (see Theorem 3.3), it is possible that negative eigenvalues λ P Σε can be approximated by eigenvalues of A

0 with negative real part

(see Theorem 3.6 and Definition (3.1) of ˜S₀´). Note here that one has to formally exclude critical spectrum of A0on the imaginary axis, when defining the pseudo-continuous spectrum Σεc. Further,

it is worth pointing out that S0 can be obtained by formal truncation of the characteristic

equa-tion after A0, i.e. neglecting the terms including delays. Analogously, one defines the following

truncated expressions of higher order: Similar to the observation above, there is a splitting of the spectral subsets with respect to the different time scales corresponding to the hierarchy of delays. Fix 1 ď k ă n and consider the scaling function Π(k)

ε : CÑ C,

Π(k)ε (a + ib) := aε´k + ib. (3.9)

Assume there exists a series of eigenvalues (λε_{), λ}ε

P Σε such that ℜ(λε₎

Ñ 0, yet ℜ_(Π(k)_ε _(λε₎₎

Ñ γ ą 0 as ε Ñ 0. Then, since ℜ (λ) „ εk, ∆ε_(λ) has the following leading

order representation: íℑ(λ)I + A0+ k´1_ÿ j=1 Ajeíσjε k´jℑ_(λ) + Ake´σkγíσk ℑ_(λ))

as ε Ñ 0. This observation motivates the following definitions.

Definition 3.5. Define the functions χ1 : Rˆ C Ñ C, χk : Rˆ Sk´1ˆ C Ñ C, k = 2, . . . , n,

χ1(ω; Y ) := det (íωI + A0+ A1Y) , χk(ω, ϕ1, . . . , ϕk´1; Y ) := det  íωI + A0+ k´1_ÿ j=1 Ajeíϕj+ AkY  ,

(39)

3.2. HIERARCHICAL SPLITTING AND ASYMPTOTIC SPECTRUM and the corresponding asymptotic spectra

S1 := " γ+ iω_{P C | Dψ P R : χ}1 ( ω, e´σ1γ´iψ)_{= 0} * , Sk := " γ+ iω_{P C | Dψ, ϕ}1, . . . , ϕk´1 P R : χk ( ω, ϕ1, . . . , ϕk´1, e´σkγ´iψ ) = 0 * . The sets S_k+ := Sk č ␣ λ_{P C | ℜ(λ) ą 0}(

are called the asymptotic continuous τk-unstable spectrum for all k = 1, . . . , n, respectively.

Sn is called the asymptotic continuous τn-spectrum. If det An = 0, additionally define ˜χk :

R_{ˆ S}k´1_{ˆ C Ñ C, k = k, . . . , n ´ 1} ˜ χk(ω, ϕ1, . . . , ϕk´1; Y ) :=det  íωJ (k+1) 1 + A (k+1) 0,1 + k´1_ÿ j=1 A(k+1)_k,1 eíϕk _{+ A}(k+1) k,1 Y  , and if k = 1, ˜χ1 : Rˆ C Ñ C, ˜ χ1(ω; Y ) := det ( íωJ1(2)+ A (2) 0,1+ A (2) 1,1Y ) , and the corresponding asymptotic continuous spectra

˜ S1 := " γ+ iω_{P C | Dψ P R : ˜}χ1 ( ω, e´σ1γ´iψ)_{= 0} * , ˜ Sk := " γ+ iω_{P C | Dψ, ϕ}1, . . . , ϕk´1 P R : ˜χk ( ω, ϕ1, . . . , ϕk´1, e´σkγ´iψ ) = 0 * . The sets ˜ S_k´ := ˜Sk č ␣ λ_{P C | ℜ(λ) ă 0}(

are called the asymptotic continuous τk-stable spectrum for all k = k, . . . , n´ 1 respectively.

A_k _{:= S}_k+_{Y ˜}_S_k´_, ₁_{ď k ă n,} A_n _{:= S}_n

are called asymptotic continuous τk-spectra. Additionally, define the corresponding spectral

subsets Σε k,ν := " λ_{P Σ}ε_{| dist}(Π(k)ε (λ), Ak ) ă ν, |ℜ(λ)| ą ν * , k = 1, . . . , n_{´ 1,}

Multiple time-scale delay systems in mathematical biology and laser dynamics