Coexistence of infinitely many large, stable, rapidly oscillating periodic solutions in time-delayed Duffing oscillators

TL;DR

This paper investigates the stability of rapidly oscillating periodic solutions in a time-delayed Duffing oscillator, revealing conditions for their stability or instability and illustrating these findings with numerical examples.

Contribution

It provides a detailed analysis of the stability and instability of large, rapidly oscillating solutions in delayed Duffing oscillators, including conditions for stabilization and destabilization.

Findings

Stable for certain parameter conditions

Unstable for other parameter conditions

Numerical simulations confirm theoretical predictions

Abstract

We explore stability and instability of rapidly oscillating solutions $x(t)$ for the hard spring delayed Duffing oscillator $$x''(t)+ ax(t)+bx(t-T)+x^3(t)=0.$$ Fix $T>0$. We target periodic solutions $x_n(t)$ of small minimal periods $p_n=2T/n$, for integer $n\rightarrow \infty$, and with correspondingly large amplitudes. Note how $x_n(t)$ are also marginally stable solutions, respectively, of the two standard, non-delayed, Hamiltonian Duffing oscillators $$x''+ ax+(-1)^nbx+x^3=0.$$ Stability changes for the delayed Duffing oscillator. Simultaneously for all sufficiently large $n\geq n_0$, we obtain local exponential stability for $(-1)^nb<0$, and exponential instability for $(-1)^nb>0$, provided that $$0 \neq (-1)^{n+1}b\,T^2< \tfrac{3}{2}\pi^2.$$ We interpret our results in terms of noninvasive delayed feedback stabilization and destabilization for large amplitude rapidly periodic solutions of the standard Duffing oscillators. We conclude with numerical illustrations of our results for small and moderate $n$ which also indicate a Neimark-Sacker torus bifurcation at the validity boundary of our theoretical results.

Equations

$$x^{\prime \prime }(t)+ax(t)+bx(t-T)+x^3(t)=0.$$

$$x^{\prime \prime }+ax+(-1{)}^{n}bx+x^3=0.$$

$$0\ne (-1{)}^{n+1}b{T}^2<\frac{3}{2}{\pi }^2.$$

$$x^{\prime \prime }+ax+bx(t-T)+x^3=0.$$

$$H=\frac{1}{2}(x^{\prime }{)}^2+\frac{1}{2}a{x}^2+\frac{1}{4}{x}^4$$

$$x(0)=A,x^{\prime }(0)=0.$$

$$x(t+p/2)=-x(t),$$

$$x^{\prime \prime }+\alpha x+x^3=u(t),$$

$$x(t+T)=(-1{)}^{n}x(t),\text{for}T=np/2,$$

$$u(t)\text{mathchar}58=\kappa \cdot (x-(-1{)}^{n}x(t-T))$$

$$x^{\prime \prime }+(\alpha -\kappa )x+(-1{)}^n\kappa x(t-T)+x^3=0$$

$$a=\alpha -\kappa ,b=(-1{)}^n\kappa .$$

$$x^{\prime }=f(x,x(t-T))$$

$$\dot{\mathbf{x}}(t)=\mathbf{f}({\mathbf{x}}_t)$$

$$\Pi ={\partial }_{\mathbf{x}}S(p,{\mathbf{x}}_0)$$

$$\dot{\mathbf{y}}(t)=f^{\prime }({\mathbf{x}}_t){\mathbf{y}}_t$$

$${\Pi }^m\text{mathchar}58$$

$${\mathbf{y}}_0\mapsto {\mathbf{y}}_{mp}$$

$$|\mu |<1$$

$$|\mu |>1$$

$$x^{\prime \prime }+(a+(-1{)}^{n}b)x+x^3=0$$

$$H=\frac{1}{2}(x^{\prime }{)}^2+\frac{1}{2}(a+(-1{)}^{n}b)x^2+\frac{1}{4}{x}^4$$

$$p=2T/n.$$

$$n>\frac{1}{\pi }T\cdot (a+(-1{)}^{n}b{)}^{1/2},$$

$$(-1{)}^{n+1}b{T}^2<\frac{3}{2}{\pi }^2.$$

$$⟺$$

$$\tilde{T}=T+kp=\nu T>0,$$

$$\underline{\nu }=\mathrm{liminf}(1+2k_n/n).$$

$$(-1{)}^{n+1}b{\underline{\nu }}^2{T}^2<\frac{3}{2}{\pi }^2.$$

$$⟺$$

Full text

[TABLE]

Fix $T>0$. We target periodic solutions $x_{n}(t)$ of small minimal periods $p_{n}=2T/n$, for integer $n\rightarrow\infty$, and with correspondingly large amplitudes. Note how $x_{n}(t)$ are also marginally stable solutions, respectively, of the two standard, non-delayed, Hamiltonian Duffing oscillators

[TABLE]

Stability changes for the delayed Duffing oscillator. Simultaneously for all sufficiently large $n\geq n_{0}$, we obtain local exponential stability for $(-1)^{n}b<0$, and exponential instability for $(-1)^{n}b>0$, provided that

[TABLE]

We interpret our results in terms of noninvasive delayed feedback stabilization and destabilization for large amplitude rapidly periodic solutions of the standard Duffing oscillators. We conclude with numerical illustrations of our results for small and moderate $n$ which also indicate a Neimark-Sacker torus bifurcation at the validity boundary of our theoretical results.

Institut für Mathematik

Freie Universität Berlin

Arnimallee 3

14195 Berlin, Germany

**

Cornell University

Ithaca,

NY 14853, USA

---

Technical University of Denmark

Nils Koppels Allé 404

2800 Kgs. Lyngby, Denmark

1. 1 Introduction and main result
2. 2 Scaling
3. 3 Floquet characteristic equation
4. 4 Wronski trace expansion
5. 5 Proof of main results
6. 6 Numerical examples

1 Introduction and main result

The Duffing oscillator [Duff1918, Kan08] is given by the special case $b=0$ of the second order pendulum equation

[TABLE]

Here we suppress time $t$ as an argument of $x=x(t)$, in absence of a delay $T>0$. For a theoretical mechanics perspective see for example [KoBr11]. For $a\geq 0=b$, the resulting integrable Hamiltonian system with energy

[TABLE]

consists of a family $x=x(t)$ of nested periodic orbits with amplitude $A>0$, i.e.

[TABLE]

The minimal period $p=p(a,A)$ is strictly decreasing to zero, for $A\nearrow\infty$, with partial derivative $p_{A}<0$. This is due to the hard spring restoring force $-ax-x^{3}$ of the Duffing oscillator. For $a<0$, in contrast, the double-well potential $\tfrac{1}{2}ax^{2}+\tfrac{1}{4}x^{4}$ of the Hamiltonian $H$ features a figure-8 pair of homoclinic lobes to the hyperbolic unstable equilibrium $x=x^{\prime}=0$, each filled with nested periodic orbits of periods $p$ bounded below. The exterior of that homoclinic pair, again, is filled with nested periodic orbits of positive Hamiltonian energy $H>0$, and with minimal periods $p\searrow 0$ for $A\nearrow\infty$. All periodic orbits with $H>0$ are odd:

[TABLE]

for all $t\in\mathbb{R}$ and any fixed real $a$. All periodic orbits are marginally stable. See also [GuHo83], [Kan08], [HiSh18] for further discussion of the classical Duffing oscillator. It is our main objective, in the present paper, to study the stabilization and destabilization of the exterior periodic orbits via a nonzero delay term $b$, in the limit of large amplitudes $A$ and for correspondingly small periods $p$.

Pyragas control [Pyr92, Pyr12, Fie&al07, Fie&al08] is a general device for noninvasive feedback control of periodic orbits. In our setting, consider the standard Duffing oscillator,

