{'fraction_non_alphanumeric': 0.047849929072007795, 'fraction_numerical': 0.022845165251635578, 'mean_word_length': 4.090698426384996, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Switching to nonhyperbolic cycles from codimension two bifurcations of equilibria of delay differential equations

In this paper we perform the parameter-dependent center manifold reduction near the generalized Hopf (Bautin), fold-Hopf, Hopf-Hopf and transcritical-Hopf bifurcations in delay differential equations (DDEs). This allows us to initialize the continuation of codimension one equilibria and cycle bifurcations emanating from these codimension two bifurcation points. The normal form coefficients are derived in the functional analytic perturbation framework for dual semigroups (sun-star calculus) using a normalization technique based on the Fredholm alternative. The obtained expressions give explicit formulas which have been implemented in the freely available numerical software package DDE-BifTool. While our theoretical results are proven to apply more generally, the software implementation and examples focus on DDEs with finitely many discrete delays. Together with the continuation capabilities of DDE-BifTool, this provides a powerful tool to study the dynamics near equilibria of such DDEs. The effectiveness is demonstrated on various models.

Using the framework of perturbation theory for dual semigroups developed in [4,5,6,7] the existence of a finite dimensional smooth center manifold for DDEs can be rigorously established [14]. As a consequence the normalization method for local bifurcations of ODEs developed in [42] can be lifted [37] rather easily to the infinite dimensional setting of DDEs. One of the advantages of this normalization technique is that the center manifold reduction and the calculation of the normal form coefficients are performed simultaneously by solving the so-called homological equation. The method gives explicit expressions for the coefficients rather than a procedure as developed in [22,23]. The critical normal form coefficients for all five generic codimension two bifurcations of equilibria of DDEs have been derived [37] and partially implemented [60] into the fully GNU Octave compatible MATLAB package DDE-BifTool [21,56].
In this paper we will perform the parameter-dependent center manifold reduction and normalization for three codimension two Hopf cases: the generalized Hopf, fold-Hopf and Hopf-Hopf bifurcations. This will allow us to initialize the continuation of codimension one bifurcation curves of nonhyperbolic equilibria and cycles emanating from the codimension two points. These are the only codimension two bifurcation points of equilibria in generic DDEs where codimension one bifurcation curves of nonhyperbolic cycles could originate. We also treat the more special transcritical-Hopf bifurcation which is frequently found in applications.
The center manifold theorem for parameter-dependent DDEs as presented in [14] assumes explicitly that the steady state exists for all nearby parameter values. However, for a generic fold-Hopf bifurcation this assumption is not satisfied. An attempt to deal with this complication has been made in [30], where it is discussed how to reduce a parameter-dependent DDE to a DDE without parameters by appending the trivial equationα = 0. However, the reduction in [30] is based on the formal adjoint approach [32] and applies specifically to DDEs, while at times it lacks consistency. Therefore we demonstrate in this paper how the reduction to the parameter-independent case can be done in the sun-star framework, enabling a rigorous approach to the existence of parameter-dependent center manifolds for a class of evolution equations that includes DDEs. This allows us to treat bifurcations of equilibria with zero eigenvalues in generic DDEs while at the same time achieving applicability of our results to other classes of delay equations. This paper is organized as follows. In Section 2 we offer a concise review of perturbation theory for dual semigroups (also called sun-star calculus), both on an abstract level as well as in application to the analysis of classical DDEs as dynamical systems. We also recall from [37] various results that are needed for the normalization.
In Section 3 we show how the theory from the previous section also applies to parameter-dependent classical DDEs by converting them into a parameter-independent system on a product state space. We again present the material in two stages: Results are first established at a more abstract semigroup level and next applied to classical DDEs depending on parameters. In particular, we define the parameterdependent local center manifold and give an explicit ODE for solutions that are confined to it.
In Section 4 we describe the general technique used to derive expressions for the normal form coefficients on the parameter-dependent center manifold in the infinite dimensional setting of classical DDEs.
Before we apply this technique to the previously mentioned codimension two bifurcations, we summarize in Section 5 relevant smooth normal forms and we list asymptotics for the codimension one cycle bifurcation curves emanating from the codimension two points as well as for the corresponding nonhyperbolic equilibria.
In Section 6 the method is then applied to the generalized Hopf (Bautin), fold-Hopf, Hopf-Hopf and transcritical-Hopf bifurcations in classical DDEs. We provide explicit expressions for all normal form coefficients necessary for the predictors of codimension one bifurcation curves. While most of the critical normal form coefficients for these bifurcations were obtained in [37], we briefly re-derive them to ensure readability.
In Section 7 we provide explicit computational formulas for the evaluation of the linear and multilinear forms used in the normal form coefficients and predictors for the simplest subclass consisting of DDEs (1) with finite many discrete delays. These formulas are actually implemented in DDE-BifTool.
In Section 8 we employ our implementation in DDE-BifTool to illustrate the accuracy of the codimension one bifurcation curve predictors through various example models displaying all aforementioned degenerate Hopf cases. A complete step-by-step walk-through of the examples, including all code to reproduce the obtained results, is provided in the Supplement.

Dual perturbation theory and classical DDEs
We begin by presenting those general elements of perturbation theory for dual semigroups that are useful for the study of classical DDEs as dynamical systems. Throughout we assume sun-reflexivity -a term that will be introduced in Section 2.1. From Section 2.4 onward, we then explain how the general results apply to classical DDEs. The standard reference for this entire section is [14], while for the underlying theory of semigroups of linear operators we recommend [19,20].

Duality structure and linear perturbation
The starting point is a C 0 -semigroup T 0 on a real or complex Banach space X. Let A 0 with domain D(A 0 ) be the infinitesimal generator (or: generator, for short) of T 0 . We denote by X the topological dual space of X and we use the prefix notation for the pairing between x ∈ X and x ∈ X, x , x := x (x).
If X is not reflexive then the adjoint semigroup T 0 is in general only weak continuous on X and A 0 generates T 0 only in the weak sense. The maximal subspace of strong continuity X := {x ∈ X : t → T 0 (t)x is norm-continuous on R + } is invariant under T 0 and we have the characterization where the bar denotes the norm closure in X . By construction the restriction of T 0 to X is a C 0 -semigroup that we denote by T 0 . Its generator A 0 is the part of A 0 in X , At this stage we again have a C 0 -semigroup T 0 with generator A 0 on a Banach space X so we can iterate the above construction. On the dual space X we obtain the adjoint semigroup T 0 with weak generator A 0 . By restriction to the maximal subspace of strong continuity X = D(A 0 ) we end up with the C 0 -semigroup T 0 . Its generator A 0 is the part of A 0 in X .
The canonical injection j : X → X defined by maps X into X . If j maps X onto X then X is called -reflexive (pronounce: sun-reflexive) with respect to T 0 . One may define an equivalent norm on X with respect to which j becomes an isometry, but this need not be assumed. However, sun-reflexivity of X with respect to T 0 will be assumed throughout.
With the abstract duality structure in place, we next turn our attention to perturbation. Let L : X → X be a bounded linear operator. Then there exists a unique C 0 -semigroup T on X that satisfies the linear integral equation where the weak Riemann integral takes values in X and the running assumption of sun-reflexivity justifies the application of j −1 . By using (3) to express the difference T − T 0 of the perturbed and the unperturbed semigroups, one proves that the maximal subspaces of strong continuity X and X are the same for T and T 0 , so there is no need to distinguish them with a subscript. In particular, X is sun-reflexive also with respect to T . On X the perturbation L appears additively in the action of A , We recover the generator A of T by considering the part of A in X . As a consequence L moves into the domain and we find D(A) = x ∈ X : jx ∈ D(A 0 ) and A 0 jx + Lx ∈ X , Ax = j −1 (A 0 jx + Lx).
For proofs of the statements so far, see [14, Appendix II.3 and Chapter III].

Nonlinear perturbation and linearization
The C 0 -semigroup T arose as a linear perturbation of the original C 0 -semigroup T 0 , so the next step is to introduce a nonlinear perturbation of T itself. In keeping with the tradition for nonlinear problems [14, Sections VII.1 and VIII.1] we only regard the case that X is a real Banach space, also see Remark 1 below. Let R : X → X be a C k -operator for some k ≥ 1 such that and consider the nonlinear integral equation Due to the nonlinearity, for a given initial condition x ∈ X one can at most guarantee existence of a maximal solution u x : I x → X of (5) on a forward time interval I x := [0, t x ) for some 0 < t x ≤ ∞ [14, Chapter VII]. The family of all such maximal solutions defines a nonlinear semiflow Σ : D(Σ) → X, that may in addition depend on parameters [14, Defs. VII.2.1 and VII. 2.9]. (For reasons discussed in Section 3, we will treat parameter dependence differently and separately. Until then, the reader can consider all parameters to be held fixed and absent in the notation.) The domain of Σ is open in [0, ∞) × X and 0 ∈ X is a stationary point of Σ, The semiflow Σ is (in fact, uniformly) differentiable with respect to the state at (t, 0) ∈ D(Σ), with the partial derivative where T is the C 0 -semigroup that satisfies (3).
Remark 1. For nonlinear problems it is customary to work on a real Banach space X. The reason is that these problems often come from concrete equations with nonlinear right-hand sides for which it is unclear if and how they can be extended to complex arguments. Consequently, if we want to analyze the linearization of Σ at 0 ∈ X using spectral theory, then it becomes necessary to complexify X and the linear operators acting on X [14, Section III.7 and last part of Section IV.2], [52,Section 1.3]. In particular, by the spectrum of A we mean the spectrum of its complexification on the complexified Banach space.

Critical local center manifolds
As in Section 2.2 we continue to assume that T 0 is a C 0 -semigroup on a real Banach space X that is sun-reflexive with respect to T 0 . In addition we assume that T 0 is eventually compact and L is a compact operator. This implies that the perturbed semigroup T defined by (3) is eventually compact as well [11,Theorem 2.8].
When considering solutions that exist for all (positive and negative) time -such as periodic orbits -it is useful to write (5) in the translation invariant form A solution of (8) is a continuous function u : I → X on some nondegenerate -possibly unboundedinterval I ⊆ R that satisfies (8) for all s, t ∈ I with s ≤ t. Naturally, u is a solution of (8) if and only if t − s ∈ I u(s) , u(t) = Σ(t − s, u(s)), ∀ s, t ∈ I with s ≤ t, where Σ : D(Σ) → X is the nonlinear semiflow from (6). The interval I is often left implicit. The general center manifold theorems from [14, Chapter IX] for equations of the type (8) apply to the particular case where T is an eventually compact C 0 -semigroup on a real, sun-reflexive Banach space. Let us therefore suppose that 0 ∈ X is a nonhyperbolic equilibrium of Σ, so the generator A of T possesses 1 ≤ n 0 < ∞ purely imaginary eigenvalues, counting algebraic multiplicities -see Remark 1. Let X 0 ⊆ X be the real center eigenspace corresponding to these eigenvalues. Then there exists a C k -smooth n 0 -dimensional local center manifold W c loc that is tangent to X 0 at the origin. Any solution u : I → X of (8) that lies on W c loc is differentiable on I and satisfies where A is the weak generator of T . We note that (9) is an identity in X .

The special case of classical DDEs
It will now be explained how the general results from Sections 2.1 to 2.3 apply to classical DDEs. We choose the nonreflexive Banach space X := C([−h, 0], R n ) as the state space, introduce a C k -smooth operator F : X → R n , and consider an equation with a finite delay 0 < h < ∞ of the forṁ with an initial condition For each t ≥ 0, the function x t : [−h, 0] → R n defined by is called the history of the unknown function x at time t. Equations of the type (DDE) will be called classical DDEs. Note that (1) is quite literally a case in point. By a solution of the initial value problem (DDE)-(IC) we mean a continuous function x : [−h, t + ) → R n for some 0 < t + ≤ ∞ that is differentiable on [0, t + ) and satisfies (DDE) and (IC). When t + = ∞ we call x a global solution.
We want to study (DDE) near an equilibrium at the origin, so assume that F (0) = 0 and split F into its linear and nonlinear parts, Here ζ : [0, h] → R n×n is a matrix-valued function of bounded variation, normalized by the requirement that ζ(0) = 0 and ζ is right-continuous on the open interval (0, h). The integral is of the Riemann-Stieltjes type, and G : X → R n is a C k -smooth nonlinear operator. It is common to denote the linear part more succinctly as so that F (ϕ) = ζ, ϕ + G(ϕ), ϕ ∈ X.
We first consider the case G = 0, whence (DDE) reduces to the linear equatioṅ In order to understand the relationship between (12) and (3) we begin by observing that the trivial with initial condition (IC) has the obvious solution Using this solution, we define the strongly continuous shift semigroup T 0 on X by We note that T 0 (h) is a compact operator, so T 0 is eventually compact. For this particular combination of X and T 0 the abstract duality structure from Section 2.1 can be constructed systematically and explicitly [14,Section II.5]. We only summarize the few facts that will be used in the sequel.
Remark 2 (Notation). For K ∈ {R, C} let K n be the linear space of column vectors and let K n be the linear space of row vectors, both over K. Elements of K n are denoted by q = (q 1 , q 2 , . . . , q n ) -commas between the entries -while elements in K n are denoted by p = (p 1 p 2 · · · p n ) -no commas between the entries. We sometimes use the pairing defined by the row-column matrix multiplication: Note that the standard Hermitian inner product between two vectors p T , q ∈ C n should be written as p · q and not as p · q.
On X : The maximal domain of strong continuity of T 0 has the representation and the duality pairing between ϕ = (c, g) ∈ X and ϕ ∈ X is On X : Switching to the dual space of (15) yields the representation and the duality pairing between ϕ = (a, ψ) ∈ X and ϕ = (c, g) ∈ X is The canonical injection (2) sends ϕ ∈ X to jϕ = (ϕ(0), ϕ), mapping X onto X . Therefore X is sun-reflexive with respect to the shift semigroup T 0 .
Next, we specify the linear and nonlinear perturbations L and R in (3) and (5), respectively, and to relate these two abstract integral equations in X to the linear and nonlinear initial value problems for (DDE). For i = 1, . . . , n we denote r i := (e i , 0) where e i is the ith standard basis vector of R n . It is conventional and convenient to introduce the shorthand and we note that wr = (w, 0) ∈ X . First we define the compact linear perturbation in (3) as where the pairing in the right-hand side is given by (10). Now (12) with (IC) is equivalent to (3) with (18) in the following sense: If T is the unique C 0 -semigroup on X satisfying (3) with (18) then is the unique global solution of (12) with (IC) and It remains to specify the nonlinear perturbation R in (5) as where G is the nonlinear operator appearing in the splitting (11). Let Σ as in (6) be the nonlinear semiflow generated by the family of maximal solutions of (5) with (19). The equivalence between (DDE)-(IC) and (5) with (19) can be formulated as follows [14, Prop. VII.6.1]. The function x ϕ : is the maximal solution of (DDE)-(IC), in the sense that any other solution necessarily exists only on a subinterval [−h, t + ) for some 0 < t + ≤ t ϕ and coincides with x ϕ there. Moreover, It is the content of (7) that generation and linearization commute: Starting with (DDE), linearization of the semiflow Σ at the stationary point 0 ∈ X yields precisely the eventually compact C 0 -semigroup T corresponding to the linearized DDE (12).

Spectral computations for classical DDEs
The eventual compactness of T implies that the spectrum of its generator A -see Remark 1 -consists entirely of isolated eigenvalues of finite algebraic multiplicity. It is clear from (18) that L is not just compact, but actually of finite rank. This implies that all spectral information about A is contained in a holomorphic characteristic matrix function ∆ : C → C n×n defined by where ζ is the real kernel from (10) [14,Sections IV.4 and IV.5]. In particular, the eigenvalues of A are the roots of the characteristic equation and the algebraic multiplicity of an eigenvalue equals its order as a root of (21). We will be concerned exclusively with simple eigenvalues, for which the geometric and algebraic multiplicities are both equal to one. Let λ ∈ C be such a simple eigenvalue of A. There exist nonzero right and left null vectors q ∈ C n and p ∈ C n of ∆(λ), The second equation is of course equivalent to p T being a nonzero right null vector of ∆ T (λ). The one-dimensional eigenspaces of A and A corresponding to λ are spanned by eigenfunctions ϕ and ϕ , respectively, with and We note that we have implicitly used -and will use consistently -the complexifications of X and of the representation (15) of X . For a simple eigenvalue λ, where the duality pairing is understood to be the complexification of (16). This nonequality implies that the eigenfunctions can be normalized to satisfy ϕ , ϕ = 1. In fact, from (16) and (22) one so this normalization can be effectuated by scaling p and q such that p∆ (λ)q = 1. Finally, it is easily seen that if µ = λ is another simple eigenvalue of A with eigenvector ψ and adjoint eigenvector ψ , then ϕ , ψ = 0, ψ , ϕ = 0.

Solvability of linear operator equations
When computing the normal form coefficients in Section 6 using (HOM), we will encounter linear operator equations of the form where z is a complex number, (w 0 , w) ∈ X is given and (v 0 , v) ∈ D(A ) is the unknown. In general, both z and the right-hand side will have a nontrivial imaginary part, so here and from here onward, it is necessary to regard systems of the form (26) as the complexification of the original operator equations.
We will however not attach additional subscripts to the operator symbols, hoping that this omission will not cause confusion. Since σ(A) consists exclusively of point spectrum, there are two situations to consider depending on whether or not z is an eigenvalue. If z is not an eigenvalue of A then z belongs to the resolvent set ρ(A) of A and (26) admits a unique solution, In order to actually find this solution, one needs a representation of the resolvent operator of A . The general result can be found in [14,Corollary IV.5.4], but here we only require a special case.
Lemma 3. Suppose that z is not an eigenvalue of A, so (26) has a unique solution (v 0 , v). If the right-hand side is represented by for some fixed vector η ∈ C n , then this solution has the representation , use the linearity of (zI − A ) −1 and apply both cases of [37,Corollary 3.4].
On the other hand, suppose that z = λ is an eigenvalue. Then (26) need not be consistent. In fact, a solution exists if and only if see [37,Lemma 3.2]). This condition is often referred to as the Fredholm solvability condition. We note that the duality pairing in (FSC) may be evaluated in concrete cases using (22) and the complexification of (17). This will be done many times in Section 6 when we apply (FSC) to specific operator equations. If z = λ is an eigenvalue and (26) is consistent, then clearly its solutions are not unique. The bordered operator inverse is used to select a particular solution in a systematic and convenient way. For the case that λ is a simple eigenvalue, it assigns the unique solution of the extended linear system to every (w 0 , w) for which (26) is consistent. The following lemma gives an explicit representation for a special case [37, Proposition 3.6 and Corollary 3.7].
In Section 6 we will use the shorthand notation for the solution in Lemma 4. We observe that the expression for ξ itself involves a bordered matrix inverse, which assigns the unique solution of the extended linear system to every y ∈ C n for which the system ∆(λ)x = y is consistent -also see Remark 2 for the notation. In practice, x = ∆ INV (λ)y can be obtained by solving the nonsingular bordered matrix system for the unknown (x, s) ∈ C n+1 that necessarily satisfies s = 0. The properties of (finite dimensional) bordered linear systems and their role in numerical bifurcation analysis are discussed more extensively in [41] and [28,Chapter 3].

