Local bifurcations in differential equations with state-dependent delay

TL;DR

This paper extends normal form algorithms for delay differential equations to include state-dependent delays, providing a framework for analyzing local bifurcations in more complex dynamical systems.

Contribution

It develops methods to analyze bifurcations in sd-DDEs based on existing algorithms for constant delays, addressing regularity issues and invariant manifolds.

Findings

Normal form algorithms can be extended to sd-DDEs.

Invariant manifolds predicted by normal forms persist in full sd-DDEs.

Partial regularity results support bifurcation analysis in sd-DDEs.

Abstract

A common task when analysing dynamical systems is the determination of normal forms near local bifurcations of equilibria. As most of these normal forms have been classified and analysed, finding which particular class of normal form one encounters in a numerical bifurcation study guides follow-up computations. This paper builds on normal form algorithms for equilibria of delay differential equations with constant delay that were developed and implemented in DDE-Biftool recently. We show how one can extend these methods to delay-differential equations with state-dependent delay (sd-DDEs). Since higher degrees of regularity of local center manifolds are still open for sd-DDEs, we give an independent (still only partial) argument which phenomena from the truncated normal must persist in the full sd-DDE. In particular, we show that all invariant manifolds with a sufficient degree of normal hyperbolicity predicted by the normal form exist also in the full sd-DDE.

Equations

$$\dot{u}(t)=F(u_t)\text{mbox},$$

$$\dot{x}(t)$$

$$x^j$$

$$f$$

$${\tau }^j$$

$$\Vert u{\Vert }_0=\max \left\{|u(t)|:t\in [-{\tau }_{\max },0]\right\}\text{mbox}.$$

$$\Vert u{\Vert }_{D^{\ell }}$$

$$\dot{u}(t)=F(u_t)\text{mbox},$$

$$u_{t}(\theta)=u(t+\theta)\mbox{\ for $t\in[0,T]$\ and $\theta\in[-\tau_{\max},0]$.}$$

$$F(u)$$

$$u^j$$

$$D^1\times D^1\ni (u,v)\mapsto \partial F(u)v\in {\mathbb{R}}^N$$

$$\dot{x}(t)$$

$$\dot{u}(t)$$

$$\partial F(u)v=$$

$$\partial [\partial F(u,v)](w,z)=$$

$$=$$

$$-u^{\prime \prime }(u(0))w(0)v(0)$$

$$-u^{\prime }(u(0))z(0)-z(u(0))\text{mbox}.$$

$$[{\tau }_{\max },0]\ni \theta \mapsto [C^{\ell }\ni v\mapsto v(\theta )]\in \mathcal{L}(C^{\ell };\mathbb{R})$$

$$\det \Delta (\lambda )$$

$${\sigma }_c=\{{\lambda }_1,\dots ,{\lambda }_{n_c}\}=\{\lambda \in \mathbb{C}:\det \Delta (\lambda )=0\}\cap \mathrm{i}\mathbb{R}\text{mbox}.$$

$$\Delta (\lambda )=\lambda I-\sum _{j=1}^m{\partial }_{j}f(0,\dots ,0){\mathrm{e}}^{-\lambda {\tau }^j}\text{mbox},$$

$${\dot{u}}_c(t)=B^{†}{\partial }_{1}h(\cdot ;u_c(t))\text{mbox},$$

$${\dot{u}}_c=A_c^1{u}_c+\sum _{j=2}^{\ell }\frac{1}{j!}{A}_c^j[{\alpha }_j]u_c^j+o(|u_c{|}^{\ell })\text{mbox}.$$

$$\Delta (\mathrm{i}\omega )q_0=\mathrm{i}\omega q_0-A[{\mathrm{e}}^{\mathrm{i}\omega \theta }{q}_0]=0\text{mbox},$$

$$Au=\sum _{j=1}^m{A}_{j}u(-{\tau }^j)$$

$$B_1^{†}u$$

$$B_2^{†}u$$

$$h_2^{11}(\theta )$$

$$L_1=\frac{1}{2\omega }\mathrm{Re}\left(p_0\left[F_{3}qq\bar{q}+F_2\bar{q}{h}_2^{20}+F_{2}q{h}_2^{11}\right]\right)\text{mbox}.$$

$$\theta \mapsto \sum _{i=1}^{n_{\max }}{q}_i{\theta }^{{\kappa }_i}\mathrm{Re}{\mathrm{e}}^{{\lambda }_i\theta }$$

$$\frac{{\partial }^{\ell }}{\partial {\delta }^{\ell }}F{\left(\delta \sum _{j=1}^{\ell }{h}_j(\theta )u_c^j\right)|}_{\delta =0}$$

$$F_{2}vv$$

$$L_1=\frac{1}{2}\mathrm{Re}\left(\frac{2-\mathrm{i}}{1+\mathrm{i}\pi /2}\right)\approx 0.0619,$$

Full text

This paper builds on normal form algorithms for equilibria of delay differential equations with constant delay that were developed and implemented in DDE-Biftool recently. We show how one can extend these methods to delay-differential equations with state-dependent delay (sd-DDEs). Since higher degrees of regularity of local center manifolds are still open for sd-DDEs, we give an independent (still only partial) argument which phenomena from the truncated normal must persist in the full sd-DDE. In particular, we show that all invariant manifolds with a sufficient degree of normal hyperbolicity predicted by the normal form exist also in the full sd-DDE.

delay, state-dependent, local bifurcation theory

Delay-differential equations (DDEs) arise frequently in models where the evolution of the system depends also on its values in the past. Typical examples arise in control (delays in feedback loops), optics (delayed feedback effects from external light reflections), mechanical engineering (effects from previous rotations in turning processes), or Earth sciences (El Niño caused by delayed feedback from waves across oceans).

The typical approach to studying DDEs is to consider them as a dynamical systems for which the state is a history segment (in our case on a bounded history interval). Several mathematical problems occur when the length of the delay depends on the state of the system, called sd-DDEs. In this case the state of the dynamical system at time $t$ does not depend smoothly on its initial condition. This makes many of the standard tools of dynamical systems theory inapplicable at first sight. In particular normal form theory requires expansion of the right-hand side to higher orders.

This paper demonstrates that normal forms can still be computed for a general class of sd-DDEs with discrete delays. We show that the computational procedure developed by Janssens, Wage, Bosschaert and KuznetsovJanssens (2010); Wage (2014); Bosschaert (2016); Bosschaert, Janssens, and Kuznetsov (2017) for DDEs with constant delays can be generalized to sd-DDEs. We also give a justification for the computed normal forms, explaining why all normally hyperbolic manifolds present in the normal form also appear in the full sd-DDE. The justification is based on an approach recently taken by Humphries et alHumphries, Calleja, and Krauskopf (2016) in a numerical bifurcation study of a prototypical sd-DDE.

I Introduction

Delay-differential equations (DDEs) are a class of differential equations where the derivative at the current time $t$ may depend on any value of the state in the past. This paper focusses on those case where the dependence is on states from a limited time interval $[t-\tau_{\max},t]$ in the past. They are a particularly common and well-studied subclass of so-called functional-differential equationsHale and Verduyn Lunel (1993); Diekmann et al. (1995). Mathematically, DDEs are dynamical systems with an infinite-dimensional phase space, since the appropriate initial value is a prescribed piece of history of the physical variable on an interval $[-\tau_{\max},0]$. A typical choice of phase space is the space of $n$-dimensional continuous functions on $[-\tau_{\max},0]$, written as $C^{0}([-\tau_{\max},0];\mathbb{R}^{n})$ with the maximum norm (short $C^{0}$). The right-hand side is given by a functional $F:C^{0}\to\mathbb{R}^{n}$. An example is $F(u)=-u(-\tau)$ for a fixed $\tau>0$ and functions $u$ close to [math] in $C^{0}$. Then one will write the differential equation $\dot{u}(t)=-u(t-\tau)$ as

[TABLE]

where the subscript $t$ indicates a time-shifted history interval. So, for a function $u:[-\tau_{\max},T]\to\mathbb{R}^{n}$ and $t\in[0,T]$, $u_{t}$ is a function on $[-\tau_{\max},0]$ defined by $u_{t}(\theta)=u(t+\theta)$.