[TABLE]

with a control term $u=u(t)$. For $u\equiv 0$ we assume that the periodic orbit $x$ possesses positive energy $H>0$, so that $x$ is odd. Then (1.4) implies

[TABLE]

i.e. for any integer multiple $n=1,2,3,\ldots$ of the minimal half-period $p/2$. In particular, any delayed linear feedback control

[TABLE]

with $\kappa\in\mathbb{R}$ vanishes, i.e. becomes noninvasive, on our target periodic orbit. See [NaUe98], [Fie&al10], for this “half-period” variant of Pyragas control. See [Schn13, SchnBo16] and the references there, for many more sophisticated symmetry-related refinements. The main point is that the linear stability of the target periodic orbit $x$ of the resulting delayed Duffing oscillator

[TABLE]

may well depend on the choice of the control amplitude $\kappa\in\mathbb{R}$ and of the delay $T=np/2$. And it does! Note that (1.8) takes the form (1.1) with

[TABLE]

Our present paper in fact establishes the simultaneous stability of an unbounded sequence of rapidly oscillating periodic solutions, alternating with an unbounded sequence of rapidly oscillating periodic solutions which are unstable. See theorem 1.1 below. Here rapidly oscillating periodic solutions are defined by minimal periods $0<p<2T$.

Our results for the second order equations (1.1), (1.8) stand in marked contrast with many results in the literature on scalar delay differential equations (DDEs)

[TABLE]

for $f(0,0)=0$ and nonlinearities $f$ which are strictly monotone in the delayed variable $x(t-T)$; see [Nuss74, Wal83, MPNu13, Wal14, MP88]. Remarkably, all rapidly oscillating periodic solutions are linearly unstable in this monotone feedback setting. Only slowly oscillating periodic solutions, where consecutive zeros only occur at distances $t>T$, and therefore minimal periods exceed $2T$, may be stable. That stability also requires negative monotone feedback, i.e. nonlinearities $f$ which are strictly decreasing in the delayed variable $x(t-T)$. For the construction of a nonmonotone example with $xf(x)<0$ and infinitely many stable slowly oscillating periodic solutions see [Vas11].

In the analysis of the scalar equation (1.10) with monotone feedback, a crucial role is played by the zero number $z$: a discrete valued Lyapunov function which essentially counts the number of sign changes of solutions over time intervals of length $T$; see [MP88, MPSe96a]. This suggests that rapidly oscillating solutions of (1.10), for which necessarily $z\geq 2$, may well decay to a slowly oscillating periodic solution, for which $z\in\{0,1\}$. See [FieMP89] for such heteroclinic orbits. The opposite direction, on the other hand, is strictly forbidden by the monotone decay of the zero number.

Once the condition on monotonicity of $f$ in the delayed variable of the scalar delay equation (1.10) is removed, stable rapid oscillations have been observed; see [IvLo99, Sto08, Sto11] and the references there. The constructions of such $f$ involve smooth approximations of step functions and prove the existence of at least one linearly stable rapidly oscillating periodic solution.

Our approach preserves the monotonicity condition on $f$ but explores second order equations, instead. For second order delay equations, the existence of several slowly and rapidly oscillating solutions has been established for the delayed van-der-Pol oscillator [KiLe17]. Techniques involved a combination of the contraction mapping theorem and interval arithmetic. Stability of such solutions was not proved, but was definitely supported by numerical evidence.

For the delayed Duffing oscillator (1.1) stability of large amplitude rapidly oscillating periodic solutions has been observed numerically, and supported by formal asymptotic expansions. See [WaCha04, HaBe12, MChB15, DaShRa17]. Similar methods have been applied by [XuChu03] towards delayed feedback control of a forced van der Pol - Duffing oscillator. Neither those (or any) numerical simulations, nor the formal methods employed so far, however, amount to a mathematical proof for the coexistence of an infinity of large stable rapidly oscillating solutions. We close this gap of mathematical rigor in sections 1-5 of the present paper. For further discussion of numerical aspects which illustrate and support – but do not prove – our results, see section 6.

Asymptotic stability and instability of periodic orbits $\mathbf{x}(t)$ of general retarded functional differential equations

[TABLE]

is governed by Floquet theory. We recall [Hale77, HaleVL93, Die&al95]. For $\mathcal{X}\mathrel{\mathop{\mathchar 58\relax}}=C^{0}([-T,0],\mathbb{R}^{N})$ as a phase space for $\mathbf{x}_{t}\in\mathcal{X},\ \mathbf{x}_{t}(\vartheta)\mathrel{\mathop{\mathchar 58\relax}}=\mathbf{x}(t+\vartheta),\ -T\leq\vartheta\leq 0$, we obtain a nonlinear (local) solution semiflow $\mathbf{x}_{t}=S(t,\mathbf{x}_{0})$. Periodic orbits $\mathbf{x}(t+p)=\mathbf{x}(t)$ with minimal period $p>0$ define fixed points $\mathbf{x}_{t}=S(p,\mathbf{x}_{t})$. Let

[TABLE]

denote the linearized time-$p$ map, along our periodic orbit. Note that $\Pi$ itself is the time-$p$ map of the linearized evolution

[TABLE]

along the periodic orbit $\mathbf{x}_{t}$, where $f^{\prime}(\mathbf{x}_{t})$ denotes the Fréchet derivative of $f$ at $\mathbf{x}_{t}$. In particular, the Arzela-Ascoli theorem implies that

[TABLE]

is compact, for $mp\geqslant T$. Therefore, the spectrum of $\Pi^{m}$, and likewise of $\Pi$ itself, consists of isolated nonzero eigenvalues of finite algebraic multiplicity possibly accumulating at the spectral value [math]. The nonzero eigenvalues $\mu$ of $\Pi$ are called Floquet multipliers of $\mathbf{x}_{t}$. Note that $\mu=1$ is a (trivial) Floquet multiplier, with eigenvector $\mathbf{y}_{0}(\vartheta)\mathrel{\mathop{\mathchar 58\relax}}=\mathbf{\dot{x}}_{0}(\vartheta)$, for $-T\leq\vartheta\leq 0$. Moreover, $\mathbf{x}_{t}$ is locally asymptotically stable if

[TABLE]

holds for all Floquet multipliers, except for the trivial Floquet multiplier $\mu=1$ which is required to be algebraically simple. If $\mathbf{x}_{t}$ possesses any Floquet multiplier

[TABLE]

outside the complex unit circle, then $\mathbf{x}_{t}$ is unstable. For brevity, we use (1.15) and (1.16) as definitions of linear stability and linear instability, respectively.

To formulate our main result, fix any delay $T>0$. Let $x_{n}(t)$ denote the unique periodic orbit of the non-delayed Duffing oscillator

[TABLE]

with strictly positive Hamiltonian

[TABLE]

and with minimal period

[TABLE]

Here we have to assume that $n\in\mathbb{N}$ satisfies

[TABLE]

in case $a+(-1)^{n}b>0$. Note that (1.6), (1.19) imply that the same periodic solutions $x_{n}(t)$ also satisfy the delayed Duffing equation (1.1).

1.1 Theorem**.**

In the setting (1.17)–(1.20) assume $T>0$ and $b\neq 0$ satisfy

[TABLE]

Moreover assume that $n\geqslant n_{0}$ is chosen large enough.

Then the periodic solution $x_{n}(t)$ of the delayed Duffing equation (1.1) is

[TABLE]

in the sense of (1.15), (1.16).

Of course, condition (1.21) of the theorem is trivially satisfied in case (1.23).

The minimal periods $p=2T/n$ of our rapid periodic solutions $x_{n}(t)$ imply that the same $x_{n}(t)$ also solve the delayed Duffing equation (1.1), with the delay $T$ replaced by a new delay