Parameter dependence and classical DDEs
In Section 3.1 we motivate our approach by explaining why the standard literature result does not apply to the problem at hand. This is most easily done at the concrete level of classical DDEs. The structure of the remaining subsections parallels that of Section 2. Namely, we first solve the problem of parameter dependence at the more abstract level of dual perturbation theory. In the final Section 3.6 we then return to classical DDEs to see how the general results apply in this special case.

Motivation
We are concerned with the situation where the right-hand side of (DDE) depends explicitly on parameters. Specifically, we considerẋ where F : X × R p → R n is C k -smooth for some k ≥ 1 with F (0, 0) = 0. We assume that at the critical parameter value α = 0 the linearization of (28) has 1 ≤ n 0 < ∞ purely imaginary eigenvalues, counting multiplicities. The goal of Section 3 is to obtain a parameter-dependent family of local center manifolds for a class of evolution equations that includes (28). In [14, Section IX.9.1] this problem is approached as follows. One augments (28) with a trivial equation for the constant parameter dynamics. This gives the system on the state space X := X × R p , with X = C([−h, 0], R n ) as before in Section 2.4. Then the right-hand side of the first equation of (29) is split as which definesG : X → R n , cf. [14, (9.7) in Section IX.9.1]. The first term on the right of (30) acts only on the X-component of the state in X, so the semigroupT on X obtained by perturbing the shift-semigroup T 0 is diagonal.
However, there is an obstruction. In order to satisfy the hypotheses of the parameter-independent center manifold theorem,G must be a pure nonlinearity on X, i.e.
The first two of these conditions are clearly fulfilled, but in general there is no reason for the third condition to be met. It does hold whenG(0, α) = 0 for all α ∈ R p in a neighborhood of zero, i.e. when the zero equilibrium of (28) persists under small parameter variations. For a generic fold-Hopf bifurcation -as well as for a generic Bogdanov-Takens bifurcation that we do not discuss here -there is no such persistence. In this article, the above difficulty is resolved by considering instead of (30) the splitting Using this splitting, (29) is written as Now both D 1 F (0, 0) as well as D 2 F (0, 0) appear in the perturbation of T 0 . As a consequence the perturbed semigroup T is no longer diagonal, but still simple enough for a complete analysis. Moreover, so the parameter-independent center manifold theorem can be applied without having to assume equilibrium persistence. Of course G =G and T =T whenever D 2 F (0, 0) = 0. Remark 5. In a first attempt we regarded the augmented system (29) as a DDE on the state space , also see [30], but we found this approach to be a bit unsatisfactory: The proofs in Sections 3.2 to 3.5 do not depend on the details of the class of equations under consideration, while those same details sometimes lead to notation that is more complicated than necessary.