There is mathematically a large difference between DDEs with constant delays and DDEs with state-dependent delays. For constant delays, a framework that poses DDEs as abstract ODE has been developed by Hale & Verduyn-LunelHale and Verduyn Lunel (1993) and Diekmann et alDiekmann et al. (1995). In this framework DDEs of the type $\dot{u}(t)=F(u_{t})$ are smooth dynamical systems on the phase space $C^{0}$. That is, the time-$t$ map $u_{0}\mapsto u_{t}$ for fixed $t$, mapping the initial condition $u_{0}\in C^{0}$ to the solution $u_{t}\in C^{0}$ at time $t$, is smooth. The smoothness of the time-$t$ map follows from the smoothness of the functional $F:C^{0}\to\mathbb{R}^{n}$.

This is in contrast to the case when the functional $F$ involves state-dependent delays. We refer to this type of DDEs as DDEs with state-dependent delays (short sd-DDEs). An example is the differential equation $\dot{u}(t)=p-u(t+u(t))$ for fixed parameter $p$, for which the functional $F$ has the form $F:u\mapsto p-u(u(0))$ (for $u$ close to $p$ and $p<0$). The derivative of the right-hand side $F$ with respect to its argument $u$ is $\partial F(u)v=-v(u(0))-u^{\prime}(u(0))v(0)$ if it exists. Thus, it is undefined for $u\in C^{0}$ that are not differentiable. This has the consequence that the standard theory from textbooksHale and Verduyn Lunel (1993); Diekmann et al. (1995) for DDEs is not applicable. The currently most practical statements (for dynamical systems theory) about the regularity of the time-$t$ map with respect to its initial value are by HartungHartung (2011) and WaltherWalther (2003). They are much more restricted, achieving at best continuous differentiability (once) of the time-$t$ map. A review by Hartung et al from 2006Hartung et al. (2006) presents a snapshot of developments regarding general existence and regularity theory. Section II.2 summarizes the most relevant results.

Applications and numerical software

In parallel to developments in the theory of sd-DDEs, computational tools have been created to help solving practical problems arising in engineering and science. The review by Hartung et alHartung et al. (2006) lists a few classical applications such as control by echo locationWalther (2002), models for cutting processes with a finite tool stiffness in directions tangential to the rotating surfaceInsperger, Barton, and Stépán (2008); Insperger, Stépán, and Turi (2007) and the electromagnetic two-body problemLuca et al. (2010). Other examples are time-delayed feedback control where the time-delay is adjusted dynamicallyPyragas and Pyragas (2011), and models for granulopoiesisCraig, Humphries, and M.C.Mackey (2016).

Two common tasks to be performed numerically in applications are initial-value problem solving (a black-box solver for sd-DDEs including neutral terms is RADAR5Guglielmi and Hairer (2001)) and numerical bifurcation analysis. Numerical bifurcation analysis tracks branches of equilibria (constant solutions of $F(u)=0$), periodic orbits (time-periodic solutions of $\dot{u}(t)=F(u_{t})$) and their bifurcations and linear stability. Equilibria of sd-DDEs are given by algebraic equations and periodic boundary-value problems can be reduced to equivalent systems of smooth algebraic equationsSieber (2012). Thus, numerical computations of these are feasible in principle and have been implemented in DDE-BiftoolEngelborghs, Luzyanina, and Roose (2002); Engelborghs, Luzyanina, and Samaey (2001); Sieber et al. . Its capabilities for sd-DDEs with discrete delays (as described in Section II.1) include:

• •

continuation of families of equilibria and computation of their stability (present since version 2.0);

• •

continuation of codimension-one bifurcations of equilibria (Hopf bifurcations and saddle-node bifurcations, present since version 2.0);

• •

continuation of periodic orbits in one parameter and computation of their stability (present since version 2.0, completed for the class of sd-DDEs with discrete delays described in Section II.1 in version 3.0);

• •

continuation of local codimension-one bifurcations of periodic orbits (saddle-node bifurcations, period doubling bifurcations and torus bifurcations, present since version 3.0);

Normal forms of local bifurcations

This paper gives the background on how direct normal form computations for codimension-one and -two bifurcations of equilibria have been added for sd-DDEs to the general sd-DDE capabilities. The procedures are based on the corresponding code and work by Kuznetsov, Janssens, Wage and BosschaertJanssens (2010); Wage (2014); Bosschaert (2016); Bosschaert, Janssens, and Kuznetsov (2017) for constant-delay DDEs. Section III reviews these recent developments for constant delays. Appendix A gives more details.

Normal form computations help classify all generic (up to codimension two) bifurcations into a finite number of well-studied cases. Thus, they help the systematic numerical exploration in applications. For example, when a Hopf bifurcation is detected, one may compute the so-called Lyapunov coefficient which determines to which side the periodic orbits branch off from the equilibrium (that is, whether the Hopf bifurcation is sub- or supercritical, or, using the terms coined in engineering, safe or dangerousThompson and Stewart (2002)). The illustrative example of a linear position control problem with state-dependent delay in Section V shows a typical scenario.

Similarly, when following a Hopf bifurcation in two parameters, one typically encounters crossings with other Hopf bifurcations (a common scenario for DDEs). At these so-called Hopf-Hopf interaction points various branches of secondary bifurcations can be expected depending on the normal form of the Hopf-Hopf interaction. Humphries et alHumphries, Calleja, and Krauskopf (2016) studied bifurcations of a scalar sd-DDE in detail. They encountered several Hopf-Hopf interactions, derived the normal form on paper, and then followed the predicted secondary bifurcations, which turned out to exist in the expected directions.

Justification of normal form expansion in sd-DDEs

The normal form of most codimension-one and -two bifurcations depends on expansion terms of order higher than one. Expansion to this degree is not immediately justifiable for sd-DDEs since the time-$t$ map of sd-DDEs is only continuously differentiable once. For ordinary differential equations (ODE), there are precise statements about the relation between the phase portraits and their bifurcations in truncated normal forms and the full dynamical system (they depend on the particular bifurcation)Guckenheimer and Holmes (1990); Kuznetsov (2004). To obtain the same statements for sd-DDEs one needs that local center manifolds near equilibria are smooth to the degree required for the expansion terms in the normal form (for example, to third order for the Hopf bifurcation). A local center manifold near an equilibrium in a (sd-)DDE has the form of a graph $h:\mathbb{R}^{n_{c}}\to C^{0}([-\tau_{\max},0];\mathbb{R}^{n})$. Here $n_{c}$ is the number of eigenvalues (counted with multiplicity) of the linearized DDE on the imaginary axis, and the domain of $h$ is a coordinate representation of the corresponding eigenspace. The smoothness requirement for $h$ refers to two things. First, each element of the center manifold has to be smooth with respect to its argument (time), so $h(u_{c})\in C^{\ell}([-\tau_{\max},0];\mathbb{R}^{n})$ (the space of $\ell$ times continuously differentiable functions). Second, the graph $h$ has to be a smooth map of its argument $u_{c}\in\mathbb{R}^{n_{c}}$. Smoothness of local center manifolds has not been proven rigorously yet for degrees greater than one. Stumpf Stumpf (2011) gives a proof of continuous differentiability of center-unstable manifolds, and shows that it attracts exponentially all those solutions that stay near the equilibriumStumpf (2015a). However, we prove in Section IV.2 that many phenomena predicted by the normal form must also be present in the sd-DDE. The statement is not as strong as its classical ODE counterpart such that the availability of numerical normal form computations provides a motivation to investigate the smoothness of local center manifolds rigorously.

II DDEs with state-dependent delays

II.1 Discrete state-dependent delays

DDE-Biftool is able to perform bifurcation analysis on a class of $n$-dimensional systems of delay differential equations with $m-1$ discrete state-dependent delays (sd-DDEs) of the following form:

[TABLE]

The integers $n\geq 1$ (physical space dimension), $m\geq 1$ (number of delays) and $n_{p}\geq 0$ (number of parameters) are arbitrary. It uses the convention that $\tau^{1}=0$ and assumes that the functions

[TABLE]

are smooth. The construction (1)–(2) permits arbitrary levels of nesting in the delayed arguments of $x$. DDE-Biftool does not require an explicit value for the maximal delay. It computes equilibria and periodic orbits such that the trajectory $x(t)$ is always compact.

In sections with theoretical considerations we may assume that $n_{p}=0$ without loss of generality by incorporating the parameters into the state (appending the equation $\dot{p}=0$ to (1) and increasing $n$ to $n+n_{p}$).

II.2 General functional differential equations (FDEs) —

Review of basic properties

Notation and assumptions on the right-hand side

In the following sections we will use the abbreviation that $C^{0}$ (or just $C$) is the space $C([-\tau_{\max},0];\mathbb{R}^{n}$) of continuous functions on the interval $[-\tau_{\max},0]$ into $\mathbb{R}^{n}$ with the norm

[TABLE]