[TABLE]

for any integer $k\geq-n/2$. Here we have abbreviated $\nu\mathrel{\mathop{\mathchar 58\relax}}=1+2k/n$. For fixed prescribed $T>0$, this replication in particular produces rapidly oscillating large amplitude solutions of (1.1), for a dense set of delays $\tilde{T}$ which are rational multiples of $T$. Consider any sequence $k=k_{n}$ and define

[TABLE]

For fixed odd $\nu$, for example, we may choose $k_{n}=n(\nu-1)/2$, to obtain constant $\tilde{T}=\nu T$, and $\underline{\nu}=\nu$.

1.2 Corollary**.**

In the setting (1.17)–(1.20), (1.24), (1.25), assume $b\neq 0<T$ and $\underline{\nu}>0$ satisfy

[TABLE]

Moreover assume that $n\geqslant n_{0}$ is chosen large enough.

Then the periodic solution $x_{n}(t)$ of the delayed Duffing equation (1.1) with delays $\tilde{T}$, instead of $T$, is

[TABLE]

in the sense of (1.15), (1.16).

To prove the corollary, we just note that

[TABLE]

Therefore theorem 1.1 applies to $\tilde{T},\tilde{n}$, and the corollary follows.

1.3 Corollary**.**

Consider the standard Duffing equation (1.5) with $u\equiv 0$. Fix $T>0$. Let $x_{n}(t)$ denote the unique periodic orbit with minimal period $p=2T/n$ and positive Hamiltonian. Let $n\geqslant n_{0}>\tfrac{1}{\pi}T\sqrt{\alpha_{+}}$ be chosen large enough, where $\alpha_{+}\mathrel{\mathop{\mathchar 58\relax}}=\max\{0,\alpha\}$. Assume the control amplitude $\kappa$ satisfies $-\kappa T^{2}<\tfrac{3}{2}\pi^{2}$.

Then the noninvasive delayed feedback control (1.7) of (1.5) makes $x_{n}(t)$

[TABLE]

The corollary follows from theorem 1.1, for $b=(-1)^{n}\kappa$ as in (1.9), because neither stabilization nor destabilization depends on the coefficient $a=\alpha-\kappa$ in (1.17), at all.

The remaining sections are organized as follows. In section 2 we scale the delayed Duffing oscillator (1.1) such that the amplitudes $A_{n}=x_{n}(0)$ are normalized to $1$. In particular, this introduces a small parameter

[TABLE]

which also regularizes the linearized delay equation. In section 3 we use oddness of $x_{n}(t)$ to introduce half-period Floquet multipliers $\mu$; the Floquet multipliers for the full period $p$ described in (1.12)–(1.14) above become $\mu^{2}$ in this notation. We then derive a characteristic equation (3.13) for $\mu$, of unbounded polynomial order, which involves

[TABLE]

The term $\sigma$ originates from the delay $T=np/2$ which amounts to $n$ half-periods $p/2$; see (1.19). We treat $\sigma$ as a free complex parameter in section 4. Note how $|\sigma|<1<|\mu|$ indicates linear instability. This motivates $\sigma$-uniform expansions of the $2\times 2$ Wronskian, and in particular its trace, which enters the characteristic equation. Evaluation of the expansions, in section 5, provides $\varepsilon$-expansions for the critical nontrivial half-period Floquet multiplier $\mu$ near $\mu=-1$, and proves theorem 1.1. The explicit expansions for $\mu$ itself, the standard Floquet-multiplier $\mu^{2}$, and the Floquet multiplier $\mu^{n}=1/\sigma$ at time $T=np/2$, are summarized in theorem 5.1. We conclude with some numerical illustrations, and a comparison with earlier results, in section 6.

For a nontechnical summary of our results, we refer to [Sah&al19]. Our result claims to accurately analyze and predict the in-/stability of unbounded infinities of rapidly oscillating periodic solutions in the delayed Duffing equation. Due to the delicate technical nature of our claims, we have aimed at complete mathematical proofs, which neither rely on ad hoc “approximations” nor resort to unwarranted “simplifications”.

Acknowledgement. This work originated at the International Conference on Structural Nonlinear Dynamics and Diagnosis 2018, in memoriam Ali Nayfeh, at Tangier, Morocco. We are deeply indebted to Mohamed Belhaq, Abderrahim Azouani, to all organizers, and to all helpers of this outstanding conference series. They indeed keep providing a unique platform of inspiration and highest level scientific exchange, over so many years, to the benefit of all participants. Original typesetting was patiently accomplished by Patricia Hăbăşescu. This work was partially supported by the Deutsche Forschungsgemeinschaft through SFB 910 project A4. Authors RHR and SMS gratefully acknowledge support by the National Science Foundation under grant number CMMI-1634664.

2 Scaling

In this section we rescale the integrable Duffing oscillator

[TABLE]

so that the periodic solution $x_{n}$ with minimal period $p=2T/n$, positive energy $H$, and amplitude $x_{n}(0)=A_{n}>0,\ x^{\prime}_{n}(0)=0$ is normalized to amplitude 1; see (1.17)–(1.19). With the abbreviation $\alpha=a+(-1)^{n}b$ of (1.9) and the scaling parameter $\varepsilon\mathrel{\mathop{\mathchar 58\relax}}=A_{n}^{-2}$ of (1.32), the rescaled solution

[TABLE]

obviously satisfies the rescaled Duffing equation

[TABLE]

with $\dot{}=\tfrac{d}{ds}\,$ and rescaled minimal period

[TABLE]

With the positive Hamiltonian energy

[TABLE]

of (2.2) on $x(\varepsilon,s)$, we immediately obtain the explicit solution

[TABLE]

Here $\mathrm{cn}(u,m)$ denotes the Jacobi elliptic function with parameter $m$, i.e. the primitive $\int^{u}_{0}(1-m\sin^{2}v)^{-1/2}dv\,.$ For $u=\pi/2$, we obtain the complete elliptic integral $K(k^{2})\mathrel{\mathop{\mathchar 58\relax}}=\mathrm{cn}(\pi/2,k)$. Indeed claims (2.7), (2.8) follow by integration of the ODE (2.6) with $x(0)=1$. Claim (2.9) follows from (2.5) and the scaling $A_{n}=1/\sqrt{\varepsilon}$.

Let $x_{*},\ p_{*}$ denote the resulting limits of $x(\varepsilon,\cdot),\;p(\varepsilon)$ at $\varepsilon=0$. Then (2.7), (2.8) imply

[TABLE]

In particular, the amplitude $A_{n}$ is of order $n$ or, in other words, $n$ is of order $1/\sqrt{\varepsilon}$.

As an aside, we remark that the appearance of elliptic functions in (2.3), i.e. of doubly periodic meromorphic functions in complex time, is not a surprise. Indeed, let $x(t)$ solve (2.3), this time with Dirichlet initial condition $x(0)=0$. Then $\xi(t)\mathrel{\mathop{\mathchar 58\relax}}=-ix(it)$ solves (2.3), with $\alpha$ replaced by $-\alpha$ and with the same initial condition for $\dot{\xi}(t)=\dot{x}(t)$. In particular solutions which pass through zero are doubly periodic, with one real and one imaginary period, and the periods relate the two opposite signs in the Duffing equation. See [Akh90] for more details, for example.

For later reference we remark that

[TABLE]

Indeed suppression of the index $*$, integration by parts, and $\ddot{x}=-x^{3},\ 4H\equiv 1$ yield

[TABLE]

which proves claim (2.12).

3 Floquet characteristic equation

In this section we linearize the scaled delayed Duffing oscillator

[TABLE]