Duality structure and linear perturbation
We work in the setting of Section 2.1. Namely, let T 0 be a C 0 -semigroup on a real or complex Banach space X that is sun-reflexive with respect to T 0 . We write K ∈ {R, C} for the underlying scalar fieldas in Remark 2. Define T 0 on X by The procedure of taking adjoints and restrictions (twice) then yields semigroups T 0 , T 0 , T 0 and T 0 on X X × K p , X X × K p , X X × K p and X X × K p . (The symbol indicates an identification via a natural isometric isomorphism.) It is straightforward to check that on and that the canonical injection j : X → X has the form where j : X → X is the canonical injection from (2). In particular, X is sun-reflexive with respect to T 0 .
As in Section 2.1 we now introduce a bounded linear perturbation L : X → X of T 0 . We let it be of the form with L : X → X and L p : K p → X bounded linear operators. Perturbing T 0 by L and T 0 by L yields C 0 -semigroups T on X and T on X, respectively. Let A and A be their generators.
Remark 6. There are at least two equivalent ways to compute T and A on X and their weak counterparts T and A on X . One approach -suggested to us by Odo Diekmann -uses integrated semigroup theory to calculate first T and next T . Then A and A are calculated, in that order.
Here we go the other way around: We start by calculating A and use it to obtain T . If desired A and T can then be found by application of (36) and its inverse. This approach is more elementary -we use only theory that was already introduced in Section 2.1 -and it yields the same outcome, as it should.
Proposition 7. The weak generator A of T has the representation where A is the weak generator of T .
Proof. According to the general theory of Section 2.1 and (4) in particular, we have Using (36) and (37) we calculate and the result follows.
Lemma 8. Let L p : X → K p be the restriction of L p to X . Then L p = L p .
Proof. We begin by noting that -strictly speaking -this statement involves two canonical identifications. Namely, let i : X → X and i p : K p → K p be the canonical injection and bijection, respectively. Then L p := L p i and we need to prove that L p i p = L p . For this it is not difficult to show that L p i p α, ϕ = L p α, ϕ , for all α ∈ K p and for all ϕ ∈ X .
For the purpose of notation, we define the integrated semigroup W for T as with on the right a weak Riemann integral of the same type as the integral in (3).
Proposition 9. The semigroup T that is weakly generated by A has the representation Proof. We define a one-parameter family S of bounded linear operators on X by It is easy to check that S is a C 0 -semigroup on X . By Lemma 8 the adjoint semigroup S (t) equals the right-hand side of (39) for all t ≥ 0. We will show that the weak generator A of T is also the weak generator of S . This will then imply that T = S . We use Proposition 7. Let C be the generator of S, so C is the weak generator of S . For any (ϕ , α) in X and any t > 0 we have We note that t −1 W (t)L p α → L p α weakly as t ↓ 0. It follows that the right-hand side converges weakly if and only if ϕ ∈ D(A ) and in that case the weak -limit equals (A ϕ + L p α, 0) = A (ϕ , α). We conclude that C = A .

Spectral theory and the center eigenspace
Let T 0 be a C 0 -semigroup on a complex Banach space X that is sun-reflexive with respect to T 0 . For the purpose of spectral theory, we explicitly take C as the underlying scalar field. In examples, X will often be a complexification of a real Banach space, see Remark 1.
We are interested in a description of the spectrum and the corresponding (generalized) eigenspaces of the generator A of T . In particular, Propositions 11 and 14 below guarantee, respectively, the existence and smooth parametrization of the parameter-dependent local center manifold in Section 3.5.
for every z in the resolvent set ρ(A ).
Proof. From Proposition 7 we have This upper triangular operator matrix has a bounded inverse if and only if both entries on its diagonal have bounded inverses, which happens if and only if z ∈ ρ(A ) and z = 0. In that case, the inverse is given precisely by the stated expression for R z (A ) := (zI − A ) −1 .
In addition we assume that T 0 is eventually compact and the perturbation L in (37) is compact. As a consequence, the spectral analysis of A reduces to an analysis of the poles of its resolvent operator [19,Corollary V.3.2], [58, Section V.10].
Proposition 11. T is an eventually compact C 0 -semigroup.
Proof. The eventual compactness of T 0 , the finite rank of I p and (34) together imply that T 0 is eventually compact. Since L p has finite rank and L is compact by assumption, it follows from (37) that L is compact, so T is eventually compact by [11,Theorem 2.8].
Theorem 12. The generalized eigenspace corresponding to λ ∈ σ(A ) is given by where Γ 0 is a bounded linear operator on X mapping into X .
Proof. Let λ ∈ σ(A ) be arbitrary. Taking residues at z = λ in (40), we obtain where δ λ := δ λ,0 is the Kronecker delta and P λ and P λ are the spectral projectors corresponding to λ for A and A . (If λ is in the resolvent set of the respective operator, then the residue -hence the spectral projector -is identically zero.) Γ λ is pointwise equal to a contour integral with an integrand in the closed subspace X of X , so Γ λ maps into X . We will now calculate the range of P λ from (41). In general, First we assume that λ = 0, so δ λ = 0. We are going to show that Together with (42) this will then prove the theorem for M λ (A ). To verify (43) let p ∈ N be the order of λ as a pole of R z (A ). For n = 1, . . . , p let B n be the coefficient of (z − λ) −n in the Laurent series for R z (A ). A small computation shows that From [58, Section V.10] we recall the relation B n+1 = (A − λI) n B 1 for all n ∈ N. Since B 1 = P λ and its range M λ (A ) is an invariant subspace of A , this relation implies that B k takes values in M λ (A ) for all k = 1, . . . , p, so the same is true for Γ λ by (44). From this it follows that (43) holds. For the remaining case λ = 0 we have δ λ = 1, so from (42) we get the direct sum The first summand equals M 0 (A ) × {0} and this gives the result.
Corollary 13. The center eigenspace X 0 corresponding to the purely imaginary eigenvalues of A is given by with dim X 0 = dim X 0 + p.
Proof. By Proposition 10 we have the disjoint union σ(A ) = (σ(A ) \ {0}) ∪ {0}. Using this and Theorem 12 we first compute the center eigenspace for A as and then we apply j −1 from (36) to both sides of this equality.
In Sections 3.4 to 3.6 we will consider nonlinear problems on a real Banach space. In this case spectral analysis must be preceded by complexification, see Remark 1 and in particular [14, last part of Section IV.2].
Proposition 14. Suppose that X = Y C is a complexification of a real Banach space Y and let Y 0 ⊆ Y be the real center eigenspace associated with X 0 . Then the real center eigenspace Y 0 associated with where Q : R p → Y is a bounded linear operator. Furthermore, ι : Y 0 → Y 0 × R p defined by ι(ψ, α) := (ψ − Qα, α) is a linear isomorphism.
Proof. X is naturally identified with Y C where Y = Y × R p . Let P Λ with range X 0 be the spectral projector on X for the spectral set Λ = Λ of all purely imaginary eigenvalues of A C . A direct generalization of [14,Corollary IV.2.19] implies that P Λ is the complexification of a projector P Y Λ on Y and the range Y 0 of P Y Λ -identified with a subspace of Y -is the real center eigenspace for A. Also, Γ 0 on X is self-conjugate by (41). Together with Corollary 13 this implies (45). It is easily verified that the linear operator ι is an isomorphism.
Remark 15. We will not make a notational distinction between the real and complex center eigenspaces, indicating both X 0 and Y 0 with X 0 and both X 0 and Y 0 with X 0 , respectively. We hope that the underlying scalar field will be clear from the immediate context.

Nonlinear perturbation
Let T 0 be a C 0 -semigroup on a real Banach space X that is sun-reflexive with respect to T 0 . Introduce a nonlinear perturbation R : X → X of the form where R : X → X is C k -smooth, satisfying We associate with T and R the integral equation We expect all nontrivial dynamics to be contained in the first component, and this is indeed the case: where α ∈ R p denotes the constant value of u p .
Proof. We use Proposition 9 and (36). For any continuous function u = (u, u p ) : I → X we compute while another computation shows that From the above together with the definition (38) of W (t) we see that (48) is equivalent to the system for −∞ < s ≤ t < ∞. The statement now follows.

Parameter-dependent local center manifolds
We consider again a C 0 -semigroup T 0 on a real Banach space X that is sun-reflexive with respect to T 0 . We also assume that T 0 is eventually compact and L in (37) is compact, so Proposition 11 implies that T is eventually compact. If furthermore the nonlinearity R satisfies (47), then all conditions are fulfilled for the application of the center manifold theory from [14, Chapter IX] to (48). Therefore, if the generator A of T has 1 ≤ n 0 < ∞ purely imaginary eigenvalues, counting algebraic multiplicities, then by Proposition 14 the real center eigenspace X 0 has dimension n 0 +p. There exists a C k -smooth local center manifold W c loc in X that is tangent at the origin to X 0 . In fact, Proposition 14 implies that W c loc is the image of a C k -smooth map where U ⊆ X 0 and U p ⊆ R p are neighborhoods of the origin. Since (46) has a zero in the second component, it follows from [14, (5.1) in Section IX.5] that C has the form where C : U × U p → X is the first component function.

Definition 17.
The image W c loc (α) := C(U, α) is a local center manifold for (49) at α ∈ U p . It is a direct consequence of the above definition that for every α ∈ U p we can parametrize W c loc (α) by coordinates on the real center eigenspace X 0 that depend C k -smoothly on α. This will be important for the discussion of the normalization method following (55) in Section 4.
Proposition 18. If α ∈ U p is sufficiently small then W c loc (α) is locally positively invariant for the semiflow generated by (49).
Proof. Let Σ and Σ be the semiflows generated by (48) and (49), respectively. By Proposition 16, and for all s in a common interval of existence I ψ,α . By [14, Theorem IX.5.3(i)] there exists δ > 0 such that if (ψ, α) ∈ W c loc and if then Σ(t, (ψ, α)) ∈ W c loc which by (51) implies that Σ(t, ψ) ∈ W c loc (α). We note that if ψ ∈ W c loc (α) then by (50) it follows that (ψ, α) ∈ W c loc . Therefore, if |α| ≤ δ 2 and ψ ∈ W c loc (α) then implies that Σ(t, ψ) ∈ W c loc (α). This is precisely local positive invariance of W c loc (α) for Σ. Next we consider a solution u : I → X of (49) that lies in W c loc (α). By Proposition 16 the function u = (u, α) : I → X is a solution of (48). Also, since u lies in W c loc (α) we see from (50) that u lies in W c loc and therefore satisfies the differential equation cf. (9). By (36) and (46) and Proposition 7 the first component of this equation gives the differential equation that is satisfied by u.
In summary, Theorem 19 (Parameter-dependent local center manifold). Let T 0 be an eventually compact C 0semigroup on a sun-reflexive real Banach space X and let T be the C 0 -semigroup on X defined by (3) where L is a compact perturbation. Suppose that the generator A of T has 1 ≤ n 0 < ∞ purely imaginary eigenvalues with corresponding n 0 -dimensional real center eigenspace X 0 . Furthermore, assume that R in (49) is C k -smooth and (47) holds. Then there exists a C k -smooth map C : U × U p → X defined in a neighborhood of the origin in X 0 × R p and such that for every sufficiently small α ∈ R p the manifold W c loc (α) := C(U, α) is locally positively invariant for the semiflow generated by (49) at parameter value α. Furthermore, any solution u : I → X of (49) that lies on W c loc (α) satisfies (52).