Similarly, for any space $D$ of functions on an interval $I\subset\mathbb{R}$ and integer $\ell>0$, we denote the subspace $D^{\ell}$ as the space of functions which have a $\ell$th derivative in $D$. Their respective norms are

[TABLE]

We also use the phrase, for example, “$f$ is $C^{\ell}$ ” for $f$ being $\ell$ times continuously differentiable in all its arguments.

Basic existence and regularity theory for solutions of sd-DDEs has been developed for differential equations in the form

[TABLE]

where $F:C([-\tau_{\max},0];\mathbb{R}^{n})\to\mathbb{R}^{n}$ is a continuous nonlinear functional Hartung et al. (2006). For a function $u:[-\tau_{\max},T]\to\mathbb{R}^{n}$ the notation $u_{t}$ refers to a time shift of $u$ back to a function on the interval $[-\tau_{\max},0]$:

[TABLE]

For the type of equations that can be treated with DDE-Biftool the functional $F$ (incorporating parameters into the state variables) has the form

[TABLE]

If the coefficient functions $f$ and $\tau^{j}$ are $\ell$ times continuously differentiable, we call such a functional $F$ a functional with $C^{\ell}$ coefficients and $m$ state-dependent discrete delays less than $\tau_{\max}$.

The general conditions on $F$ to ensure existence and regularity of solutions vary between different papers. A set of conditions that covers functionals $F$ with discrete state-dependent delays and $C^{\ell}$ coefficients and satisfies the assumptions in many fundamental papers is mild differentiability. Consider a continuous functional $F:D\to\mathbb{R}^{N}$ for some $N\geq 1$ and some $D$ that is a subspace of $C^{0}(I;\mathbb{R}^{N})$ for some interval $I\subset\mathbb{R}$. For mild differentiablity of $F$ we require the following two conditions.

1. (S1)

The functional $F$ is continuously differentiable when restricted to the subspace $D^{1}$. We denote its derivative by $\partial F:D^{1}\to{\cal L}(D^{1};\mathbb{R}^{N})$. 2. (S2)

The map

[TABLE]

can be extended continuously to the space $D^{1}\times D$.

We put the argument $v$ of $\partial F$ outside of the bracket to emphasize that $\partial F$ is linear in $v$. Since $\partial F:D^{1}\times D\to\mathbb{R}^{N}$ is continuous, we can apply the definition for mild differentiability recursively, treating the pair $(u,v)\in D^{1}\times D$ as the single argument of $\partial F$. This leads naturally to the definition that a functional $F:D\to\mathbb{R}^{N}$ is $\ell$ times mildly differentiable if

1. (S3)

$\partial F:(D^{1}\times D)\to\mathbb{R}^{N}$ is $\ell-1$ times mildly differentiable.

Scalar illustrative example

An illustrative example is the sd-DDE

[TABLE]

This corresponds to the choice $f(x,y,p)=p-y$ and $\tau^{2}(x,p)=-x$ in (6)–(7) (using letters $x$ and $y$ in the arguments of $f$ instead of superscripts to avoid confusion with powers), where we keep $p=-\pi/2$ fixed for illustration initially. So, $F$ is a functional with $2$ delays and $C^{\infty}$ coefficients. The first two derivatives of this functional $F$ are

[TABLE]