which results from our original Duffing setting (1.1) via the normalization (1.32), (2.2). Here and below we write $t$ for the rescaled time $s$ again. Note the large resulting time delay $\varepsilon^{-1/2}T$. Linearization along the periodic orbits $x(\varepsilon,t)$, as in (1.11)–(1.13), leads to the linear nonautonomous delay equation

[TABLE]

Here $x^{2}(t)$ possesses minimal period $p(\varepsilon)/2$, because $x(t)$ is odd. Therefore it makes sense to study half-period Floquet multipliers, i.e. nonzero eigenvalues $\mu$ of

[TABLE]

instead of the full period $p$ in (1.12). To deal with the large time delay

[TABLE]

in (3.2), according to (2.5), we can then insert

[TABLE]

for any Floquet eigenvector $y_{t}$ of $\Pi$, and solve

[TABLE]

Here we have stealthily introduced

[TABLE]

The crucial idea in our rigorous treatment of (half-period) Floquet multipliers $\mu$, now, is to regard $\sigma$ as an independent complex parameter in the second order ODE (3.6), to be recombined with $\mu$ only later.

We rewrite (3.6) as a $2\times 2$ system

[TABLE]

and denote the linear evolution of (3.8) by the $2\times 2$ Wronski matrices $W=W(\varepsilon,\sigma,t_{1},t_{0})$ as

[TABLE]

With the subsequent abbreviation

[TABLE]

we therefore have a half-period Floquet multiplier $\mu$ if, and only if

[TABLE]

Since the trace of the right hand side of (3.8) vanishes identically, we observe

[TABLE]

Therefore (3.11) is equivalent to the Floquet characteristic equation

[TABLE]

We conclude by collecting a few special properties of $W$ at $\sigma=1$ and at $\varepsilon=0$. At parameter $\sigma=1$, we rewrite (3.1), (3.6) as the original Duffing oscillator

[TABLE]

The two columns $z_{\iota}=(y,\dot{y})^{T}$ of $W=(z_{1}\,\vdots\,z_{2})$ are given by their initial conditions $z_{1}=(1,0)^{T}$ and $z_{2}=(0,1)^{T}$, at $t=0$. We reintroduce the general amplitude $A$ of (rescaled) $x=x(\varepsilon,A,t)$ via the initial condition

[TABLE]

and denote the partial derivative $\partial_{A}x(\varepsilon,A,t)$ at (rescaled) $A=1$ by $x_{A}$. Then

[TABLE]

both solve the linearization (3.2) with the appropriate initial conditions, and we obtain the explicit expression

[TABLE]

for the Wronskian at $\sigma=1$.

At $\varepsilon=0$, the linearization (3.6) and the Floquet characteristic equation (3.13) become independent of the complex “parameter” $\sigma=(-\mu)^{-n}$, even though the delay $\varepsilon^{-1/2}T$ in (3.2) becomes unbounded. Indeed, in-/stability of $x(\varepsilon,A,t)$ at $A=1$ is only decided by $\mid\sigma\mid\leq 1$, for any $\varepsilon>0$. This justifies the notation

[TABLE]

when we address the limit $\varepsilon=0$ of original amplitudes $A_{n}=\varepsilon^{-1/2}\nearrow\infty$, below.

3.1 Proposition**.**

Assume $\sigma=1$ in the setting (3.14) – (3.19) above, and let $p=p(\varepsilon,A)$ denote the minimal period of (rescaled) $x=x(\varepsilon,A,t)$, with partial derivative $p_{A}\mathrel{\mathop{\mathchar 58\relax}}=\partial_{A}p$ at $A=1$. Then at half-period we obtain the Wronski matrix

[TABLE]

At $\,\varepsilon=0$, and independently of $\alpha$ and $\sigma$, the partial derivative $p_{A}$ of the minimal period $p_{*}=p(0,A)$ at $A=1$ is given by

[TABLE]

as detailed in (2.11). For the Wronskians at $\varepsilon=0,\ x=x_{*}(t)$ we obtain

[TABLE]

[TABLE]

Proof.

To prove claim (3.20), we just note that

[TABLE]

at half-periods, by oddness of the periodic solutions $x$; see (1.4). Differentiation of (3.24) with respect to $t$, at $t=0$ and $A=1$, and subsequent insertion into (3.19) with $t=\tfrac{1}{2}p$, proves that the second column of (3.20) follows from the initial condition on $z_{2}$. Differentiation of (3.24) with respect to $A$, at $A=1$, shows

[TABLE]

Insertion of $t=0$, of $\dot{x}=0$ at $\tfrac{1}{2}p$, and of $x_{A}=1$ at $t=0$, provides the upper left entry $x_{A}(\varepsilon,1,\tfrac{1}{2}p)=-1$ of claim (3.20). To determine the remaining lower left entry $\dot{x}_{A}$ we differentiate (3.25) with respect to $t$, at $t=0$, and obtain

[TABLE]

Here we have used the initial condition $\dot{x}_{A}(0)=0$, oddness (3.24) to replace $\ddot{x}(\tfrac{1}{2}p)=-\ddot{x}(0)$, and the ODE (3.14) with $x(0)=1$ to evaluate $\ddot{x}(0)$. We also suppressed the arguments $A=1$ and $\varepsilon$. This proves claim (3.20).

To prove claims (3.21) – (3.23), we observe that

[TABLE]

solves the pure cubic Duffing oscillator (3.14), for $\varepsilon=0$, if and only if $x(t)$ does. In particular, the scaling invariance (3.27) and the initial conditions (3.16) imply

[TABLE]

where $x_{*}(t)=x(0,1,t)$ as in (3.19). In particular, the minimal periods at $\varepsilon=0$ are given by

[TABLE]

which proves claim (3.21). For the partial derivative $x_{A}=\partial_{A}x(0,A,t)$ at $A=1$ in the Wronskian (3.18) at $\varepsilon=0$, (3.28) implies

[TABLE]

at $x=x_{*}(t)$. Insertion into (3.18) yields (3.22). Since $\textrm{det}\,W_{*}(t,0)=1$, by (3.12), this allows us to calculate

[TABLE]

Here we have used (3.20) and (3.21) to evaluate $W_{*}(\tfrac{1}{2}p_{*},0)$, and (3.22) to evaluate the unimodular inverse $W_{*}(t,0)^{-1}$. Performing the matrix multiplication proves the remaining claim (3.23), and the proposition. ∎

4 Wronski trace expansion

In the characteristic equation (3.13), the trace of the half-period Wronski matrix

[TABLE]

plays the decisive role for Floquet in-/stability. Here $W(\varepsilon,\sigma,t,0)$ accounts for the linearization (3.6), (3.8) with complex parameter $|\sigma|\leqslant 1$ along the periodic solution $x=x(\varepsilon,t)$ of (3.14) with minimal period $p=p(\varepsilon)$. Note analyticity of $x,p$ and $W$ in all variables. We collect the relevant expansions at $\sigma=1$ and at $\varepsilon=0$, as follows.

4.1 Proposition**.**

Denoting partial derivatives by indices we have

[TABLE]

at $\varepsilon=0$, for all $\sigma,t$. For the trace of $W$ we obtain the expansion

[TABLE]

with the $\sigma$-independent limiting value

[TABLE]

for the analytic function $\uptau=\uptau(\varepsilon,\sigma)$.

Proof.

With the notation $W(t)=W(\varepsilon,\sigma,t,0)$, $\ \dot{W}(t)=\partial_{t}W(\varepsilon,\sigma,t,0)$, and $W_{*}$ to indicate $\varepsilon=0$, we first recall the matrix ODE

[TABLE]

Here the $2\times 2$ matrix $\mathcal{A}=\mathcal{A}(\varepsilon,\sigma,t)$ abbreviates the right hand side of (3.8),

[TABLE]