The special case of parameter-dependent classical DDEs
In this section we will formulate a corollary of Theorem 19 that applies specifically to the classical parameter-dependent DDE (28). As in Section 2.4 our starting point is (32) with F = 0, in the unknown (x, µ) with initial condition (ϕ, α) in the state space X := X × R p where X := C([−h, 0], R n ). So, we interpret the first component of (53) as a DDE but the second component as an ODE. By comparison with (13) it is clear that the solution of the initial value problem for (53) defines with T 0 the eventually compact shift semigroup on X from (14) and I p the identity on R p . We specify the perturbations L and L p in (37) as where D 1 F (0, 0) and D 2 F (0, 0) are as in (31). Then L is of finite rank, so certainly it is compact. Also, we choose the nonlinear perturbation R in (46) as where G is defined by the splitting in (31). Then (33) implies that the conditions in (47) hold.
Corollary 20 (Parameter-dependent local center manifold for DDEs). Consider the particular case of the classical DDE in (28), where F : X ×R p → R n is C k -smooth for some k ≥ 1 with F (0, 0) = 0. Let T be the C 0 -semigroup on X corresponding to the linearization of (28) at 0 ∈ X for the critical parameter value α = 0. Suppose that the generator A of T has 1 ≤ n 0 < ∞ purely imaginary eigenvalues with corresponding n 0 -dimensional real center eigenspace X 0 . Then there exists a C k -smooth map C : U × U p → X defined in a neighborhood of the origin in X 0 × R p and such that for every sufficiently small α ∈ R p the manifold W c loc (α) := C(U, α) is locally positively invariant for the semiflow generated by (28) at parameter value α.
Furthermore, if the history x t associated with a solution of (28) exists on some nondegenerate interval I and x t ∈ W c loc (α) for all t ∈ I, then u : I → X defined by u(t) := x t is differentiable and satisfies where A is the weak generator of T and the operator G : X × R p → R n defined by (31), is the nonlinear part of F . 4 Normal forms on the parameter-dependent center manifold The normalization technique described in this section goes back to [8]. In [42] it was applied to obtain expressions for the critical normal form coefficients of all generic codimension one and two bifurcations of equilibria in ODEs, also see [43, §8.7]. In this context, these expressions are independent of the (finite) dimension of the phase space and they involve only critical eigenvectors of the Jacobian matrix and its transpose as well as higher order derivatives of the right-hand side at the critical equilibrium. These properties make them suitable for both symbolic and numerical evaluation.
In [46] the same technique was applied to parameter-dependent normal forms to start the continuation of nonhyperbolic cycles emanating from generalized Hopf, fold-Hopf and Hopf-Hopf bifurcation points of ODEs. The resulting predictors were implemented in the freely available software package MatCont, a MATLAB toolbox for continuation and bifurcation analysis of finite dimensional dynamical systems. This makes it possible to verify transversality conditions and to initialize the continuation of the nonhyperbolic cycles mentioned above. A similar switching problem for iterated maps was solved earlier in [27].
In [37] the normalization technique was lifted to an infinite dimensional setting, providing explicit expressions for the critical normal form coefficients of generic codimension one and two equilibrium bifurcations in DDEs. These expressions were partially implemented [60] in the software DDE-BifTool. This package can be considered as the DDE equivalent of MatCont in command line mode.
In this section we extend the normalization method from [37] to include parameters. Suppose 0 ∈ X is a stationary state of (28) at the critical parameter value 0 ∈ R p and assume there are n 0 ≥ 1 eigenvalues on the imaginary axis, counting algebraic multiplicities. Let P 0 be the corresponding real spectral projector on X, so the range X 0 of P 0 is the real n 0 -dimensional center eigenspace. Corollary 20 applies to give a parameter-dependent local center manifold W c loc (α) for (28). We allow for the introduction of a new parameter β defined in a neighborhood of 0 ∈ R p such that α = K(β) for some locally defined C k -diffeomorphism K : R p → R p that is to be determined below, up to a certain order. If u : I → X with u(t) := x t ∈ W c loc (α) is as in Corollary 20, then u is differentiable on I and satisfies where R encodes the nonlinear part of F as in (54). Choose a basis Φ of X 0 and let H : R n0 × R p → X be a locally defined C k -smooth parametrization of W c loc (α) with respect to Φ and in terms of the new parameter β, see the remark following Definition 17. For every t ∈ I we define v(t) ∈ R n0 as the coordinate vector of P 0 u(t) with respect to Φ. Then v : I → R n0 satisfies a parameter-dependent ordinary differential equation of the forṁ where the C k -smooth vector field on the right has been expanded up to some sufficiently large -but finite -order. The multi-indices ν and µ have lengths n 0 and p, respectively. We assume that (56) is in parameter-dependent normal form, up to a certain order. Since H parametrizes W c loc (α), with both u and v depending on the parameter, although this is left implicit in the notation. Substituting the above relation into (55) produces the homological equation withv given by the parameter-dependent normal form (56). The unknowns in (HOM) are H, K and the coefficients g νµ from (56). They are determined, up to a certain order, by the assumption that (56) is in normal form. For r, s ≥ 0 with r + s ≥ 1 we denote by D r 1 D s 2 F (0, 0) : X r × [R p ] s → R n the mixed Fréchet derivative of order r + s, evaluated at (0, 0) ∈ X × R p , with the understanding that at most one of the factor spaces X r or [R p ] s is absent if either r = 0 or s = 0. We expand the nonlinearity R as where ϕ (r) := (ϕ, . . . , ϕ) ∈ X r and α (s) := (α, . . . , α) ∈ [R p ] s . The mappings H and K can be expanded as Substituting (56)-(58) into (HOM), collecting coefficients of terms v ν β µ from lower to higher order and solving the resulting linear systems, one can solve recursively for the unknown coefficients g νµ , H νµ and K µ by applying the Fredholm alternative and taking ordinary or bordered inverses, as explained in Section 2.6.

Predictors for normal forms on center manifolds
Starting from this point, we will focus exclusively on two-parameter DDEs, i.e. we will have p = 2. This enables the initialization of codimension one equilibrium and cycle bifurcations emanating from the generalized Hopf, fold-Hopf, Hopf-Hopf and transcritical-Hopf codimension two bifurcation points. The corresponding steady state and cycle bifurcation curves can then be continued in two parameters using the continuation capabilities of DDE-BifTool, providing an unfolding of the singularity. In this section we list known asymptotics for codimension one nonhyperbolic cycles emanating from generalized Hopf, fold-Hopf and Hopf-Hopf bifurcations obtained in [46]. The asymptotics are derived from the Poincaré normal forms which are obtained by considering near-identity changes of coordinates generated by homogeneous polynomial functions without using time reparametrization. These normal forms are therefore ready to be used in conjunction with the homological equation (HOM) where all the time derivatives are assumed to be with the same unit of time. Following the same method as in [46], we also derive asymptotics for codimension one nonhyperbolic cycles emanating from the transcritical-Hopf bifurcation. We also provide asymptotics for the codimension one equilibrium bifurcations emanating from the degenerate Hopf points.

Neimark-Sacker curve
A Hopf bifurcation in the amplitude system (69), i.e. when the trace of the Jacobian matrix of the amplitude system vanishes but the determinant of this matrix is positive, corresponds to a Neimark-Sacker bifurcation in the original system. Let ρ = and suppose that g 011 (0) Re(g 110 (0)) < 0. Then by the implicit function theorem, we obtain the second-order predictor for the Neimark-Sacker curve. It follows from (69) that the period of the cycle is approximated to the same order by Here (z 0 , β 1 , β 2 ) are as in (70) and Remark 22. Notice that the cubic terms g 210 (0) and g 300 (0) do not appear in the predictor. However, these coefficients are needed to determine the stability of the two-dimensional torus, see [43, Section 8.5] and (106).

Fold and Hopf curves
For the approximation to the fold and Hopf curves it is sufficient to consider the second order terms in (z 0 , ρ) in the amplitude system in (69). We obtain three equilibrium solutions given by It follows that the fold curve is approximated by and the Hopf curve by for | | small.

Hopf curves
There are two semi-trivial equilibria of the amplitude system of (73). Translating to the original system provides the Hopf bifurcation curves

Transcritical-Hopf bifurcation
A majority of papers in which fold-Hopf bifurcations in DDEs are studied, deals with models where the steady state remains fixed under variation of parameters. In this case the unfolding is not given by (67) anymore and we have to consider the smooth local normal form The bifurcation analysis can be carried out similar to the fold-Hopf case, see [29] and [40]. An alternative approach is presented in [62]. In contrast with the fold-Hopf bifurcation, there are in general two Neimark-Sacker bifurcation curves. Furthermore, the fold bifurcation curve becomes a transcritical bifurcation curve, and meets the Hopf bifurcation curve transversally.
Under the assumption of the same transversality condition as in the fold-Hopf bifurcation we introduce new parameters (β 1 (α), β 2 (α)) to obtain the truncated normal form with b 1 (0) = 0. Letting z 1 = ρe iψ and separating the real and imaginary parts yields the three dimensional system (76)

Hopf and transcritical bifurcations curves
The transcritical bifurcation curve in the normal form is obtained by substituting ρ = 0 in the amplitude system of (76). Then β 2 is unrestricted and z 0 = −β 1 /g 200 (0). The transcritical bifurcation curve is therefore given by To obtain a predictor for the Hopf bifurcation curve we truncate (76) to the second order. We obtain a trivial equilibrium (z 0 , ρ) = (0, 0), a semi-trivial equilibrium (z 0 , ρ) = (− β1 g200(0) , 0) and a nontrivial equilibrium It follows that the Hopf bifurcation curves are approximated by