Note how the second derivative includes differentiation of the first derivative with respect to $v$ according to our convention such that it has $4$ arguments (generally, the $\ell$th derivative will have $2^{\ell+1}$ arguments). We reserve the notation $\partial^{j}F(u)$ for the usual $j$-linear form. The above expressions show that the $\ell$th derivative of $F$ depends on the lowest $\ell$ derivatives of $u$, on the lowest $(\ell-1)$ derivatives of the deviation $v$ and $w$, and only on the values of $z$. So, $\partial^{1}F$ is continuous in $C^{1}\times C^{0}$ and $\partial[\partial F]$ is continuous in $(C^{2}\times C^{1}\rangle\times(C^{1}\times C^{0})$. Moreover, the map $u\mapsto\partial F(u,\cdot)$ is continuous as a map, mapping $u\in C^{1}$ into the space ${\cal L}(C^{1};\mathbb{R})$ of linear functionals from $C^{1}$ into $\mathbb{R}$, but not as a map into the space ${\cal L}(C^{0};\mathbb{R})$ of linear functionals from $C$ into $\mathbb{R}$. The reason for this discontinuity is the second term $-v(u(0))$: the map

[TABLE]

is only continuous in $\theta$ if $\ell\geq 1$. Mild differentiability of second order requires that $(u,v)\mapsto\partial[\partial F(u,v)](\cdot,\cdot)\in{\cal L}(C^{2}\times C^{1};\mathbb{R})$ is continuous, which is the case for the right-hand side $F$ in example (II.2).

The example illustrates that the assumptions of mild differentiability permit dependence of the delays on the state. We note that for varying $p$, we have to include the equation $\dot{p}=0$. The combined system also satisfies mild differentiability to all orders. Equation (II.2) has an equilibrium at $u=p$, which loses its stability in a Hopf bifurcation at $p=-\pi/2$. We will use the above example (II.2) to illustrate various technical assumptions and difficulties in the following sections. For example, the form of the first derivative of $F$ in (II.2) implies that $F$ is not locally Lipschitz continuous in $C^{0}$.

Basic results on solutions of sd-DDEs

Successive differentiation and application of the chain rule imply that functionals $F$ with discrete delays and $C^{\ell}$ coefficients (in the form of (6)–(7)) satisfy assumptions (S1–S3) up to the order $\ell$. Thus, all of the following basic results apply to this class of sd-DDEs with discrete delays.

WaltherWalther (2003, 2004) proved that initial value problems (IVPs) have a unique solution $u$ for all times $t$, or the solution blows up in finite time, if the initial value $u_{0}$ lies in the manifold ${\cal M}_{F}=\{u\in C^{1}:u^{\prime}(0)=F(u)\}\subset C^{1}\mbox{.}$ Moreover, for times $t$ before blow-up the map ${\cal M}_{F}\ni u_{0}\mapsto u_{t}\in{\cal M}_{F}$ is continuously differentiable. Thus, sd-DDEs generate a $C^{1}$ semiflow (time-$t$ maps) in suitable open subsets of ${\cal M}_{F}$ (for example, in a sufficiently small neighborhood of equilibria or periodic orbits). Hence, Walther’s result immediately implies that the principle of linearized stability applies with respect to perturbations in ${\cal M}_{F}$, in particular to equilibriaStumpf (2015b) and periodic orbits. This basic existence result requires only first-order mild differentiability (a slightly weaker version of them, since continuity of $F$ in $C^{0}$ is not needed Walther (2003); Hartung et al. (2006)). Krisztin Krisztin (2003) proved that the unstable manifold of equilibria is a $C^{\ell}$ graph for $\ell$ times mildly differentiable right-hand sides, using a slightly different (possibly equivalent) definition of mild differentiability for orders greater than $1$. Based on Walther’s semiflow results, StumpfStumpf (2011); Stumpf et al. (2016) proved the existence and attractivity of $C^{1}$ local center-unstable and center manifolds near equilibria. Alternative proofs are given by KrisztinKrisztin (2006a, b). Furthermore, the assumptions (S1–S3) imply that periodic boundary-value problems are equivalent to finite-dimensional smooth systems of algebraic equations for a sufficiently large number of first Fourier coefficientsSieber (2012). This equivalence permits us to perform a classical Lyapunov Schmidt reduction near equilibria $u_{*}$ for which the characteristic matrix $\Delta(\lambda)\in\mathbb{C}^{n\times n}$, defined by $\Delta(\lambda)q=\lambda q-\partial F(u_{*})[\theta\mapsto q\exp(\lambda\theta)]$ has a single pair of roots on the imaginary axis. Consequently, the classical Hopf bifurcation theorem about a family of periodic orbits branching off from $u_{*}$ is validEichmann (2006); Sieber (2012), including formulas determining criticality of the Hopf bifurcation. More generally, the reduction of periodic boundary value problems to smooth algebraic equations implies that all objects computed by DDE-Biftool depend as expected on parameters and the right-hand side such that they can be computed using standard numerical discretizationsSieber (2012). This includes branches of periodic orbits in parameter-dependent systems, the variational problems for folds, period doublings and torus bifurcationsSieber (2013). Statements about periodic orbit families branching off at period doublings and resonant torus bifurcations (in resonance tongues, first computational demonstrations for DDEs were for an El-Ninõ modelTziperman et al. (1998); Keane, Krauskopf, and Postlethwaite (2015, 2016); Krauskopf and Sieber (2014)) follow in a similar way from a Lyapunov-Schmidt reduction as the Hopf bifurcation statement.

III Normal form computations in DDEs with constant delays — Review

Recent work by Kuznetsov, Janssens, Wage and BosschaertJanssens (2010); Wage (2014); Bosschaert (2016); Bosschaert, Janssens, and Kuznetsov (2017) has developed and implemented expressions for the normal form coefficients of local bifurcations in DDEs with constant delays. For discrete delays, this corresponds to the case where the delay functions $\tau^{j}$ in (7) are all constant (e.g., parameters) independent of the state. Their procedure follows closely the methods originally developed for ODEsKuznetsov (1999) (and is in principle applicable to other abstract ODEsvan Gils et al. (2013)). They assume that the DDE $\dot{u}(t)=F(u_{t})$ has an equilibrium at $u_{*}$. For our notation we assume $F(0)=0$, and denote the first derivative of the right-hand side $F:C^{0}\to\mathbb{R}^{n}$ in [math] by $A=\partial F(0)\in{\cal L}(C^{0};\mathbb{R}^{n})$.

III.1 Linear stability and center manifold

The matrix $\Delta(\lambda)\in\mathbb{C}^{n\times n}$ defined by $\Delta(\lambda)q=\lambda q-A[\theta\mapsto q\exp(\lambda\theta)]$ for $q\in\mathbb{C}^{n}$ is called the characteristic matrix. We assume that the characteristic equation

[TABLE]

has $n_{c}$ roots (including multiplicity) on the imaginary axis:

[TABLE]

For the type of functionals $F$ that DDE-Biftool treats, $\Delta(\lambda)$ is given by

[TABLE]

where for constant delays the $\tau^{j}$ are parameters, while for state-dependent delays, the $\tau^{j}$ are evaluated at the equilibrium [math]. The corresponding eigenvectors are in $C^{\infty}$, and have the form $\theta\mapsto q\exp(\lambda\theta)$. The generalized eigenvectors (also in $C^{\infty}$ if present) have the form $\theta\mapsto\sum_{j=0}^{j_{\max}}q^{j}\theta^{j}\exp(\lambda\theta)$, where $j_{\max}+1$ is the length of the Jordan chain and $q^{0},\ldots,q^{j_{\max}}$ are in $\mathbb{C}^{n}$. Let $B=\{b_{1},\ldots,b_{n_{c}}\}$ be a basis of real functions of the linear center subspace $U_{c}=\operatorname{span}B$ of $\dot{u}=Au_{t}$ in $C^{0}$, and let $B^{\dagger}:C^{0}\to\mathbb{R}^{n_{c}}$ be such that $B^{\dagger}B=I$ in $\mathbb{R}^{n_{c}}$ and $BB^{\dagger}$ is a spectral projection onto $\operatorname{span}B$ (see (29)–(30) in the Appendix for a concrete expression based on the resolvent formalism).

Center manifold for constant delays

For DDEs with constant discrete delays ($\tau^{j}=\mathrm{const}$ in (4)) the time-$t$ map $C\ni u_{0}\mapsto u_{t}\in C$ is as smoothDiekmann et al. (1995); Hale and Verduyn Lunel (1993) as the right-hand side $f:\mathbb{R}^{n\times m}\mapsto\mathbb{R}^{n}$ in (1). The reason is that, for those $f$, the right-hand side as a map $F:C^{0}\to\mathbb{R}^{n}$ is smooth. Hence, in a ball $B_{r}(0)$ around [math] with sufficiently small radius $r$ a smooth center manifold of dimension $n_{c}$, $h:B_{r}(0)\subset\mathbb{R}^{n_{c}}\to C^{0}$ exists.

More precisely, let us assume that the right-hand side coefficient function $f$ in (1) is at least $\ell$ times continuously differentiable. Then we can find a radius $r>0$ such that the invariant graph $h:B_{r}(0)\subset\mathbb{R}^{n_{c}}\to C^{\ell}$ is $\ell$ times differentiableDiekmann et al. (1995); Hale and Verduyn Lunel (1993). We write the graph as $h(\theta;u_{c})$, putting the argument of the function $h(u_{c})$ in $C^{\ell}$ first. For any initial condition $u_{0}(\theta)=h(\theta;u_{c}^{0})$ ($u_{c}^{0}\in B_{r}(0)$) on the graph, $u_{t}(\theta)$ equals $h(\theta;u_{c}(t))$, where

[TABLE]

and $u_{c}(0)=u_{c}^{0}$, as long as $|u_{c}(t)|\leq r$.

III.2 Normal form computation

Assuming that the right-hand side $F$ and the center manifold $h$ are smooth up to a desired order $\ell$ (as is the case for constant delays), it is known that the flow on the local center manifold can be brought into a normal form up to order $\ell$, such that the flow on the center manifold $\dot{u}_{c}=B^{\dagger}\partial_{1}h(\cdot;u_{c})$ has a given expansion

[TABLE]

Equation (10) is an ODE for $u_{c}\in\mathbb{R}^{n_{c}}$. All derivatives up to order $\ell$ of the remainder $o(|u_{c}|^{\ell})$ are smaller than the corresponding derivatives of the lower-order terms for all small $|u_{c}|$. All of the $j$-linear coefficients $A_{c}^{j}$ depend only on the type of equilibrium (which local bifurcation?), except for the still-to-be-determined normal form parameters $\alpha_{j}$ at each order $j>1$. The linear coefficients $A_{c}^{1}$ are uniquely determined by $B$ and $B^{\dagger}$: $A_{c}^{1}=B^{\dagger}B^{\prime}$, where $B^{\prime}$ is the derivative of $B$ with respect to the space variable $\theta$. There exists a $C^{\ell}$-smooth coordinate change in $\mathbb{R}^{n_{c}}$ that transforms the ODE (9), describing the semiflow of the DDE restricted to its local center manifold $h$, into Equation (10) (this is called smooth local equivalence).

Normal form computations are concerned with the computations of these unknown coefficients $\alpha_{j}$ and, if desired, the expansion coefficients $h_{j}(\theta)=\partial^{j}_{2}h(\theta;0)$ of the center manifold. Inputs are the expansion coefficients $F_{j}=\partial^{j}F(0)$ (also $j$-linear forms) of the right-hand side of the DDE, and the general parametric normal form expansion coefficients $A_{c}^{j}[\cdot]$, which depend on the type of the bifurcation investigated (e.g., Hopf bifurcation and degenerate Hopf bifurcation in the example in Section V). The procedure for computing the coefficients $\alpha_{j}$, as outlined for ODEs by KuznetsovKuznetsov (1999), and adapted to DDEs recentlyJanssens (2010); Wage (2014); Bosschaert (2016); Bosschaert, Janssens, and Kuznetsov (2017), is summarized in Section A in the appendix.

The invariance of $h$ gives at each order a linear system of equations for the expansion coefficients $h_{j}(0)$ of the center manifold at $\theta=0$. The system depends also linearly on $\alpha_{j}$ (if at order $j$ a normal form coefficient is present). The coefficients of the linear system for $h_{j}(0)$ and $\alpha_{j}$ depend only on $A$ (same as $F_{1}$), the linear part of $F$. At each order $j$, the coefficient $\alpha_{j}$ is determined by the Fredholm alternative as the unique value for which the linear system is solvable for $h_{j}(0)$.

III.3 General example — Hopf bifurcation

A typical result of the procedure is the normal form coefficient $L_{1}$ (which would be the real part of $\alpha_{3}$, divided by $\omega$) for the Hopf bifurcationWage (2014), as implemented in DDE-BiftoolSieber et al. ; Wage (2014); Bosschaert (2016); Bosschaert, Janssens, and Kuznetsov (2017). Suppose the linearized DDE $\dot{u}=\partial F(0)u_{t}=Au_{t}$ has a purely imaginary eigenvalue pair $\pm\mathrm{i}\omega$, with the eigenvector $q=q_{0}\mathrm{e}^{\mathrm{i}\omega\theta}$ and its complex conjugate $\bar{q}=\bar{q}_{0}\mathrm{e}^{-\mathrm{i}\omega}$. That is,

[TABLE]

and $\pm\mathrm{i}\omega$ are the only roots of $\det\Delta(\cdot)$ on the imaginary axis. For notational convenience one chooses as basis $B=h_{1}$ of the center subspace of $C^{0}$ the vectors $\{q,\bar{q}\}$, thus using complex notation instead of, for example, $\{\operatorname{Re}q,\operatorname{Im}q\}$. The projection $B^{\dagger}$ is given by the normalized adjoint eigenvector $p$ for $\mathrm{i}\omega$ and its complex conjugate $\bar{p}$. The general expression for adoint eigenvectors is given by Diekmann et alDiekmann et al. (1995). For the particular case, where the linear functional $A$ has the form

[TABLE]

(as arising in problems treatable with DDE-Biftool) and the critical spectrum consists of simple eigenvalues $\pm\mathrm{i}\omega$, the projection is of the form

[TABLE]

The $C^{1\times n}$ vector $p_{0}$ is given by $p_{0}\Delta(\mathrm{i}\omega)=0$ and (after normalization) $p_{0}\Delta^{\prime}(\mathrm{i}\omega)q_{0}=1$. At order $2$ the linear system for the coefficients of the center manifold is regular (thus, $\alpha_{2}$ is empty). Solving it yields

[TABLE]

(the remaining coefficient is $h_{2}^{02}=\bar{h}_{2}^{20}$). At order $3$, there is a single complex coefficient ($\alpha_{3}\in\mathbb{C}$ of which the real part is the coefficient $\omega L_{1}$) such that:

[TABLE]

If the coefficient $L_{1}$ is non-zero the Hopf bifurcation is non-degenerate (subcritical if $L_{1}>0$, supercritical if $L_{1}<0$).

IV Extension to DDEs with state-dependent delays

Several observations about the normal form reduction imply that at least the computational procedure can be extended to DDEs with state-dependent delays (sd-DDEs).

The procedure described in section III.2 requires the expansion coefficients $F_{j}$ of the nonlinearity $F$ up to the desired order (often at least $3$). However, we observe that the derivatives are applied only to deviations that are expansion coefficients of the center manifold, $(\theta,u_{c})\mapsto h_{j}(\theta)u_{c}^{j}$, where $\theta$ is the history variable and $u_{c}$ is the deviation along the center manifold. At each order $j$, the unknown coefficient $h_{j}(\theta)$ is a solution of the linear ODEs (35) (see Appendix) with constant coefficients and an inhomogeneity that is a linear combination of $h_{k}(\theta)$ from lower orders ($k<j$). The basis of the linear center subspace (called $B$ in the previous section and equal to $h_{1}$) consists of functions of the form of a finite sum

[TABLE]

of some length $n_{\max}$ with $n_{\max}$ non-negative integer powers $\kappa_{i}$ of $\theta$ (possibly, some $\kappa_{i}=0$), and complex exponents $\lambda_{i}$. Therefore the ODE (35) defining the coefficients $h_{j}(\theta)$ implies that all center manifold expansion coefficients have the form (12). Hence, they are smooth in $\theta$ such that the functional $F$ can be differentiated in the equilibrium in the direction of $\sum_{j=1}^{\ell}h_{j}(\theta)u_{c}^{j}$ for all $\ell$ and all $u_{c}\in\mathbb{R}^{n_{c}}$.

The derivative of expressions of the form (12) is known analytically such that a user routine computing the directional derivative

[TABLE]

can rely on all derivatives of the argument of $F$ with respect to $\theta$. Similarly, finite-difference approximations of the derivative with respect to $\delta$ are known to converge. Both approaches are experimentally supported in the current development version of DDE-BiftoolSieber et al. . Section V will illustrate their use for a position control problem.

IV.1 Illustration for Hopf bifurcation in sd-DDE (II.2)

For the example $\dot{x}(t)=p-x(t+x(t))$ the characteristic matrix $\Delta(\lambda)$ of the linearization in the equilibrium $x_{*}=p$ has the form $\Delta(\lambda)=\lambda-\mathrm{e}^{\lambda p}$, which has a Hopf bifurcation with critical eigenvalue $\mathrm{i}\omega=\mathrm{i}$ at $p=-\pi/2$. Thus, the right eigenvector is $q(\theta)=\mathrm{e}^{\mathrm{i}\theta}$, and the left eigenvector $p$ will be scaled such that $p(0)\Delta^{\prime}(\mathrm{i})q(0)=1$. Thus, $p_{0}=1/(1+\mathrm{i}\pi/2)\approx\mathtt{0.2884-0.4530i}$. The second and third directional derivatives of $F(u)=p-u(u(0))$ in [math] along a fixed direction $v$ are

[TABLE]

The mixed derivatives $F_{2}q\bar{q}$ and $F_{3}qq\bar{q}$ can be constructed from directional derivatives using the polarization identity (DDE-Biftool’s implementation uses this approach). Following the procedure for the general Hopf normal form in Section III.3 we compute $h_{2}^{20}(\theta)=(0.4+0.8\mathrm{i})\mathrm{e}^{2\mathrm{i}\theta}$ and $h_{2}^{11}(\theta)=-4$ (constant), resulting in a Lyapunov coefficient

[TABLE]

which indicates that the Hopf bifurcation is subcritical (dangerous) for this example.

IV.2 Smoothness of coefficients

A combination of previous results provides an immediate partial justification for the normal forms computed with the procedure given by Kuznetsov et alJanssens (2010); Wage (2014); Bosschaert (2016); Bosschaert, Janssens, and Kuznetsov (2017) and summarized in Section III. First of all, trajectories of sd-DDEs become more regular over time. This effect is well known for DDEs with constant delays, but also holds for sd-DDEs. The general proof requires the precise definition of order-$\ell$ mild differentiability. We formulate the the statement here for DDEs with discrete state-dependent delays.

Proposition IV.1** (Smoothness for large times)**

Assume that $F$ is a functional with $C^{\ell}$ coefficients and $m$ discrete state-dependent delays (of the form (6)–(7)) less than $\tau_{\max}$. Let $u(t)$ with $t\in[-\tau_{\max},T]$ be a solution of $\dot{u}(t)=F(u_{t})$ with $u_{0}\in C^{1}$ and $u_{0}^{\prime}(0)=F(u_{0})$. Then $u_{t}\in C^{\ell}$ if $t\geq\ell\tau_{\max}$. The $\ell$th derivative $u^{(\ell)}$ satisfies a (differential) equation of the form

[TABLE]

where $F^{\ell}$ has $C^{0}$ coefficients and $m_{\ell}=(m+1)^{\ell-1}m$ discrete delays less than $\ell\tau_{\max}$.

Proof

We show this statement (inductively). For $\ell=1$ the statement follows from the differential equation with $F^{1}=F$ ($f^{1}=f$ and $m=m_{1}$). Assume that we have for $t\geq\ell\tau_{\max}$

[TABLE]

where $u^{j}=u(t-\tau^{j}_{\ell}(u^{1},\ldots,u^{j-1}))$ and all $\tau_{\ell}^{j}\leq\ell\tau_{\max}$ (for $\ell=1$, $\tau_{1}^{j}=\tau^{j}$ for $j=1,\ldots,m$). Thus, for $t\geq(\ell+1)\tau_{\max}$ $u_{t}(\theta)$ is $C^{1}$ for all $\tau\in[-\ell\tau_{\max},0]$. Consequently, the right-hand side of (14) is differentiable with respect to time for $t>(\ell+1)\tau_{\max}$ (and, hence, the left-hand side). Its derivative is

[TABLE]

For $j=1$ the above expression (16) for $V^{j}$ equals $\dot{u}(t-\tau_{\ell}^{1})=\dot{u}(t)$. We replace $\dot{u}(t-\tau_{\ell}^{j})$ in (16) with $F^{1}(u_{t-\tau_{\ell}^{j}})$ such that

[TABLE]

where for $k=1,\ldots,m_{1}$

[TABLE]

We see that the right-hand side in (15) is a functional $F^{\ell+1}$ of the same form as $F^{\ell}$, but where $f^{\ell+1}$ has $m_{\ell}+m_{1}m_{\ell}$ arguments such that we have $m_{\ell}+m_{1}m_{\ell}$ delays. Those delays are $\tau_{\ell}^{1}$,…, $\tau_{\ell}^{m_{l}}$ and for $j=m_{\ell}+(i-1)m_{\ell}+k$ ($i=1\ldots,m_{\ell}$, $k=1,\ldots,m_{1}$)

[TABLE]

which are all less than $(\ell+1)\tau_{\max}$. Hence, $u^{(\ell+1)}$ exists for $t>(\ell+1)\tau_{\max}$ and satisfies $u^{(\ell+1)}(t)=F^{\ell+1}(u_{t})$. (End of proof of Proposition IV.1)

Since $F^{\ell}(0)=0$, and the coefficients $f^{j}$ and $\tau^{k}_{j}$ are still at least $C^{1}$ for all $j\leq\ell$ (we have differentated only $\ell-1$ times), we have for all $u_{0}\in C^{1}$ sufficiently close to [math] that

[TABLE]

for $t\geq\ell\tau_{\max}$ and all $j\leq\ell$ and some constant $C(t)>0$.

A local center-unstable manifold $h$ is exists and is continuously differentiable for functionals $F$ with $C^{1}$ coefficients and discrete state-dependent delays, according to Stumpf Stumpf (2011). Consequently, if $F(0)=0$ and the critical spectrum $\sigma_{c}$ of $\dot{u}=\partial F(0)u_{t}$ is not empty, a continuously differentiable local center manifold $h$ exists, too (applying the standard local center manifold theorem to the ODE with $C^{1}$-smooth coefficients that one obtains by restricting the sd-DDE onto its local center-unstable manifold, see also Stumpf’s or Krisztin’s argumentsKrisztin (2006a, b); Stumpf et al. (2016)). A simple backwards extension and Proposition IV.1 permit us to conclude that all elements of the local center manifold $h$ are in $C^{\ell}$:

Lemma IV.2** (Smoothness on center manifold)**

Assume that $F$ is a functional with $C^{\ell}$ coefficients and discrete state-dependent delays (of the form (6)–(7)), with $F(0)=0$, a center subspace $\operatorname{span}B$ of $\Delta(\lambda)=\lambda I-\partial F(0)[\theta\mapsto\exp(\lambda\theta)]$ of dimension $n_{c}$ and a continuously differentiable local center manifold $h:B_{r}(0)\subset\mathbb{R}^{n_{c}}\to C^{1}$, defined in a ball $B_{r}(0)$ of radius $r>0$ in $\mathbb{R}^{n_{c}}$, for $\dot{u}(t)=F(u_{t})$.

Then there exists a constant $C>0$ and a radius $r_{\ell}>0$ such $h(\cdot;u_{c})\in C^{\ell}$ and $\|h(\cdot;u_{c}^{0})\|_{\ell}\leq C\|h(\cdot;u_{c}^{0})\|_{0}$ for all $u_{c}\in B_{r_{\ell}}(0)$.

Proof Let $L\geq 0$ be the Lipschitz constant for the right-hand side of the ODE on the center manifold $\dot{u}_{c}=B^{\dagger}\partial_{1}h(\cdot;u_{c})$ on $B_{r}(0)$ (if necesssary, choose $r$ sufficiently small such that $L$ exists). Thus, for all $u_{c}^{0}\in B_{r_{\ell}}(0)$ with $r_{\ell}<r\exp(-\ell\tau_{\max}L)$ the solution of $\dot{u}_{c}=B^{\dagger}\partial_{1}h(\cdot;u_{c})$ starting from $u_{c}(0)=u_{c}^{0}$ does not leave $B_{r}(0)$ for times $t$ with $|t|\leq\ell\tau_{\max}$. Thus, the flow map $U_{c}:[-\ell\tau_{\max},\ell\tau_{\max}]\times B_{r_{\ell}}(0)\ni(t,u_{c}^{0})\mapsto u_{c}(t)\in B_{r}(0)$ is well defined. However, this implies that, for every $u_{c}^{0}\in B_{r_{\ell}}(0)$, $h(\cdot;u_{c}^{0})$ is the solution of the DDE $\dot{u}=F(u_{t})$ starting from $h(\theta;U_{c}(-\ell\tau_{\max};u_{c}^{0}))$. Consequently, by Proposition IV.1, $h(\cdot;u_{c}^{0})$ is in $C^{\ell}$. The relation between the $\|_{\cdot}\|_{\ell}$-norm and the $\|\cdot\|_{0}$-norm follows then from estimate (17) and the Lipschitz constant for $U_{c}(-\ell\tau_{\max};\cdot)$. (End of proof of Proposition IV.2)

Consequently, we can expand at least $F$ in the expression $F(h(u_{c}))$, which is present in the normal form expansion. Humphries et alHumphries, Calleja, and Krauskopf (2016) used this fact to demonstrate for their example how one can expand a sd-DDE near an equilibrium up to order $\ell$ such that all terms of order $j\leq\ell$ are $j$-linear (and have, thus, constant delays). The remainder term is of order $o(\|u_{t}\|_{0}^{\ell})$ and has state-dependent delays. One incurs delays of length up to $\ell\tau_{\max}$ such that we have the following statement, generalizing the approach of Humphries et al:

Lemma IV.3** (Expansion with longer delays)**

Let $F$ be a functional with $C^{\ell}$ coefficients and $m$ discrete state-dependent delays $\tau^{1}$,…,$\tau^{m}$ (of the form (6)–(7)). Let $u_{0}\in C^{1}$ be sufficiently small with $u_{0}^{\prime}(0)=F(u_{0})$. Then the segments $u_{t}$ solving $\dot{u}(t)=F(u_{t})$ satisfy after time $\ell\tau_{\max}$ a sd-DDE of the form

[TABLE]

The $j$-linear functionals $F_{j}$ and the remainder map $C([-\ell\tau_{\max},0];\mathbb{R}^{n})$ into $\mathbb{R}^{n}$. The expansion products $(u_{t})^{j}$ have delays that are sums $\tau^{k_{1}}+\ldots+\tau^{k_{j}}$, where $\{k_{1},\ldots,k_{j}\}\subseteq\{1,\ldots,m\}$ and all delays are evaluated at $u=0$.

Proof

Since after time $t\geq\ell\tau_{\max}$ the solution $u_{t}$ is $\ell$ times continuously differentiable, we can expand the functional $F$ in the equilibrium [math] and in the direction of $u_{t}$ to order $\ell$ using its classical differentiability when restricted to $C^{\ell}$:

[TABLE]

In expansion (19) the $j$-form $\partial^{j}F(0)$ is continuous only on functions in $C^{j-1}$. To keep track of this dependence on the derivatives of $u_{t}$, we include the derivatives explicitly into the multi-linear arguments in (19). To get an expansion that depends on $u_{t}\in C^{0}([-\ell\tau_{\max},0];\mathbb{R}^{n}$ (no derivatives, but longer history), we recursively replace derivatives $u^{(j)}(t)$ by $F^{j}(u_{t})$ (as obtained in Proposition IV.1), followed by expansions of $F^{j}(u_{t})$. A functional $F^{j}:C([-j\tau_{\max},0];\mathbb{R}^{n})\to\mathbb{R}^{n}$ generates also a map $F^{j}_{j+k}$ from $C([-(j+k)\tau_{\max},0];\mathbb{R}^{n})$ into $C([-k\tau_{\max},0];\mathbb{R}^{n}$ for any $k\geq 0$ via $F^{j}_{j+k}(u_{t})(\theta)=F^{j}(u_{t+\theta})$. The subscript $j+k$ indicates the length of the time interval that arguments of $F^{j}_{j+k}$ should have. Thus, after the first replacement of $u^{(\nu)}_{t}$ by $F^{\nu}_{\nu+1}(u_{t})$, we have that for $t\geq\ell\tau_{\max}$, $u$ satisfies

[TABLE]

At subsequent expansions terms from lower orders will change expansions at higher orders. It remains to be shown inductively that eventually all derivatives disappear except for the remainder, and that the length of the history segments $u_{t}$ does never exceed $\ell\tau_{\max}$.

Let us make the inductive assumption that a history segment $u_{t}^{(j)}$ of length $k\tau_{\max}$ shows up at order $(j+1)k\leq\ell$. In the first inductive step we have $k=1$, $j\in\{1,\ldots,\ell-1\}$ and orders at which derivatives of $u_{t}$ appear from $2$ to $\ell$. When replacing $u_{t}^{(j)}$ by $F^{j}_{(j+1)k}(u_{t})$ the history interval increases to $(j+1)k$. Then $F^{j}_{(j+1)k}(u_{t})$ has to be expanded up to order $\left\lceil\ell/((j+1)k)\right\rceil$ ($\lceil r\rceil$ is the lowest integer greater or equal than $r$). In this expansion, we have $\nu$-linear forms containing derivatives of $u_{t}$ up to order $\nu-1$. A derivative of order $i\leq\nu-1$ shows up for orders of $u_{t}$ greater or equal than $(i+1)(j+1)k$.

Hence, a term $u^{(j)}_{t}$ at order $(j+1)k\leq\ell$ creates new $i$th derivative terms ($i\geq 1$) only at order greater or equal than $(i+1)(j+1)k$ such that the recursion must terminate. (We restrict to orders less or equal than $\ell$.) Also, the length of the history interval of the new $i$th derivative term is $(j+1)k\tau_{\max}$, which is less than $\ell\tau_{\max}$, since $(j+1)k\leq\ell$ by inductive assumption.

(End of proof of Lemma IV.3)

We combine the result of Lemma IV.2 with Lemma IV.3 to sharpen the estimate for solutions of the FDE $\dot{u}(t)=F(u_{t})$ starting on the local center manifold: $u_{0}=h(\cdot;u_{c}^{0})$ with $u_{c}^{0}\in B_{r_{\ell}}(0)$. Then the remainder term is also of order $o(\|u_{t}\|_{0}^{\ell})$ (since Lemma IV.2 provides an estimates for $\|h(\cdot;u_{c})\|_{\ell}$ in terms of $\|h(\cdot;u_{c})\|_{0}$:

[TABLE]

Since $u_{t}=h(\cdot;u_{c}(t))$, we may also also replace the remainder by $o(|u_{c}(t)|^{\ell})$. The truncated DDE (20) (dropping the remainder term) has only constant delays. Hence, the semiflow and local center manifold $h_{\mathrm{trunc}}$ of the truncated DDE (20) are smooth, and can, thus, be transformed into normal form with the procedure described in Section III.2. Since this normal form transformation up to order $\ell$ is independent of terms of order $o(|u_{c}|^{\ell})$ and keeps these terms at order $o(|u_{c}|^{\ell})$, we have that for $u$ on the local center manifold $h$ of the non-truncated sd-DDE $\dot{u}(t)=F(u_{t})$, the center component $u_{c}=B^{\dagger}u_{t}$ satisfies an ODE equal to the normal form of the truncated DDE (20) except for a different remainder (still of order $o(|u_{c}|^{\ell})$). The result has the form (compare (10))

[TABLE]

where all coefficients $\alpha_{j}$ are identical to those of the normal form of the truncated DDE (20). However, in contrast to the constant-delay DDE, only the first derivative of the remainder $O(|u_{c}|^{\ell})$ is guaranteed to be small for all small $u_{c}$, but not the higher-order derivatives. This was also demonstrated numerically by Humphries et alHumphries, Calleja, and Krauskopf (2016) for their example. Any phenomenon predicted by the normal form that persists under perturbations of size $o(|u_{c}|^{\ell})$ will also be present in the sd-DDE. This includes all periodic orbits and their changes of stability.

Normally hyperbolic invariant manifolds

For some bifurcations the normal form of the truncated system may predict the presence of, for example, invariant tori that branch off along torus bifurcation curves, away from strong resonances ($1:1$ to $1:4$, seeKuznetsov (2004)). Their degree of normal hyperbolicity is proportional to their distance from the torus bifurcation in the truncated system. Our perturbation (the remainder term $o(|u_{c}|^{\ell})$) is $C^{1}$ small in a ball around [math], but not guaranteed to be $C^{j}$ small compared to lower order terms (with $j>1$), except in [math], because the local center manifold has not been proven to be smooth. Hence, close to the torus bifurcation the invariant tori may be altered by the remainder term. However, the region around the torus bifurcation where the invariant tori are not sufficiently normally hyperbolic shrinks as we approach the neighborhood of [math] if the remainder term decreases faster than the normal hyperbolicity. This is the case if one chooses $\ell$ sufficiently large. For example, Humphries et alHumphries, Calleja, and Krauskopf (2016) indeed reported invariant tori branching off from the torus bifurcation near the Hopf-Hopf interactions as predicted by the normal form. In their paper the authors compared for their example the results from the direct normal form expansion for the sd-DDE as explained in general in Section IV to the results from the constant-delay DDE as constructed via Lemma IV.3 and found agreement up to numerical round-off errors.

V Illustration - position control

A good example suitable for illustration of simple nonlinear behaviour introduced by state-dependence of the delay is the position control problem discussed by Walther Walther (2002) (see also review Hartung et al. (2006)). A mover aims to control its position $x$ relative to an obstacle using linear position feedback (see Figure 1).

We assume that the controlled motion is free of inertia such that (in non-dimensionalized quantities)

[TABLE]

In (22) $k$ is the linear control gain, $x_{0}$ is the reference position that the mover aims to maintain, $x_{\mathrm{est}}$ is the mover’s estimate of the current position, and $\tau_{0}$ is a processing or reaction delay in the control loop. Even if the estimate $x_{\mathrm{est}}(t)$ is perfect (equal to $x(t)$), the equilibrium $x_{0}$ of the controlled system (22) will be linearly unstable if $k\tau_{0}>\pi/2$. If the mover estimates the current position by sending out a signal and measuring the traveling time for the reflected signal then an additional state-dependent delay is introduced. Let $s(t)$ be the time that the reflected signal, arriving at the mover time $t$, needed since leaving the mover, and let $c$ be the signal traveling speed. Then

[TABLE]

The mover estimates its current position via

[TABLE]

Let us introduce the reference travel time $s_{0}=\frac{c}{2}x_{0}$ corresponding to the reference position $x_{0}$. The full equation of motion is

[TABLE]

The differential equation for $s$ follows from (22) and (23) via Baumgarte regularization: we rewrite (23) in the form $g(t)=0$ (where $g(t)=cs(t)-x(t-s(t))-x(t)$), and then replace it by the condition $\frac{\mathrm{d}}{\mathrm{d}t}g(t)=-\gamma g(t)$, re-arranged for $\dot{s}(t)$. Every orbit of (25)–(26) that is periodic or lies on a local center manifold with internal contraction rate less than $\gamma$ satisfies also the algebraic constraint (23). When writing system (25)–(26) in the general form $\dot{u}=F(u_{t})$, the right-hand side of (25)–(26) corresponds to a functional $F$ with the form ($u=(u_{1},u_{2})^{T}=(x,s)^{T}$)

[TABLE]

Equilibria and periodic orbits computed in this illustration had their $s(t)$ component in the range $[s_{\min},s_{\max}]$ with $s_{\min}\geq 0$ and $s_{\max}<10$ in the parameter ranges used for figures 2 and 3. Hence, we may set $\tau_{\max}=10$ and treat $F$ as a functional from $C([-\tau_{\max},0];\mathbb{R}^{2})$ to $\mathbb{R}^{2}$.

For our demonstration we fix $k=1$, $c=2$ and $\gamma=1$ in non-dimensionalized quantities. We vary $\tau_{0}$ and $s_{0}$ in a two-parameter bifurcation study. The system has one constant delay $\tau_{0}$ and two state-dependent delays. In the notation of DDE-Biftool the function $f:\mathbb{R}^{2\times 4}\times\mathbb{R}^{2}\to\mathbb{R}^{2}$ has the time-dependent arguments $u(t-\tau_{j})=[x(t-\tau_{j}),s(t-\tau_{j})]^{T}$ for $j=1,2,3,4$, and the parameters $(\tau_{0},s_{0})$, where

[TABLE]

The system (25)–(26) has a unique equilibrium at $u_{*}=(x_{*},s_{*})=(cs_{0}/2,s_{0})$. As part of the principle of linearized stability proved by WaltherWalther (2003) comes the description for how to compute stability (which is implemented in DDE-Biftool): “freeze” the state-dependent delays at the values in the equilibrium, and then compute the linearization of the corresponding DDE with constant delays Hartung and Turi (2001); Cooke and Huang (1996); Hartung et al. (2006). For the position control problem this procedure gives a algebraic relation between the parameter values at which Hopf bifurcations occur:

[TABLE]

The Hopf bifurcation that forms the boundary of the stability region in the $(\tau_{0},s_{0})$-plane is the curve for $\omega_{0}^{+}$, shown in Figure 2 (right panel) as a green dashed/solid curve.

As expression (27) is still implicit, the curve in Figure 2 was computed with DDE-Biftool. The standard Hopf bifurcation theorem can be applied to sd-DDEsEichmann (2006); Sieber (2012) such as system (25)–(26). Hence, a family of periodic orbits branches off from the Hopf bifurcation. Near the equilibrium the stability of periodic orbits can be predicted using the expression (11) for $L_{1}$ as implemented by Kuznetsov et alWage (2014); Bosschaert (2016); Bosschaert, Janssens, and Kuznetsov (2017). This was rigorously proven using a Lyapunov-Schmidt reduction for periodic boundary value problemsSieber (2012). Its value along the Hopf curve is shown in the left panel of Figure 2. The value of $L_{1}$ crosses zero at $s_{0}\approx 4.02$, $\tau=1.05$. There the Hopf bifurcation is degenerate and the second Lyapunov coefficient is $L_{2}\approx-1.9\times 10^{-3}$. This implies that the family of periodic orbits exists to the right and is stable where the Hopf curve is solid in Figure 2. The family of periodic orbits is unstable and exists to the left, before folding in a fold of periodic orbits to the right where the Hopf curve is dashed in Figure 2.

VI Conclusion

As this paper shows, expressions for normal form coefficients for constant-delay DDEs can be generalized to sd-DDEs. The mathematical justification is only partially complete, but for many phenomena it is already clear how they persist when the truncation is removed. The complete justification requires smoothness for the local center manifold. Krisztin has provisional resultsKrisztin (2006a) that show how his proof for smooth unstable manifolds of equilibriaKrisztin (2003) can be extended to local center manifolds. Ideally, the general result for persistence of compact normally hyperbolic manifolds should in some sense be adapted to sd-DDEs in the following form. Consider a sd-DDE of the form

[TABLE]

where $F_{\mathrm{c}}:C^{0}\to\mathbb{R}^{n}$ is smooth and $\dot{u}(t)=F_{\mathrm{c}}(u_{t})$ has a compact overflowing invariant normally hyperbolic (say, stable) manifold ${\cal M}_{0}$. If we also assume that $F_{\mathrm{sd}}$ has a sufficiently small Lipschitz constant with respect to the space of Lipschitz continuous functions $C^{0,1}$ (and is mildly differentiable up to order $\ell$), then (28) should also have a compact overflowing invariant normally stable manifold ${\cal M}$. The smoothness of ${\cal M}$ should only be restricted by the spectral gap in the exponential dichotomy on ${\cal M}_{0}$.

Acknowledgements.

J.S. gratefully acknowledges the financial support of the EPSRC via grants EP/N023544/1 and EP/N014391/1. J.S. has also received funding from the European Union’s Horizon 2020 research and innovation programme under Grant Agreement number 643073.

Appendix A Details of normal form expansion for local bifurcations of DDEs

This appendix gives a few additional details for the computation of coefficients in the normal form procedure of Section III.

The linear DDE $\dot{u}=Au_{t}$

Recall that the characteristic matrix is denoted by $\Delta(\lambda)\in\mathbb{C}^{n\times n}$, which has $n_{c}$ eigenvalues on the imaginary axis (counting multiplicity). Let $B=\{b_{1},\ldots,b_{n_{c}}\}$ be a basis of the linear center subspace $U_{c}=\operatorname{span}B$ of $\dot{u}(t)=Au_{t}$. A spectral projection $P_{c}$ onto the space $U_{c}$ is given by residue of the resolvent $R(\lambda)$:

[TABLE]

where the curve integral is taken around the critical spectrum $\sigma_{c}$. The resolvent $R(\lambda)$, mapping $C^{0}$ into $C^{1}$ is defined as the unique solution $x\in C^{1}$ of

[TABLE]

which is

[TABLE]

We define $B^{\dagger}:C^{0}\ni x\mapsto x_{c}\in\mathbb{R}^{n_{c}}$, where $x_{c}\in\mathbb{R}^{n_{c}}$ is the unique vector of coordinates such that $Bx_{c}=P_{c}x$. Thus, $B^{\dagger}B$ is the identity in $\mathbb{R}^{n_{c}}$, and $BB^{\dagger}=P_{c}$.

Center manifold expansion

The semiflow of the DDE, restricted to the center manifold $\{u\in C^{0}:u(\theta)=h(\theta;u_{c}),u_{c}\in\mathbb{R}^{n_{c}}\mbox{\ small}\}$, introduced in Section III, satisfies the ODE in $\mathbb{R}^{n_{c}}$

[TABLE]

The invariance of graph of the manifold

[TABLE]

under the DDE $\dot{u}=F(u_{t})$ implies

[TABLE]

Let us introduce expansions for $F$ and $h(\theta;\cdot)$ up to order $\ell$ in the point $u=0$ (for $F$) and $u_{c}=0$ (for $h(\theta,\cdot)$):

[TABLE]

The first-order coefficient $F_{1}$ of $F$ is the linear operator $A$, the first-order coefficient $h_{1}(\theta)$ of the manifold graph is $B(\theta)$. The coefficients $h_{j}$ for $j>1$ are only determined up to conjugacy of the flow on the center manifold to order $j$. A different choice of $h_{j}$ corresponds to a different, but conjugate, ODE for $u_{c}$. For example, requiring $B^{\dagger}h_{j}[u_{c}]^{j}=0$ for all $j>1$ and all $u_{c}\in\mathbb{R}^{n_{c}}$ would determine $h_{j}$ uniquely in combination with the invariance (32)–(33).

Determining systems for coefficients $h_{j}(0)$ and $\alpha_{j}$

However, the approach proposed by KuznetsovKuznetsov (1999) and taken in DDE-Biftool’s normal form extensionJanssens (2010); Wage (2014); Bosschaert (2016); Bosschaert, Janssens, and Kuznetsov (2017) is to choose the expansion coefficients $h_{j}$ such that the ODE (31) on the center manifold for $u_{c}$ is already in normal form:

[TABLE]

In (34) the matrix $A_{c}^{1}=B^{\dagger}\circ[\partial/\partial\theta]\circ B=B^{\dagger}\circ B^{\prime}\in\mathbb{R}^{n_{c}\times n_{c}}$ is the projection of the linear DDE on the eigenspace for the spectrum $\sigma_{c}$ on the imaginary axis. For higher orders $j>1$ the coefficients $A_{c}^{j}$ are given except for a finite number of to-be-determined normal form coefficients $\alpha_{j}$. We use square brackets to indicate that $A_{c}^{j}$ is a given map depending linearly on $\alpha_{j}$ and $j$-linearly on $u_{c}$. The coefficient $\alpha_{j}$ may be empty (for example, $\alpha_{1}$ is always empty). Inserting the expansions for $h$, $F$ and $\dot{u}_{c}$ into the invariance equation (33) gives at order $j$ a $n$-dimensional inhomogeneous constant-coefficient differential equation for each coefficient of the symmetric $j$-form $h_{j}(\theta)$:

[TABLE]

where

[TABLE]

is a known function determined by orders lower than $j$ (it is not present for orders $1$ and $2$. Let us denote the solution $h_{j}$ of the affine ordinary differential equation (35) by

[TABLE]

The above expression indicates that the solution is linear in $h_{j}^{0}=h_{j}(0)$ (its initial value), $\alpha_{j}$ and $R_{j}$, and $j$-linear in $u_{c}$. If the basis $B$ consists only of eigenvectors (eigenvector $b_{i}$ for eigenvalue $\lambda_{i}$), then $A_{c}^{1}$ is diagonal, and $H_{j}(\theta)=\exp(\lambda_{i}\theta)h_{j,\nu}^{0}$ for coefficients $h_{j,\nu}$ of the $j$-form $h_{j}(\theta)$. In this case the ${n+j\choose j}$ differential equations for the ${n+j\choose j}$ coefficients $h_{j,\nu}$ of the $j$-form $h_{j}(\theta)$ decouple. The initial conditions $h_{j}^{0}$ are determined by the invariance at $\theta=0$, (32):

[TABLE]

The second sum is taken over multi-indices $\nu\in\operatorname{ind}(k,j)$. The set $\operatorname{ind}(j,k)$ is the set of $k$-tuples of positive integers summing up to $j$. Inserting the differential equation for $h_{j}$ and its solution $H_{j}$ at $\theta=0$ results in an affine equation for $h_{j}^{0}$ and $\alpha_{j}$ (the homological equation):

[TABLE]

One can determine $h_{j}^{0}$ and $\alpha_{j}$ for each $j$ by comparing coefficients of this $j$-form in $u_{c}$. For orders $j$, for which the square coefficient matrix $L_{h,j}$ is regular, the normal form coefficient $\alpha_{j}$ is not present (since all terms at this order are non-resonant). If the matrix $L_{h,j}$ is singular with kernel dimension $d_{j}$, then the dimension of $\alpha_{j}$ is $d_{j}$ and the dependence of $A_{c}^{j}$ on $\alpha_{j}$ is such that $[L_{h,j},-L_{\alpha,j}]$ has full rank. Thus, there is a unique coefficient $\alpha_{j}$, for which (36) is solvable for $h_{0}^{j}$. The solution $h_{0}^{j}$ is not unique, but can be made unique, for example, by forcing it to be orthogonal to the nullspace of $L_{h,j}^{T}$; see the referencesJanssens (2010); Wage (2014); Bosschaert (2016); Bosschaert, Janssens, and Kuznetsov (2017).