and $W(0)=\mathrm{id}\,$. Independence of $\mathcal{A}$ on $\sigma$, at $\varepsilon=0$, proves claim (4.2); see also (3.22).

To show claim (4.3) we recall that (3.20) of proposition 3.1 shows

[TABLE]

at $\sigma=1$. For $\varepsilon=0$, where $x=x_{*}$ and $W$ are independent of $\sigma$, insertion of $t=0$ in (3.23) shows that the trace in (4.7) coincides with $\ddot{x}-x=-x^{3}-x=-2$ at $t=0$, likewise. By analyticity, this proves (4.3) with analytic $\uptau$.

To calculate $\uptau$ at $\varepsilon=0$, as required in (4.4), we first determine $\uptau(0,\sigma)$ at $\sigma=1$. By (4.3), we simply have to calculate the mixed partial derivative

[TABLE]

Differentiation of (4.5) at $\varepsilon=0,\ \sigma=1$ yields the inhomogeneous linear equation

[TABLE]

with initial condition $W_{\varepsilon\sigma}(0,1,t,0)=0$ at $t=0$. Indeed $W_{\sigma}\equiv 0\equiv\mathcal{A}_{\sigma}$ at $\varepsilon=0$, by (4.2) and (4.6). This prevents the terms $\mathcal{A}_{\varepsilon}W_{\sigma}$ and $\mathcal{A}_{\sigma}W_{\varepsilon}$ from appearing in (4.9). Solving (4.9) by variation-of-constants, suppressing indices $*$ as usual, and successively invoking (3.23), (4.6), (3.22), (2.12), we obtain

[TABLE]

This proves claim (4.4) at $\sigma=1$.

To show that the remaining values of $\uptau(0,\sigma)$, for $\sigma\neq 1$, do not depend on $\sigma$, we first differentiate (4.3) with respect to $\varepsilon$, at $\varepsilon=0$. It is sufficient to prove that $W_{\varepsilon}(0,\sigma,\tfrac{1}{2}p_{*},0)$ and $\dot{W}(0,\sigma,\tfrac{1}{2}p_{*},0)$ are affine linear in $\sigma$. The term $\dot{W}$ is actually independent of $\sigma$, by (4.2). For $W_{\varepsilon}(0,\sigma,t,0)$ we obtain, analogously to (4.9), that

[TABLE]

with $W_{\varepsilon}=0$ at $t=0$. Since the lower left entry $-(a+(-1)^{n}b\,\sigma)-6xx_{\varepsilon}$ of $\mathcal{A}_{\varepsilon}$ is the only nonzero entry of $\mathcal{A}_{\varepsilon}$, and since $\mathcal{A}$ and $W$ are independent of $\sigma$ at $\varepsilon=0$, the variations-of-constants formula shows that $W_{\varepsilon}(0,\sigma,\tfrac{1}{2}p_{*}0)$ is indeed affine linear in $\sigma$. This proves claim (4.4), in general, and the proposition. ∎

5 Proof of main results

The proof of our main results is based on the Floquet characteristic equation (3.13). Inserting the trace expansion (4.3), (4.4) this equation becomes

[TABLE]

The algebraically double trivial solution $\mu=-1$, for $\varepsilon=0$, suggests to explore an expansion

[TABLE]

of the half-period Floquet multiplier $\mu$. Since $\sigma=(-\mu)^{-n}$ also depends on $\eta$, and the power $n$ itself grows like $1/\sqrt{\varepsilon}$, the characteristic equation (5.1) is quite implicit in the scaled exponent $\eta$. One main tool in our analysis will be the nontrivial limits $\eta_{*}$ and $\sigma_{*}$ of $\eta$ and $\sigma$, respectively, for $\varepsilon\searrow 0$.

We prove our main result, theorem (1.1), in two steps. Based on expansions for $\mu$ and $\sigma=(-\mu)^{-n}$, we first consider the case of small $T>0$. In that case, theorem 5.1 provides a quantitative expansion of the leading un-/stable half-period Floquet multiplier $\mu$, in terms of $\sigma\mathrel{\mathop{\mathchar 58\relax}}=(-\mu)^{-n}$. In proposition 5.2 we then extend the resulting in-/stability to larger $T$. In fact, we assert that half-period Floquet multipliers $\mu$ cannot cause $\sigma$ to cross the unit circle $|\sigma|=1$, and thus cannot change the in-/stability result of theorem 5.1 qualitatively, as long as the crucial condition

[TABLE]

remains valid.

5.1 Theorem**.**

Uniformly for bounded $|\sigma|$, for small $0<T<T_{0}$, and for small $0\leq\varepsilon<\varepsilon_{0}\mathrel{\mathop{\mathchar 58\relax}}=\delta_{0}T_{0}^{2}$, we obtain an analytic expansion $(-\mu)^{-n}=\sigma=\sigma(\sqrt{\varepsilon},T)\in\mathbb{R}$ for the nontrivial half-period Floquet multiplier $\mu$. For $\varepsilon=0$, the $T$-expansion reads

[TABLE]

Because $|\sigma|<1$ indicates instability, and $|\sigma|>1$ indicates stability, this confirms the in-/stability claims of theorem 1.1 for small $T$ and $\varepsilon$; see (1.22), (1.23).

Proof.

We proceed in four steps. First we address the quadratic Floquet characteristic equation in the form (5.1), with analytic $\uptau=\uptau(\varepsilon,\sigma)$ from proposition 4.1. In step 2 we insert the analytic expansion of proposition 4.1 for $\sigma=(-\mu)^{-n}$, in terms of $\sqrt{\varepsilon},T,\eta$, where $\mu=-1+\sqrt{\varepsilon}\eta$ as in (5.2). This eliminates $\sigma$. In step 3, we solve the remaining equation for $\eta=\eta(\sqrt{\varepsilon},T)$, by the implicit function theorem. Insertion of $\eta_{*}=\eta(0,T)$ and expansion with respect to $T$ will complete the proof, in step 4.

Step 1: Floquet characteristic equation

Insertion of $\mu=-1+\sqrt{\varepsilon}\eta$ into (5.1) and some cancellations yield

[TABLE]

This uniformly quadratic equation for $\eta$, and analyticity of $\uptau=\uptau(\varepsilon,\sigma)$, guarantee that $\eta$ remains bounded, a priori, uniformly for bounded $|\sigma|$ and small $\varepsilon\geq 0$. Solving for $\sigma-1$, we obtain the equation

[TABLE]

which is still implicit in $\sigma$ via the trace term $\uptau=\uptau(\varepsilon,\sigma)$.

Step 2: Expansion of $\sigma$

Insertion of $\mu=-1+\sqrt{\varepsilon}\eta$ in $\sigma=(-\mu)^{-n}$ provides

[TABLE]

with the abbreviations

[TABLE]

Note that the auxiliary function $h(z)$ is entire. Replacing the integer $n$ by $2T/(p(\varepsilon)\sqrt{\varepsilon})$, as in (3.4), and expanding the logarithm, we obtain

[TABLE]

Step 3: Elimination of $\sigma$

Equating expression (5.6) with (5.7), (5.9), and cancelling out the trivial multiplier case $\eta=0$, we obtain the implicit equation

[TABLE]

for $\eta$. The somewhat messy analytic expression $\Phi$ is given explicitly as

[TABLE]

where $\sigma$ in $\uptau=\uptau(\varepsilon,\sigma)$ again has to be replaced by (5.7). For our purposes, however, it is sufficient to insert $\varepsilon=0$ and note that

[TABLE]

Therefore we can solve (5.10) for $\eta=\eta(\sqrt{\varepsilon},T)$, near $\varepsilon=T=0$, by the implicit function theorem. In particular, we obtain

[TABLE]

to leading order in $T$.