Coefficients of parameter-dependent normal forms
Using the method outlined in Section 4, we derive here the coefficients needed for the predictors of the nonhyperbolic equilibria and cycles emanating from generalized Hopf, fold-Hopf, Hopf-Hopf and transcritical-Hopf bifurcations, see Section 5. While doing so, we also obtain the critical normal form coefficients, which were first derived in [37]. Recall, that we conside only two-parameter DDEs, so that p = 2.
For the derivation of the coefficients in this section it is sufficient to expand and truncate the nonlinearity R and the parameter-mapping K in (57) and (58), respectively, as follows Here u ∈ X, while α, β ∈ R 2 , and B, A 1 , J 2 , C, B 1 , D, C 1 and E are the standard multilinear forms arising from the expansion of F (u, α) (or G(u, α)). For example, These forms are R n -valued on real arguments, while they are linearly extended ('complexified') to C n -valued ones on complex-valued arguments. Finally, we introduce We assume in all situations that ϕ 0 = 0 is a steady state of (28) at the critical parameter value α 0 = 0 ∈ R 2 . Explicit formulas to compute the multilinear forms for the simplest DDE (1) are given in Section 7.

Generalized Hopf bifurcation
Since the eigenvalues (59) are simple, there exist eigenfunctions ϕ and ϕ such that The eigenfunctions ϕ and ϕ are explicitly given by (22) and (23) with q ∈ C n and p ∈ C n satisfying Any point y ∈ X 0 in the real critical eigenspace can be represented as where z = ϕ , y . Therefore, the homological equation (HOM) can be written as where β = (β 1 , β 2 ),ż is given by the normal form (61) and H admits the expansion For the predictors derived in Section 5.1 we need the truncated parameter-dependent normal forṁ where Re(c 1 (0)) = 0.

Critical normal form coefficients
We start by calculating the critical normal form coefficients following [37]. Collecting the coefficients of the quadratic terms z 2 and zz in the homological equation yields two nonsingular linear systems: They are solved using Lemma 3 to give For the cubic terms, the system corresponding to z 3 is also nonsingular, On the other hand, the system corresponding to z 2z is singular, The Fredholm solvability condition (FSC) requires that and from Lemma 4 we then obtain the unique solution satisfying ϕ , H 2100 = 0 as We continue by collecting the coefficients corresponding to the fourth-order terms z 3z and z 2z2 in the homological equation. The corresponding nonsingular systems may be solved using Lemma 3 and the fact that Re(c 1 (0)) = 0. For H 2200 this easily gives but for H 3100 the solution is a bit more subtle. The linear system is The critical normal form coefficient c 2 (0) is calculated by applying (FSC) to the singular linear system corresponding to the fifth-order term z 3z2 in the homological equation. This gives The second Lyapunov coefficient is now given by 2 (0) = 1 ω0 Re(c 2 (0)).

Parameter-related coefficients
Next we derive the parameter-related coefficients that provide a linear approximation to the parameter transformation. Following [46] and [2] we first expand the eigenvalue λ(α) and c 1 (α) in the normal form (60) in the original parameters α and truncate to fourth order, The parameters α and β are related via The homological equation (HOM) becomes where H admits the expansion and R is given by (78). Notice that the coefficients H jkµ with µ = (00) in (84) coincide with those in the expansion (81). Collecting the coefficients of the terms α and zα in the homological equation yields the systems where µ = (10), (01) and v 10 = (1, 0) T , v 01 = (0, 1) T . We solve these systems using Section 2.6. By the first part of Lemma 3 the first system has the (constant) solutions

Using Lemma 4 we obtain the solutions
for the second equation. To determine γ 2µ we first collect the coefficients corresponding to the z 2 α and zzα terms in the homological equation. We obtain the equations Lemma 3 implies that solutions of the first two equations are given by Applying (FSC) to z 2z α terms in the homological equation results in

Hopf and LPC predictors
Now we are ready to specify the predictors for the original parameter-dependent DDE (28). To approximate the Hopf parameter values α and the corresponding equilibrium, we merely substitute β from (65) into (82), and then put the result together with z = 0 into the expansion (84).
To approximate the LPC parameter values, we substitute β from (63) into (82). The cycle period is approximated by (64) with (83). To obtain a predictor for the periodic orbit in the phase space, we set z = e iψ into (84) using the obtained α values. Truncating to the second order in then yields

Fold-Hopf bifurcation
Since the eigenvalues (66) are simple, there exist eigenfunctions ϕ 0,1 and ϕ 0,1 satisfying as well as the mutual normalization condition The eigenfunctions ϕ 0,1 and ϕ 0,1 can be explicitly computed using (22) and (23) with q 0 ∈ R n , q 1 ∈ C n , p 0 ∈ R n and p 1 ∈ C n satisfying Any point y ∈ X 0 in the real critical eigenspace can be represented as where z 0 = ϕ 0 , y and z 1 = ϕ 1 , y . Therefore, the homological equation (HOM) can be written as where z = (z 0 , z 1 ,z 1 ), β = (β 1 , β 2 ) andż is given by the normal form (68). Here, the mapping H admits the expansion and the functions K and R are as in (78) and (79), respectively.

Critical normal form coefficients
We start by computing the critical normal form coefficients following [37]. Collecting the quadratic terms z 2 0 , z 2 1 , z 0 z 1 and z 1z1 we obtain one nonsingular and three singular linear systems. By (FSC) the singular systems are consistent if and only if This yields the three quadratic normal form coefficients. The corresponding solutions may then be obtained using Lemmas 3 and 4. Namely, For the four remaining cubic normal form coefficients, we collect the coefficients of the resonant terms This yields four singular linear systems. As before, by (FSC) these systems are consistent if and only if

Hopf, fold, and Neimark-Sacker predictors
To approximate the fold and Hopf curves and their corresponding equilibria, one should substitute the expressions for β and the equilibrium coordinates given in Section 5.2.2 into the expansions (86) and (79).

Hopf-Hopf bifurcation
Since the eigenvalues (71) are simple, there exist eigenfunctions ϕ 1,2 and ϕ 1,2 , satisfying the mutual normalization conditions The eigenfunctions ϕ 1,2 and ϕ 1,2 can be explicitly computed using (22) and (23) with q 1,2 ∈ C n and p 1,2 ∈ C n such that both as well as Any point y ∈ X 0 in the real critical eigenspace can be represented as where z 1 = ϕ 1 , y and z 2 = ϕ 2 , y . Therefore, the homological equation (HOM) can be written as where z = (z 1 ,z 1 , z 2 ,z 2 ), β = (β 1 , β 2 ) andż is given by the normal form (72). The mapping H admits the expansion and the functions K and R are as in (78) and (79), respectively.
Collecting the coefficients of the quadratic terms |z 1 | 2 , z 2 1 , z 1 z 2 , |z 2 | 2 , z 1z2 and z 2z1 in the homological equation, we obtain six nonsingular linear systems. By Lemma 3 their solutions are The desired cubic critical normal form coefficients are obtained by collecting the coefficients of the resonant cubic terms z 1 |z 1 | 2 , z 1 |z 2 | 2 , |z 1 | 2 z 2 and |z 2 | 2 z 2 in the homological equation. This leads to four singular linear systems. By (FSC) these systems are solvable if and only if

Hopf and Neimark-Sacker predictors
To approximate the Hopf curves and their corresponding equilibria, one should substitute the expressions for β and the equilibrium coordinates given in Section 5.3.1 into the expansions (92) and (79).

Transcritical-Hopf bifurcation
Compared with the fold-Hopf bifurcation in Section 6.2, the eigenvalues, eigenfunctions, the homological equation, and the functions H, K and R remain unchanged. It is only the ODE on the center manifolḋ z that changes to the normal form (75). Furthermore, also the critical normal form coefficients for the transcritical-Hopf bifurcation remain the same as for the fold-Hopf bifurcation. Therefore, we proceed only with the parameter-related equations.

Neimark-Sacker predictors
The predictors for the Hopf and transcritical bifurcation curves, as well as those for the Neimark-Sacker bifurcation curves (including the cycle periods), can be easily obtained using the asymptotics from Sections 5.4.1 and 5.4.2. In particular, to approximate the periodic orbits along the Neimark-Sacker curves, we substitute z 1 = e iψ and (77) into (86). This gives the following linear approximations: 7 Computation of derivatives for discrete DDEs All predictors described in the previous sections are implemented in version 3.2a of DDE-BifTool for models of the type (1). The discrete DDE (1) is a particular instance of (28) with h = τ m and where the linear evaluation operator Ξ : X → R n×(m+1) is defined by with the convention τ 0 := 0. In particular, by the chain rule, with M j := D 1,j f (0, 0) ∈ R n×n the partial derivative of f at the origin with respect to its jth state argument. So, if (1) has a equilibrium at the origin for α = 0, then the linear part of the splitting (11) at α = 0 is precisely the right-hand side of the above equation. Therefore ζ : [0, h] → R n×n must be such that Hence ζ has jump discontinuities M j at the points τ j for j = 0, . . . , m and is constant otherwise. So, in this case the characteristic matrix (20) is given by The multilinear forms appearing in (78) can be expressed in terms of the derivatives of the function f : R n×(m+1) × R p → R n from (1). For r, s ≥ 0 with r + s ≥ 1 the mixed derivative of order r + s of f at (0, 0) is an (r + s)-linear form on [R n×(m+1) ] r × [R p ] s , with the understanding that at most one factor may be absent in case r = 0 or s = 0. Let Q, Q 1 , . . . , Q r be matrices in R n×(m+1) and let α, α 1 , . . . , α s be vectors in R p . Then this derivative acts as where the multidimensional sum runs over 1 ≤ j 1 , . . . , j r ≤ n, 0 ≤ k 1 , . . . , k r ≤ m, 1 ≤ 1 , . . . , s ≤ p.

Examples
In this section we will demonstrate the correctness of the normal form coefficients and the accuracy of the predictors in four different models. We do this twofold. Firstly, by comparing the predictors in parameter-space with the computed in DDE-BifTool bifurcation curves, and, secondly, by performing simulations near the bifurcation point under consideration. The simulation is done either with the build-in routine dde23 of MATLAB or with the Python package pydelay [24]. The latter gives significant speed performance when considering simulation over longer time intervals. This usually is the case when one wants to demonstrate the existence of stable invariant manifolds. Since in this section only the main results are given, we provide details (including simulation results) in the Supplement. Furthermore, the source code of the examples has been included into the DDE-BifTool software package. This will hopefully provide a good starting point when considering other models.

Generalized Hopf bifurcation in a coupled FHN neural system with delay
In [65] the following system is considered Here (u 1 , u 2 ) is the completely synchronous solution of the three coupled FitzHugh-Nagumo (FHN) neuron system where α, β measure the synaptic strength in self-connection and neighborhood-interaction, respectively. The parameters b and are assumed to be positive such that 0 < b < 1 and 0 < 1. The function f is a sufficiently smooth sigmoid amplification function and τ > 0 represents the time delay in signal transmission. For the derivation of (100) from the system (101), as well as for stability conditions of the completely synchronous solution, we refer to [65]. In that article a generalized Hopf point was analyzed using the traditional formal adjoint method and the two-step center manifold reduction, see [31]. Numerical simulations where made to confirm their results. For this (β, α) are taken as the unfolding parameters and the parameters b = 0.9, ε = 0.08, c = 2.0528, d = −3.2135, τ = 1.7722 are fixed. The sigmoid amplification function f (u) = tanh(u) is used.  According to [65], a generalized Hopf point is present at the origin with the parameter values (β, α) = (1.9, −0.9710). We took this point and calculated its stability and the corresponding normal form coefficients. Although we do confirm that the point under consideration is a Hopf point, the first Lyapunov coefficient does not vanish and we conclude that the point cannot be a generalized Hopf point. However, the simulation in [65] do suggest a generalized Hopf point for nearby parameter values. Therefore we continued the Hopf point in (β, α). Then a generalized Hopf point is located at (β, α) = (1.9, −1.0429) with negative second Lyapunov coefficient 2 (0) = −15.6733, indicating the existence of a stable steady state inside a unstable cycle, which in turn is located inside a stable cycle. We remark that the second Lyapunov coefficient found in [65] is positive. This contradicts the simulation of the dynamics made in the same article. Indeed when the second Lyapunov coefficient is positive a time-reversal must be taking into account when considering the bifurcation diagram in the case the second Lyapunov coefficient is negative, see [43]. Then the situation of a stable steady state inside a stable cycle (separated by an unstable cycle) does not occur.
Using the predictors from Section 5.1 the fold and LPC bifurcation curves emanating from the generalized Hopf point. In Figure 1 the resulting bifurcation diagram is shown.

Fold-Hopf bifurcation of the Rose-Hindmarsh model with time delay
In [48] a Rose-Hindmarsh model [33,34] with time delay in the self-feedback process, which takes the form is considered. Here x represents membrane potential, y represents a recovery variable, z denotes the adaption current, and a, b, c, d > 0, S and χ are real constants. The external current I app and r are control parameters, and τ denotes the synaptic transmission delay. The constants a, b, c, d, χ and r are fixed. Let (x , y , z ) be a steady state of (102), then The conditions for a fold-Hopf bifurcation have been derived in [48] analytically. Indeed, let S be arbitrary and set Then a fold-Hopf bifurcation occurs when where k = 0, 1, 2, . . . and In [48] the parameters values are fixed. It follows that a fold-Hopf bifurcation is located at x = 0.1308, S = −0.57452592, τ = 5.768830916, and I app , (y , z ) given by (104) and (103), respectively. To unfold the singularity the 'parameters' (x , S) are used, see [48]. Here we will take the more natural unfolding parameter (I app , S). Calculating the stability with DDE-BifTool gives the eigenvalues 0.001 ± 1.0081i, −0.000 + 0.000i.
Since s = 1 and θ(0) < 0 , a global bifurcation curve or invariant tori are present for parameters sufficiently close to the bifurcation, see [43, page 342]. However, since the sign of e(0) is positive the tori are unstable. Thus according to our analysis the simulated torus in [48] cannot be attributed to the fold-Hopf bifurcation.  Thus the sign of s and θ(0) remain unchanged. However, since the sign of e(0) is negative, there is a time reversal to take into account. Therefore, we expect a stable torus to be present for nearby parameter values. Using the predictors from Section 5.2, we successfully continued the fold, Hopf, and Neimark-Sacker bifurcation curves emanating from the point, see Figure 2.

Hopf-Hopf and generalized Hopf bifurcations in Active control system
Active control system is used to control the response of structures to internal or external excitation. The mathematical model with time delay can be described as follows [50] Here x(t) is the displacement of the controlled system, m > 0 is the mass, c and k are the damping and the stiffness, respectively, τ is the time delay represented in the relative displacement feedback loop and in the relative velocity feedback loop, u and v are feedback strengths, respectively, andf represents the external excitation. Let t * = k/mt, ζ = c/2m m/k, g u = u/k, g v = v/m m/k and f (t) =f (t)/k . Then equation (107) becomes where the asterisks are omitted for simplicity. Following [16] and [50] we consider the case when f is replaced by a nonlinear position time delay feedback given by βx 3 (t − τ ), see also [63]. As in [16] we fix the parameters g u = 0.1, g v = 0.52, β = 0.1 and take ζ and τ as control parameters. Letẋ(t) = y(t), then we obtain ẋ(t) = τ y(t), Here the delay is scaled by using the transformation of time t → t/τ . In this way the delay can treated as an ordinary parameter. The trivial equilibrium undergoes a Hopf-Hopf bifurcation at the parameter values see [16] for the derivation. Using DDE-BifTool we manually construct the Hopf-Hopf point and compute its stability and normal form coefficients. We obtain the eigenvalues 0.0000 ± 4.5275i and −0.0000 ± 7.6449i. The quadratic critical normal form coefficients are g 2100 (0) = −0.0915 + 0.1214i, g 1110 (0) = 0.2151 + 0.3876i, g 1011 (0) = −0.3084 + 0.4096i, g 0021 (0) = 0.1813 + 0.3268i.

From
(Re g 2100 (0))(Re g 0021 (0)) = −0.0166 < 0, we conclude that this Hopf-Hopf bifurcation is of 'difficult' type, see [43]. Furthermore, since the quantities are such that θ < 0, δ < 0, θδ > 1 it follows that we are in case VI. We continue the Neimark-Sacker and Hopf bifurcation curses emanating from the Hopf-Hopf point using the predictors from Section 5.3. In Figure 3 a close-up is given near the Hopf-Hopf point comparing the computed curves with the predictors in parameter space.
Using the detection capabilities of DDE-BifTool one additional Hopf-Hopf point and three generalized Hopf points are located on the continued Hopf branches. The normal form coefficients of the second Hopf-Hopf point are such that (Re g 2100 (0))(Re g 0021 (0)) = 1.7331e−04 > 0   (109) in an active control system with time delay given by (108). There are two supercritical Hopf curves (blue) and two Neimark-Sacker curves (yellow). We see that the predictors (dotted) give good approximations near the codimension two point. and θ ≥ δ > 0, θδ > 1.
We conclude that we are in case I of the 'simple' type, see [43, page 360]. Therefore, no stable invariant two-dimensional torus is predicted for nearby parameter values, only two stable period orbits expected. Using the predictors from Sections 5.1 and 5.3 we can easily continue the codimension one cycle bifurcations from the located degenerate Hopf points, showing complicated bifurcation diagram in Figure 4.

Transcritical-Hopf bifurcation in Van der Pol's oscillator with delayed position and velocity feedback
In [3] a generalization of Van der Pol's oscillator with delayed feedback is considered. Here g ∈ C 3 satisfies the conditions g(0, 0) = 0, gẋ(0, 0) = a and g x (0, 0) = b. The linearization of equation (110) around the trivial solution x = 0 gives From which we obtain the characteristic equation then the characteristic equation has a simple zero and a pair of purely imaginary roots λ = ±iω 0 . Here ω 0 and τ 0 are defined by We set the function g to  To analyze the system with DDE-BifTool we set y(t) =ẋ(t) and transform the time with t → t/τ to obtain the two-component system Here we introduced the unfolding parameters (µ 1 , µ 2 ) := (b − 1, τ − τ 0 ) to translate the singularity to the origin. One immediately sees that the trivial equilibrium (ẋ, x) = (0, 0) is an equilibrium for all parameter values (µ 1 , µ 2 ). Therefore, the parameter-dependent normal form for the generic fold-Hopf cannot be used here. Instead the normal form for the transcritical-Hopf bifurcation must be used. Using DDE-BifTool we compute the stability and the normal form coefficients. The leading eigenvalues are 0.000 + 0.000i and −0.000 + 2.4539i, where 2.4539 ≈ ω 0 τ 0 , see (112). Furthermore, the normal form coefficients are such that g 011 (0) × Re (g 110 (0)) = 0.4241 × Re (−0.1337 + 0.2672i) < 0.
Therefore, there are two Neimark-Sacker bifurcation curves predicted, see Section 5.4. Using the predictors from Section 5.4 we continue the transcritical, Hopf and Neimark-Sacker bifurcation curves emanating from the transcritical-Hopf bifurcation point. In Figure 5 the bifurcation diagram is shown.

Concluding remarks
We have provided explicit formulas for the normal form coefficients needed to initialize codimension one equilibrium and nonhyperbolic cycle bifurcations emanating from generalized Hopf, fold-Hopf, Hopf-Hopf and transcritical-Hopf points in DDEs. Applications to four different models were given, confirming the correctness of the derivation of the normal form coefficients and the asymptotics. A paper providing a second-order predictor for the homoclinic orbits emanating from the generic and transcritical codimension two Bogdanov-Takens bifurcations in DDEs, along the lines of [45], is in preparation.
Our proof of the existence of a smooth parameter-dependent center manifold is given in the general context of perturbation theory for dual semigroups (sun-star calculus). Consequently the applicability of this result extends beyond classical DDEs, although here we did restrict to the case of an eventually compact C 0 -semigroup on a sun-reflexive state space. It follows that the results from Sections 3.2 to 3.5 are valid as well for other classes of delay equations such as renewal equations (also known as Volterra functional equations) and systems of mixed type [11].
Furthermore, in [15,59] the technique was used to calculate the critical normal form coefficients for Hopf and Hopf-Hopf bifurcations occurring in neural field models with propagation delays. For these models sun-reflexivity is lost, which is typical for delay equations in abstract spaces or with infinite delay. However, it is often possible to overcome this functional analytic complication, so dual perturbation theory can still be employed successfully [12,13,38,59]. It has also been used in the context of semilinear hyperbolic systems [47].
It is demonstrated -at a formal level -in [55] that the normalization technique described in Section 4 still works for DDEs with state-dependent delays. However, as already mentioned in Remark 21, for DDEs the nonlinearity generally does not depend differentiably on the delay parameters. Therefore, in the case of state-dependent DDEs it is generally not possible to justify differentiation of the nonlinearity with respect to the state, let alone to rely on higher order smoothness. So, as far as we know there is still no proof of the validity of the normalization technique for state-dependent DDEs.
Returning to the setting of classical DDEs, the most obvious next challenge is to derive normal forms for bifurcations of periodic orbits by generalizing [9,10,44] to DDEs. The resulting formulas can then be implemented in DDE-BifTool to facilitate numerical bifurcation analysis of periodic orbits in classical DDEs.
supplementary materials for: Switching to nonhyperbolic cycles from codimension two bifurcations of equilibria of delay differential equations