Step 4: Expansion of $\sigma_{*}$

Reinsertion of $\eta=\eta(\sqrt{\varepsilon},T)$ in (5.7)–(5.9) provides $\sigma=\sigma(\sqrt{\varepsilon},T)=\sigma(\sqrt{\varepsilon},T,\eta(\sqrt{\varepsilon},T))$. To leading order in $T$, expansion (5.13) implies

[TABLE]

Here we have substituted (4.4) for $\uptau_{*}\,$. This proves claim (5.4) and the theorem. ∎

5.2 Proposition**.**

As in theorem 1.1, (1.21), assume $0\neq(-1)^{n+1}b\,T^{2}<\tfrac{3}{2}\pi^{2}$. Then there exists a continuous function $\varepsilon_{0}=\varepsilon_{0}(T)$ such that the linear in-/stability (1.22), (1.23) does not depend on $(\varepsilon,T)$, for $0<\varepsilon<\varepsilon_{0}(T)$. Therefore in-/stability coincides with the claims of theorems 5.1 and 1.1.

Proof.

Our plan of proof is the following. It is sufficient to show the claim for $\varepsilon=0$. Extension to small $\varepsilon>0$, by the implicit function theorem applied to (5.10), then proves in-/stability as claimed in the proposition. To address $\varepsilon=0$ we will show below that $|\sigma_{*}|=1$ is impossible, under assumption (1.21), except for the simple trivial half-period Floquet multiplier $\mu=-1,\ \eta=0$. Then, the total multiplicity of solutions $\eta$ of the characteristic equation (5.11) with $|\sigma_{*}|\leq 1$, i.e. of

[TABLE]

cannot change, for increasing $T>0$, as long as (1.21) is not violated. Indeed, (5.15) is analytic in all variables, and therefore the total algebraic multiplicity of strictly unstable $|\sigma_{*}|<1$, alias $\mathrm{Re}\ \eta_{*}>0$ in (5.16) below, remains unchanged during that homotopy of $T$. For $\varepsilon=0$ and for small $0\leq\varepsilon<\varepsilon_{0}(T)$, alike, this will extend the results of theorem 5.1 from small $T>0$ to all $T$ satisfying assumption (1.21), as claimed in theorem 1.1.

To carry out this plan, consider the homotopy of $T>0$ in the unstable case $(-1)^{n}b\,T^{2}>0$ first. For $\varepsilon=0$ and small $T>0$, recall that expansion (5.14) revealed the only unstable Floquet multiplier $|\sigma_{*}|<1$ to be real, and to be given by the unique algebraically simple root $\eta=\eta_{*}$ of (5.10). In particular, that root remains simple and, for $\varepsilon=0$, extends to the full range of $T$ by our homotopy. We already mentioned how the implicit function theorem extends that instability to small $\varepsilon>0$.

In case $(-1)^{n}b\,T^{2}<0$, we address stability for small $\varepsilon>0$, indirectly. Suppose, to the contrary, that for some admissible $b,T$ there exist subsequences $n\rightarrow\infty$ with $(-1)^{n}b\,T^{2}<0$, and corresponding solutions $\eta_{n}$ of (5.10) at $\varepsilon_{n}=A^{-2}_{n}\searrow 0$ such that $|\sigma_{n}|\leq 1$ for $\sigma_{n}\mathrel{\mathop{\mathchar 58\relax}}=(-\mu_{n})^{-n}$. Since $\eta_{n}$ remain uniformly bounded, by (5.5), we can pass to convergent subsequences $\eta_{n}\rightarrow\eta_{*}\,,\ \sigma_{n}\rightarrow\sigma_{*}\,$. By continuity, $\eta_{*}$ solves (5.10) at $\varepsilon=0$ and $|\sigma_{*}|\leq 1$. For $|\sigma_{*}|<1$, this contradicts our homotopy result at $\varepsilon=0$.

After these preparations it only remains to address the stability boundary $|\sigma_{*}|=1$, for $\varepsilon=0$. In that limit, (5.7)–(5.9) imply

[TABLE]

Here $\eta_{*}$, in view of (5.5), satisfies

[TABLE]

with $\uptau_{*}$ from (4.4). Substitution of (5.16) for $\sigma_{*}$ leads to the transcendental equation

[TABLE]

with the algebraically simple trivial solution $\eta_{*}=0$.

We show, indirectly, that (5.18) cannot possess any other purely imaginary solutions $\eta_{*}=i\omega\neq 0$. Indeed any such solution would require \mathrm{sin}\,\big{(}\tfrac{2T}{p_{*}}\omega\big{)}=0\,, i.e.

[TABLE]

to annihilate the imaginary part in (5.18). To annihilate the real part then requires

[TABLE]

For even $k$ and $\omega\neq 0$, this is impossible. Hence $k$ must be odd. Substitution of (4.4) for $\uptau_{*}$ implies

[TABLE]

Insertion of (5.21) in the square of (5.19) finally requires

[TABLE]

for some odd integer $k$. But this contradicts our assumption (1.21) and the proposition is proved. ∎

It is worth noting how the first Hopf instability, at $k=1$ and $(-1)^{n+1}b\,T^{2}=\tfrac{3}{2}\pi^{2}$, determines the exponent $\eta_{*}=i\omega$ in (5.19) above. In fact, $\varepsilon^{-1/2}T=np/2$ in (3.4), odd $n$, and $y(t-\varepsilon^{-1/2}T)=-\sigma y(t)$, in (3.5), (3.7), then suggest a minimal period $q=np$ for the pair $(x(t),y(t))$. For $n>4$, this indicates a torus bifurcation at a rational rotation number, with subharmonic $1\mathrel{\mathop{\mathchar 58\relax}}n$ resonance.

More generally, our analysis (5.16)–(5.22) of nonzero purely imaginary exponents $\eta$ indicates a sequence of delay-induced torus bifurcations, which increasingly destabilize large amplitude rapidly oscillating periodic solutions $x_{n}$ of the delayed Duffing oscillator. The destabilizations originate from $n=\infty$, as $(-1)^{n+1}b\,T^{2}$ successively increases through the values $\tfrac{3}{2}(k\pi)^{2}$ for odd integer $k$ and large odd $n$. See figures 6.5, 6.6 below for illustrations of the case $k=1,\ n=33$.

Proof of theorem 1.1. With all tools at hand, we can now summarize the proof of our main result as follows. In section 2, we have rescaled the unique periodic orbits $x_{n}$ of the delayed Duffing equation (1.1), (2.1) with large amplitude $A_{n}=x_{n}(0)>0$ and rapid minimal period $p_{n}=2T/n$, to become solutions of (2.3), (2.4) with amplitude 1, small parameter $\varepsilon\mathrel{\mathop{\mathchar 58\relax}}=A^{-2}_{n}$, and rescaled minimal period $p(\varepsilon)=2TA_{n}/n$ of order 1. The advantage of (2.3), (2.4) was that the unwieldy limit of large amplitudes $A_{n}$ in (1.1), (2.1), became a regular perturbation of order $\varepsilon$. The disadvantage was the appearance of a large time delay $T/\sqrt{\varepsilon}$. In section 3 we have derived an expansion for the associated half-period Wronski matrix $W$ of the linearized rescaled delayed Duffing equation (3.8) of (2.3), along those periodic orbits; see proposition 3.1. The large rescaled time delay, however, caused the appearance of a term $\sigma\mathrel{\mathop{\mathchar 58\relax}}=(-\mu)^{-n}$ in the Floquet characteristic equation (3.13). Up to the very end, we treated $\sigma$ as just a complex coefficient in our analysis of instability, i.e. for $|\sigma|\leq 1$. Section 4 provided an expansion, in terms of $\varepsilon$ and $\sigma$, of the Wronski trace $\mathrm{tr}\,W(\varepsilon,\sigma,\tfrac{1}{2}p(\varepsilon),0)$; see proposition 4.1.