March 19, 2019
In this supplement we provide walkthroughs of the examples given in Section 8 with DDE-BifTool 1 [21]. These walkthroughs enable other researchers to reproduce the results obtained in the main text.
Additionally, we will show the code used for simulation near the bifurcations points under consideration. Either using the build-in routine dde23 from MATLAB [53] or the Python package pydelay 2 [24]. Other DDE models, undergoing one of the degenerate Hopf bifurcations treated in this paper, can easily be studied by making minor modifications to the given code.
The focus will be on the initialization and continuation of the various codimension one equilibrium and cycles bifurcation curves emanating from the degenerate Hopf points and on simulation near the bifurcation points. For a complete overview of the capabilities and functionality for DDE-BifTool, we refer to the online tutorials files and also the manual and the references therein.
All code has been included into the DDE-BifTool package version 3.2a on the SourceForge repository and can be executed without the need to copy and paste. Note that the code is tested on MATLAB 2018b and GNU Octave 4.2.2. Different results may occur with other versions of MATLAB and GNU Octave.

S1 Generalized Hopf bifurcation in a coupled FHN neural system with delay
A completely synchronous solution of the three coupled FitzHugh-Nagumo (FHN) neuron system is given by the system (S1) see Section 8.1 and [48]. As before, we fix the parameters b = 0.9, ε = 0.08, c = 2.0528, d = −3.2135, τ = 1.7722, and take for f : R → R the sigmoid amplification function f (u) = tanh(u). The parameters (β, α) are used to unfold the singularity.
Remark 23. This demonstration can be found in the directory demos/tutorial/VII/FHN relative to the main directory of the DDE-BifTool package.