At that stage it became possible to solve the full characteristic equation (5.1), with reinserted $\sigma=(-\mu)^{-n}$, in terms of the rescaled exponent $\eta\mathrel{\mathop{\mathchar 58\relax}}=(\mu+1)/\sqrt{\varepsilon}$ for the nontrivial half-period Floquet multiplier $\mu$. In fact, the implicit function theorem provided an expansion $\eta=\eta(\sqrt{\varepsilon},T)$, although limited to small $\varepsilon,T>0$. See (5.10), (5.13). In theorem 5.1, this proved the qualitative claims of theorem 1.1 by a quantitative expansion (5.14), for small $\varepsilon,T>0$.

The full qualitative claims of theorem 1.1, for all $0\neq(-1)^{n+1}b\,T^{2}<\tfrac{3}{2}\pi^{2}$ as required in assumption (1.21), were only established in proposition 5.2. In particular it followed from the homotopy to small $T$, there, that the unstable dimensions of the original periodic orbits $x_{n}$ with $n$ large and $(-1)^{n}b>0$ are all equal to 1, given by a simple real half-period Floquet multiplier $\mu<-1$. For $(-1)^{n}b<0$ satisfying assumption (1.21), in contrast, stability prevailed. This proves the main theorem 1.1. $\bowtie$

6 Numerical examples

In this section we numerically investigate the stability and instability of the rapidly oscillating periodic solutions $x_{n}(t)$ of the delayed Duffing equation (1.1) with parameters $a=0,\,b=1$. We recall that theorem 1.1 predicts asymptotic stability, for “sufficiently large” odd $n$, and instability, for even $n$. For “sufficiently small” time delays $T>0$, more specifically, theorem 5.1 predicts an expansion

[TABLE]

of the real Floquet exponent $\tilde{\eta}=-\tfrac{1}{T}\log|\sigma|$, which determines stability; see (5.4). Let us illustrate those theoretical predictions.

To determine the amplitudes $A_{n}$ and the periodic solutions $x_{n}(t)$ in (1.17) with minimal period $p_{n}=2T/n$, we proceed as indicated in section 2. We briefly summarize these results in the original variables, prior to rescaling (2.2).

We first recall the invariant Hamiltonian (1.18) with parameters $a=0,\,b=1$ to be

[TABLE]

Solving (6.2) for $x^{\prime}\equiv dx/dt$, and separating variables, determines the minimal period $p_{n}$ of the periodic orbit $x_{n}(t)$ of amplitude $A_{n}$ to be

[TABLE]

The invariant Hamiltonian $H$ of the periodic orbit $x_{n}(t)$ can be evaluated at $t=0$, where $x_{n}=A_{n}$ and $x^{\prime}_{n}=0$, to be

[TABLE]

Replacing $H$ in (6.3) by (6.4) yields

[TABLE]

Precision values of the amplitudes $A_{n}$ are obtained by numerical solution of the implicit integral equation (6.5), for any specific value of $n=1,2,3,\ldots$ and any time delay $T$. To obtain $A_{n}$, the elliptic integral in (6.5) is numerically evaluated by the Python-based function quad. Then fsolve is called to determine the amplitude $A_{n}$ satisfying (6.5). The routine fsolve is a function wrapper around MINPACK’s hybrd and hybrj algorithms. These algorithms, in turn, are based on Powell’s hybrid method [Pow70], which combines Newton’s method and the steepest descent method. As an initial guess for $A_{n}$ in fsolve, for any chosen values of $n$ and $T$, we use the approximation in [DaShRa17], Eq. 4, for the exact amplitude. To double check, we have also solved (6.5) for $A_{n}$ by explicit inversion of the series expansion for the elliptic integral, to degree 9, using the symbolic mathematica package and precision evaluation of the Gamma function value $\Gamma(\tfrac{1}{4})$. This provided the approximations given below.

To determine the exact solutions $x_{n}(t)$ of (1.17) with minimal periods $p_{n}=2T/n$ and amplitude $A_{n}$ we recall section 2. Indeed the exact solutions $x_{n}(t)$ are expressed by elliptic integrals similar to (6.5). This leads to the Jacobi elliptic cosine function cn,

[TABLE]

Here $A_{n}\,,\,\omega_{n}$ , and $0<m_{n}<1$ are the amplitude, angular frequency, and Jacobi elliptic modulus, respectively; see [Akh90]. The three parameters are related to each other through the equations

[TABLE]

see also [Rand94], for example. Here $A_{n}>\sqrt{2}$ for odd $n$, to ensure $H>0$. In particular the single parameter $A_{n}$ , as determined above, fully describes the exact solution (6.6) with prescribed minimal period $p_{n}=2T/n$. The amplitudes $A_{n}$, derived from (6.3) numerically, can therefore be substituted into (6.7) to obtain the reference periodic solutions $x_{n}(t)$ of (1.1).

Next, we numerically integrate the initial value problem for the delayed Duffing equation (1.1) using the dde23 package. For the numerical integrations, we use Pydelay [Flun11], which is a Python library for DDEs. The code of dde23 is based on the Bogacki-Shampine method [BoSh89] which, in turn, implements the 2(3) Runge-Kutta method. All plots in the present work fix the maximal step size at $\Delta t=10^{-4}$.

As initial conditions we consider the Jacobi elliptic history functions

[TABLE]

of amplitude $A$, for $-T\leq t\leq 0$. See (1.11) for the notation $x_{t}$ . Here $\omega$, $m$ are again defined via (6.7), once the initial amplitude $A$ is chosen. In particular, initial amplitudes $A$ close to the amplitudes $A_{n}$ of the periodic solutions $x_{n}(t)$ indicate initial histories $x_{0}$ close to the periodic histories $(x_{n})_{0}$ in function space $\mathcal{X}=C^{0}([-T,0],\mathbb{R}^{2})$.

For delay $T=0.6$, figure 6.1 compares a numerical solution (blue) of the delayed Duffing equation (1.1) with the reference periodic solution $x_{n}(t)$ (red) for $n=1$. The figure shows time history (a) and phase plane (b). The green curve denotes the initial history function (6.8) with initial amplitude $A=3.7$. The black curve denotes the final state of the history function. The simulated solution (blue) approaches the reference periodic solution (red) for large times $t>20$. Locally, but neither for small $n=1$ nor for the large initial deviations $|A-A_{n}|$ tested here, this is predicted by asymptotic stability of the periodic solution $x_{n}(t)$ for odd $n$, according to theorem 1.1.

Again for delay $T=0.6$, figure 6.2 shows the time histories of two numerical solutions of the delayed Duffing equation (1.1) with initial history functions (6.8) and initial amplitudes $A=1.42$ (blue) and $A=12.29$ (light purple), respectively. The amplitudes $A_{n}$ of the reference periodic solutions $x_{n}(t)$ for $n=1$ (red) and $n=2$ (teal), respectively, are $A_{1}=6.29721145...$ and $A_{2}=12.30144591494...\,$. Figure 6.2 shows how the simulated (blue, light purple) solutions for both initial amplitudes approach the same reference (red) periodic solution $x_{1}(t)$. Also note how the simulated solution with initial amplitude $A=12.29$, quite close to the periodic amplitude $A_{2}=12.30144591494...$, actually diverges from the reference (teal) periodic solution $x_{2}(t)$. Again, this confirms the asymptotic stability of the periodic solution $x_{n}(t)$ for $n=1$, and instability for $n=2$, as predicted by theorem 1.1. The global feature of heteroclinicity from $x_{2}$ to $x_{1}$, manifested by the light purple orbit, is beyond our present scope, of course.