S1.1 Generate system files
Before we start to analyze the system with DDE-BifTool, we first create a system file. This file contains the definition of the system (S1), the standard derivatives needed for calculation of the eigenvalues and eigenvectors, the continuation of bifurcation points and cycles, and also the multilinear forms, see Section 7, used for the calculation of the coefficients of the critical and parameter-dependent normal forms. Alternatively, one can only supply the system itself, see Listing S2. Then finite difference is used to approximate the derivatives. However, this is less efficient and accurate, and therefore not recommended. A separate script gen_sym_FHN.m is used to create a system file. The most important parts of this script are listed and discussed below. pkg load symbolic end % % Create parameter names as strings and define fixed parameters % The demo has the parameters | beta | , | alpha | and | tau | parnames ={ ' beta ' , ' alpha ' , ' tau ' }; b = sym (0.9 , 'r ') ; epsilon = sym (0.08 , 'r ') ; c = sym (2.0528 , 'r ') ; d = sym ( -3.2135 , 'r ') ; % % Create symbols for parameters , states and delays states % | par | is the array of symbols in the same order as parnames . % Due to the following two lines we may , for example , % use either beta or par (1) to refer to the delay . syms ( parnames {:}) ; % create symbols for beta , alpha and tua par = cell2sym ( parnames ) ; % now beta is par (1) etc % % Define system using symbolic algebra % create symbols for u1 ( t ) u1 (t -tau ) , u2 ( t ) , u2 (t -tau ) syms u1 u1t u2 u2t du1_dt = -u1^3/3+( c + alpha ) * u1^2+ d * u1 -u2 +2* beta * tanh ( u1t ) ; du2_dt = epsilon *( u1 -b * u2 ) ; % % Differentiate and generate code ( multi -linear forms ) [ fstr , derivs ]= dde_sym2funcs (... The variable ddebiftoolpath is directed to the DDE-BifTool main folder, which should have been extracted somewhere on the computer. Here a path relative to the current working directory is used. Note that although we only use the parameters (β, α) as unfolding parameters, in the current version of DDE-BifTool, we also need to include the delay(s) in the list of parameters. After running the script, the function dde_sym2funcs creates two system files sym_FHN_mf.m and sym_FHN.m. The first file sym_FHN_mf.m implements the higher order derivatives as multilinear forms, as explained in Section 7, and therefore the file we will solely be using. The second file sym_FHN.m uses directional derivatives to implement the higher order derivatives. The directional derivatives approach formally allows the use of state-dependent delays, see [55]. Although both approaches yields (up to rounding errors) identical normal form coefficients, multilinear forms are much faster.

S1.2 Loading the DDE-BifTool package
Now that a system file is created we continue with DDE-BifTool to analyze (S1). The code in the following sections highlights the import parts of the file FHN.m. DDE-BifTool consists of a set of MATLAB routines. Thus, to start using DDE-BifTool, we only need to add DDE-BifTool directories to the search path.
Listing S1: Add DDE-BifTool scripts to the search path There are four subdirectories added to the search path: ddebiftool Containing the core files of DDE-BifTool.
ddebiftool_extra_psol An extension for enabling continuation of periodic orbit bifurcations for delay-differential equations with constant or state-dependent delay.
ddebiftool_extra_nmfm An extension for normal form computation.

S1.3 Set parameter names
The following code allows us to use ind.beta instead of remembering the index of the parameter β in the parameter array, and similarly for the other parameters. In this way, fewer mistakes are likely to be made and the code is easier to read.

S1.4 Initialization
Next, we set up the funcs structure, containing information about where the system and its derivatives are stored, a function pointing to which parameters are delays, and various other settings.
Listing S2: Define system without a system file Inspecting the output of the funcs handle gives. hjac : @ ( ord ) eps^(1/(2+ ord ) ) s y s _ u n fo l d in g _ pa r a me t e rs : [] tp_del : 0 sys_deri_provided : 0 sys_dirderi_provided : 0 The output shows that no derivative file is supplied. In this case, the derivatives are calculated using finite-difference approximations with the function dde_dirderiv. Again, we do not recommend using the latter approach. However, it can be useful for debugging purposes.

S1.5 Stability and normal form coefficients of the generalized-Hopf point
We manually specify a steady-state at the generalized-Hopf point found in [48] and calculate its stability. The eigenvalues confirm that the point under consideration is indeed a Hopf point. Next, we convert the steady-state point to a Hopf point and calculate the normal form coefficients with the function nmfm_genh, which implements the coefficients derived in Section 6.1. % % Calculate critical normal form coefficients hopf = p_tohopf ( funcs , stst ) ; method = df_mthod ( funcs , ' hopf ') ; hopf . stability = p_stabil ( funcs , hopf , method . stability ) ; genh = p_togenh ( hopf ) ; genh = nmfm_genh ( funcs , genh ) ; The normal form coefficients are stored in the genh.nmfm structure.

S1.6 Continue Hopf point
Since the simulations in [48] do indicate a generalized-Hopf point nearby, we continue the Hopf point. The continued branch hbr is shown in Figure S1a.

S1.7 Detect bifurcation points
To detect bifurcation points on the Hopf branch, we use the function LocateSpecialPoints.
[ hbr_wbifs , hopftests , hc2_indices , hc2_types ]=... L oc ateSpecialPoints ( funcs , hbr ) ; The MATLAB console shows the following output.  This function uses the file nmfm_POfold_from_genh_init.m in which the predictor, see Sections 5.1 and 6.1.3, is implemented. Using this function, a branch with three initial corrected cycles is created which is continued in the standard way, see Figure S1b.

S1.9 Calculate predicted periodic orbits
To compare the computed parameter values and periodic orbits on the branch lpcbr with the predictor, we again use the function C1branch_from_C2point, but with the additional argument predictor set to 1 and step to an interval of ε-values. Now the cycles are left uncorrected.

S1.10 Bifurcation diagram
The following code produces the bifurcation diagram presented in the main text, see Figure 1, and has been reproduced here in Figure S3. The figure was exported with the MATLAB and GNU Octave compatible package matlab2tikz, see [26].

S1.11 Plot comparing computed and predicted periodic orbits
Lastly, we create a plot to compare the computed and predicted periodic orbits. The resulting plot is shown in Figure S4. Note that the cycles shown have different underlying parameter values. Nonetheless, we see that the cycles are in good agreement.

S1.12 Simulation with MATLAB
Next, we simulate the dynamics near the generalized Hopf point. For this, we take a point in each of the three regions as shown in Figure S3. The following code, from the file FHN_simulation.m, uses the MATLAB function dde23.

S2 Fold-Hopf bifurcation in the Rose-Hindmarsh model with time delay
In [48] a Rose-Hindmarsh model [33,34] with time delay in the self-feedback process, Remark 24. This demonstration can be found in the directory demos/tutorial/VII/RH relative to the main directory of the DDE-BifTool package. Here, we omit the code to generate a system file. The system file sym_RH_mf.m has been generated with the script gen_sym_RS.m. Also, we assume that the DDE-BifTool package has been loaded as in Listing S1. The code in Sections S2.1 to S2.10 highlights the important parts of the file RH.m.

S2.1 Set parameter names and funcs structure
As in the previous example, we set the parameter names and define the funcs structure.

S2.2 Stability and normal form coefficients of the fold-Hopf point
We construct a steady-state at the fold-Hopf point and calculate its stability.  Figure S5: Simulation near the generalized Hopf point in the system (S1). In (a) we see a stable steadystate corresponding to a point in region I. When we enter region II, the stability of the steady-state is lost and a stable cycle appears, as seen in (b). In region III, there is a stable steady-state inside a stable cycle. This is confirmed in (c) and (d). In (c), the initial point of the orbit in blue is just outside the unstable cycle and converges to the stable cycles. The initial point of the orbit in red is just inside the unstable cycle and converges to the stable steady-state. In (d), the time series of these orbits are shown in the (t, u 1 (t)) plane. The function bifurcationvalues calculates (x , I app , τ ) according to the formulas as given in Section 8.2.
The MATLAB console shows the following output. We have a zero eigenvalue and a pair of purely imaginary eigenvalues. Furthermore, the remaining eigenvalues have negative real parts. Next, we calculate the normal form coefficients and the transformation to the center manifold with the function nmfm_zeho, which implements the coefficients as derived in Section 6.2. For this we need to set the argument free_pars to the unfolding parameter (I app , S). These coefficients will be used to start the continuation of the various branches emanating from the fold-Hopf point. The MATLAB console shows the following output. Since s > 0 and θ(0) < 0 global bifurcations or invariant tori are present. However, since the sign of e is positive the stability of the invariant tori will be unstable for nearby parameter values. It follows that the torus observed in [48] by simulations does not originate from the fold-Hopf point under consideration.

S2.3 Set bifurcation parameter range and step size bounds
Before continuing the various branches emanating from the fold-Hopf point, we create the variable brpars containing parameter bounds and maximal stepsizes.

S2.4 Continue NS, Hopf and fold branch
As in the previous example, we use the function C1branch_from_C2point to start to continue the branches emanating from the fold-Hopf point. Figure S6 is created using similar code as in Listing S3. We remark that even when there would be stable tori present for nearby parameter values, the window in which these tori would exist is quite small. Indeed, the parameter values would have to be below the Hopf curve and above the Neimark-Sacker curve to the left of the fold-Hopf point.

S2.5 Stable invariant tori
We change the parameters r = 1.4 and S = −8, while keeping the other fixed parameters as in (S3). Using the formulas given in Section 8.2, we calculate x , I app and τ . The coefficients s and θ(0) reveal that we are in case III of the fold-Hopf bifurcation, see [43, page 342].
Since the sign of e(0) is negative, there is a time reversal to take into account. Therefore, we expect stable tori to be present for nearby parameter values.

S2.6 Adjusting bifurcation parameter range
We adjust the variable brpars to reflect the current situation.

S2.7 Detect special points on the Hopf branch
Since the code to continue the Neimark-Sacker, Hopf and fold curves is identical to the code in Section S2.4, we continue with detecting bifurcations on the Hopf branch. The Hopf points on the branch hbr_wbifs will contain the normal form coefficients L1 and L2. These will be used to visualize the criticality of the Hopf points (sub or super) in the bifurcation diagram.

S2.8 Predictors
As in the previous example, we obtain predictors for the various branches simply by setting the argument predictor to 1 and the argument step to a range of ε-values when calling the function

S2.9 Bifurcation diagram
We plot the obtained curves and the predictor for the Neimark-Sacker and Hopf curve with the following code.

S2.10 Plots comparing computed and predicted periodic orbits
We create a plot to compare the computed and predicted periodic orbits.   The resulting plot is shown in Figure S8a. To compare the computed and predicted periods, we use the following code. The resulting plot is shown in Figure S8b.

S2.11 Simulation with pydelay
In this section we simulate the dynamics near the fold-Hopf point. The following code, from the file RH_simulation_torus.py, uses the Python package pydelay [24].
import numpy as np from pydelay import dde23 import matplotlib . pyplot as plt from mpl_toolkits . mplot3d import Axes3D   We will simulate the dynamics in region I and II, see Figure S7, where a stable periodic orbit and stable two-dimensional torus, respectively, should be present. For a point in region I, we will take the unfolding parameter values (I app , S) = (−18.886177304147466, −8.044197084548104).
As an initial condition, we use the constant history function with values of the location of the fold-Hopf point. We integrate the DDE on the time interval t ∈ [0, 10.000] with pydelay. In Figure S9 the time series of the components x,y and z and the orbit in (x, y, z)-space of the last 10 time steps are shown, clearly indicating a stable orbit. Next, we simulate the dynamics in region II. Therefore, we adjust the unfolding parameter values to (I app , S) = (−18.902420391705071, −8.045234985422740).
Furthermore, we increase the integration interval to t ∈ [0, 20.000]. Keeping the history function the same, we plot the last 1000 time steps, see Figure S9. We conclude that the dynamics near the fold-Hopf point are as predicted in [43]. x y (c) (d) S3 Hopf-Hopf and generalized Hopf bifurcations in Active control system In Section 8.3 we considered the following active control system ẋ(t) = τ y(t), which is used to control the response of structures to internal or external excitation, see [50]. The function f is substituted by βx 3 (t − 1) and the parameters We conclude that this Hopf-Hopf bifurcation is of 'difficult' type, since (Re g 2100 )(Re g 0021 ) = −0.0166 < 0, see [43]. Furthermore, the quantities θ = Re g 1011 Re g 0021 = −1.7009, δ = Re g 1101 Re g 2100 = −2.3517 are such that θ < 0, δ < 0, θδ > 0. It follows that we are in case VI.

S3.3 Set bifurcation parameter range and step size bounds
Before continuing the various branches emanating from the transcritical-Hopf point, we create the variable brpars containing parameter bounds and maximal stepsizes.

S3.5 Predictors
For comparison in the bifurcation diagram, we obtain predictors for the various branches by setting the argument predictor to 1 and step to a range of ε-values when calling the function C1branch_from_C2point.

S3.6 Bifurcation diagram
We plot the computed curves and the predictors for the Neimark-Sacker and Hopf curves with the following code.

S3.8 Detect special points on the Hopf branch
Using the detection capabilities from DDE-BifTool via the function LocateSpecialPoints, we detect one additional Hopf-Hopf point and three generalized Hopf points. % % Detect codimension two points on hopf_br1 [ hopf_br_wbifs , hopftests , hc2_indices , hc2_types ]=... L oc ateSpecialPoints ( funcs , hbr (1) ) ; The MATLAB console shows the following output. There are seven bifurcations detected on the Hopf branch hbr(1): four Hopf-Hopf bifurcations and three generalized Hopf bifurcations. However, by inspecting the parameters of the detected Hopf Hopf points suggest that there are only two distinct Hopf-Hopf points, which are connected with the same Hopf branch. By comparing the location of the underlying points confirms the premise. We are therefore not interested in continuing the Hopf curves emanating from these points. What we are interested in is continuing the Neimark-Sacker and limit point of cycle curves emanating from the second Hopf-Hopf point and the three generalized Hopf points, respectively. The normal form coefficients of the second Hopf-Hopf point are such that (Re g 2100 )(Re g 0021 ) = 1.7331e − 04 > 0 and θ ≥ δ > 0, δθ = 4 > 1.
We conclude that we are in case I of the 'simple' type, see [46, page 360]. Therefore, no stable invariant two-dimensional torus is predicted for nearby parameter values. We only expect to find two stable periodic orbits.

S3.9 Continuing third Neimark-Sacker bifurcation and limit point of cycle bifurcation curves
It turns out that we only need to continue one of the Neimark-Sacker bifurcation curves emanating from the second Hopf-Hopf point. Indeed, the other Neimark-Sacker bifurcation curve is given by

S3.10 Bifurcation diagram
We plot the computed degenerate Hopf points, the limit point of limit cycle curves, and the Neimark-Sacker curves with the following code.  (2) Figure S12 shows the resulting bifurcation diagram. There, we see that two generalized Hopf points are connected by the same limit point of cycles bifurcation curve.

S3.11 Simulation near Hopf-Hopf point with pydelay
In this last Section, we simulate the dynamics near the manually constructed Hopf-Hopf point at parameter values (109). As remarked before, the unfolding of the Hopf-Hopf point hoho is of 'difficult' type case VI. The normal form coefficients predict a stable invariant two-dimensional torus. Furthermore, this torus undergoes a bifurcation in which a three-dimensional torus is born. Since the DDE under investigation (108) does not contain any terms of order higher than three, the results for the normal form remain valid for the system. To confirm the unfolding in Figure 3, we fix the delay see Figure S13a. The cross-sections in Figures S13b and S14d are generated with the following Python code using pydelay, see also the file acs_simulation.py import numpy as np import pylab as pl from pydelay import dde23 from matplotlib import colors    (1) pl . scatter ( x_cros , y_cros , s =0.8 , color = ' royalblue ') pl . xlabel ( ' $x$ ') pl . ylabel ( ' $y$ ') pl . show () In Figure S13b there are two dots, corresponding to a stable period orbit. Crossing the curve ns1_br, a stable two-dimensional torus branches off, see Figure S14a. The torus still exists at ζ = ζ 3 , as seen in Figure S14b. Then, at ζ = ζ 4 , only slightly smaller than ζ 3 , a three-dimensional torus is observed, see Figure S14c. Lastly, Figure S14d shows the three-dimensional torus near the curve where the torus blows up.
Remark 26. This demonstration can be found in the directory demos/tutorial/VII/vdpo relative to the main directory of the DDE-BifTool package. Here, we omit the code to generate a system file. The system file sym_vdpo_mf.m has been generated with the script gen_sym_vdpo.m. Also, we assume that the DDE-BifTool package has been loaded as in Listing S1. The code in Sections S4.1 to S4.8 highlights the important parts of the file vdpo.m.

S4.1 Set parameter names and funcs structure
We set the parameter names and define the funcs structure.

S4.2 Stability and normal form coefficients of the transcritical-Hopf point
We construct a steady-state at the transcritical-Hopf point and calculate its stability.
% % Construct transcritical -Hopf bifucation point stst = dde_stst_create ( 'x ' ,[0;0] , ' parameter ' ,[0 0 1]) ; % Calculate stability method = df_mthod ( funcs , ' stst ') ; stst . stability = p_stabil ( funcs , stst , method . stability ) ; stst . stability . l1 The MATLAB console shows the following output. We have a zero eigenvalue and a pair of purely imaginary eigenvalues. Furthermore, the remaining eigenvalues have negative real parts. Next, we calculate the normal form coefficients and the transformation to the center manifold with the function nmfm_zeho, which implements the coefficients as derived in Sections 6.2 and 6.4. For this we need to set the argument free_pars to the unfolding parameter (µ 1 , µ 2 ). These coefficients will be used to start the continuation of the various branches emanating from the transcritical-Hopf point. The normal form coefficients are such that g 011 Re(g 110 ) = 0.4241 Re(−0.1337 + 0.2672i) < 0.
Therefore, there are two Neimark-Sacker bifurcation curves predicted, see Section 5.4.

S4.3 Set bifurcation parameter range and step size bounds
Before continuing the various branches emanating from the transcritical-Hopf point, we create the variable brpars containing parameter bounds and maximal stepsizes.

S4.4 Continuing Hopf, transcritical and Neimark-Sacker bifurcation curves
We use the function C1branch_from_C2point to continue the various branches emanating from the singularity.

S4.5 Detect special points on the Hopf branches
We continue with detecting bifurcations on the Hopf branches. The Hopf points on the branch hbr_wbifs(i)(i = 1, 2) will contain the normal form coefficients L1 and L2. These will be used to visualize the criticality of the Hopf points (sub or super) in the bifurcation diagram.

S4.6 Predictors
For comparison in the bifurcation diagram we obtain predictors for the various branches by setting the argument predictor to 1 and step to a range of ε-values when calling the function C1branch_from_C2point.

S4.7 Bifurcation diagram
We plot the obtained curves and the predictors for the Neimark-Sacker, Hopf, and transcritical curves with the following code.  Figure S15 shows the resulting bifurcation diagram.

S4.8 Plot comparing computed and predicted periodic orbits
Lastly, we create a plot to compare the computed and predicted periodic orbits.  Figure S14: Simulation with pydelay illustrating the branching of a three-dimensional torus from a two-dimensional torus. We refer to the text for further description.  Figure S15: Bifurcation diagram near the transcritical-Hopf bifurcation in the delayed Van der Pol's oscillator given by (113). There are two supercritical Hopf curves (blue), two subcritical Hopf curves (red), two Neimark-Sacker curves (yellow), and one transcritical curve (green). We see that the predictors (dotted) give a good approximation for nearby values. (c) (d) Figure S17: In (a) and (b), the periodic solution is shown to be present for parameter values in region III of Figure S15. In (c), the torus present in IV of Figure S15 is shown. In (d), a cross-section of the torus with y(t) = −0.009 in the phase-space (x, x(t − τ 0 ), y) is taken.