To test the expansion (6.1) of theorem 5.1 for the Floquet exponent $\tilde{\eta}$ of periodic solutions $x_{n}(t)$, we track the Hamiltonian (6.2) numerically. See figures 6.3 and 6.4 for illustrations, as detailed below. Tracking the Hamiltonian eliminates the lack of convergence in phase, which is due to the trivial Floquet exponent $\tilde{\eta}=0$. Indeed, let $H_{n}>0$ denote the time-independent Hamiltonian on $x_{n}(t)$; see (6.4). Then $|H(t)-H_{n}(t)|$ indicates the distance of our numerical solution $x(t)$ for the delayed Duffing equation (1.1) from the reference periodic orbit $\{x_{n}(t)\,|\,t\in\mathbb{R}\}$, as a set, rather than the distance from any particular point $x_{n}(t)$ on that orbit.

We track the time-dependent Hamiltonian $H=H(t)$ as it exponentially converges to (orange curves), or diverges from (teal curves), the stationary limit $H_{n}$. According to theorem 5.1, this occurs for odd $n$ (orange) and even $n$ (teal), respectively. Let us be a little more specific. For odd $n$ (orange), we start with initial amplitudes $A$ slightly below $A_{n}$ and observe convergence to $H_{n}$. For even $n$ (teal), we start with initial amplitudes $A$ slightly below $A_{n}$ and observe convergence to $H_{n-1}$.

The slope of $\log|H(t)-H_{n}|$, asymptotically with respect to time $t\rightarrow\pm\infty$, then coincides with the Floquet exponent $\tilde{\eta}$. This determines the instability or stability of the periodic solution $x_{n}(t)$, depending on the positive or negative sign of the slope.

Figures 6.3 and 6.4 show examples of the simulated Hamiltonian $H(t)$ and of the constant reference Hamiltonian $H_{n}$ . The figures confirm that periodic solutions $x_{n}(t)$ are stable for odd $n$ (orange curves), while even $n$ are unstable (teal curves).

In figure 6.3, with delay $T=0.3$, precision amplitudes $A_{n}$ of the stable orange reference periodic solutions are $A_{1}=12.41931822569...$ and $A_{11}=135.97083402978303460...$; the amplitudes of the unstable teal periodic solutions are $A_{2}=24.69151341060282...$ and A_{12}=148.32106281755626611...\. Initial conditions are $A=4.68$ and $A=130$ for the orange curves converging to $H_{1}$ and $H_{11}$ , respectively. The teal curves diverging from $H_{2}$ and $H_{12}$ start from $A=24.68$ near $A_{2}$, and $A=147$ near $A_{12}$ .

In figure 6.4, with delay $T=0.9$, amplitudes $A_{n}$ of the stable orange reference periodic solutions are $A_{27}=111.25102887868052589...$ and $A_{51}=210.13193020360773942...$; amplitudes of the unstable teal periodic solutions are $A_{28}=115.35833191723956861...$ and A_{52}=214.24522922435665376...\. Initial conditions are $A=110.1$ and $A=209.1$ for the orange curves converging to $H_{27}$ and $H_{51}$ , respectively. The teal curves diverging from $H_{28}$ and $H_{52}$ start from $A=115.2$ near $A_{28}$ , and $A=214$ near $A_{52}$ .

Lower plateaus in the logarithmic plots indicate residual relative numerical errors of our numerical simulations. The local relative error tolerance of dde23 is $10^{-3}$. All simulations support the theoretically predicted Floquet exponent $\tilde{\eta}=(-1)^{n}T/3$ of (6.1), which corresponds to the slopes $\sim\pm 0.1$ for $T=0.3$, in figure 6.3 (c), (d), and slopes $\sim\pm 0.3$ for $T=0.9$ in figure 6.4 (c), (d). Slopes were determined by least square fits (gray lines). Note how the slopes only depend on the even/odd parity, but not on the value, of $n$, asymptotically for large $n$. Given that our original expansion (5.4) was limited to “sufficiently small” $T$ and “large enough” $n$, we are rather surprised at such quantitative agreement far from those limits.

Our final figures are testing for the conjectural torus bifurcation at the critical boundary

[TABLE]

of assumption (1.21) in theorem 1.1, for odd $n$ and $b=1$. We check for oscillatory $H(t)$ at $T=T_{\textsf{crit}}+0.1\,$.

In figure 6.5 we plot the relative deviations $(H(t)-H_{n})/H_{n}$ for $n=33$. On the top (orange), we start at an initial amplitude of $A=31.1$, slightly above the amplitude $A_{33}=31.021414799836585...$ of the reference rapid periodic solution $x_{33}(t)$. We clearly observe a loss of stability of the solution $x_{33}(t)$, which we asserted to be stable for $T<T_{\textsf{crit}}$ and large enough $n$. In fact $H(t)-H_{33}$ increases to an asymptotically periodic, sinusoidal oscillation. This suggests a supercritical Neimark-Sacker bifurcation [Ioo79, KuSa08] of $x_{33}(t)$ to a stable 2-torus $\mathbb{T}^{2}$, near $T=T_{\textsf{crit}}$ . Our remark following the proof of proposition 5.2 predicts that the bifurcation occurs near $1\mathrel{\mathop{\mathchar 58\relax}}n=1\mathrel{\mathop{\mathchar 58\relax}}33$ resonance. This may result in an asymptotic oscillation of $x(t)$ itself which is still $1\mathrel{\mathop{\mathchar 58\relax}}33$ subharmonic, or possibly quasiperiodic.

On the bottom (teal), we start at an initial amplitude $A=31.9$ slightly below the reference amplitude $A_{34}=31.91443613945749...\,$. Real instability of $x_{34}(t)$ persists to dominate, and we observe an asymptotic decay to the same sinusoidal periodic oscillation of $H(t)-H_{33}$ , as in the previous case. This further attests to the presence of a stable invariant 2-torus $\mathbb{T}^{2}$, which causes the asymptotically $1\mathrel{\mathop{\mathchar 58\relax}}33$ resonant subharmonic, or possibly quasiperiodic, oscillation of $x(t)$, further examined in figure 6.6.

For post-transient times $300\leq t\leq 311$, we examine a close-up on the sinusoidal oscillations of $-0.04\leq(H(t)-H_{33}))/H_{33}\leq+0.04$ and the subharmonic, or possibly quasiperiodic, fluctuation of $x(t)$ in the peak region $30.7\leq x\leq 31.35$. See the top and bottom rows (a), (b) and (c), (d) of figure 6.6, respectively. The orange graphs, in the left column, and the teal graphs, in the right column, refer to the same initial conditions as in figure 6.5. The top row clearly indicates convergence of both solutions to the same invariant 2-torus $\mathbb{T}^{2}$, with identical sinusoidal oscillations of $H(t)-H_{33}$ up to a phase shift. In particular the asymptotic periods $\sim 2T$ coincide, right and left. The sinusoidal character of $H(t)-H_{33}$ indicates that our delay parameter $T$ is close to the actual bifurcation point, where the invariant 2-torus $\mathbb{T}^{2}$ bifurcates from the destabilizing rapidly periodic reference solution $x_{33}(t)$.

The bottom row of figure 6.6 shows a slow sinusoidal fluctuation, over slow periods $\sim 2T$, of the amplitudes of the rapid oscillations of $x(t)$ of minimal “periods” near 2T/n=0.23925...\. This indicates the $1\mathrel{\mathop{\mathchar 58\relax}}n=1\mathrel{\mathop{\mathchar 58\relax}}33$ subharmonic, or possibly quasiperiodic, flow on the invariant 2-torus $\mathbb{T}^{2}$, and agrees well with our remark following the proof of proposition 5.2.