Delay-induced homoclinic bifurcations in modified gradient bistable systems and their relevance to optimisation

TL;DR

This paper investigates how delay-induced homoclinic bifurcations in modified gradient systems can reorganize phase space, potentially aiding in global optimization by enabling trajectories to explore multiple local minima.

Contribution

It introduces a class of delayed gradient systems with two-well potentials and demonstrates how delay can induce bifurcations that reorganize attractors, offering insights into delay's role in optimization.

Findings

Delay induces a chain of homoclinic bifurcations.

Bifurcations eliminate local attractors.

Manifolds reorganize, removing barriers between minima.

Abstract

Nonlinear dynamical systems with time delay are abundant in applications, but are notoriously difficult to analyse and predict because delay-induced effects strongly depend on the form of the nonlinearities involved, and on the exact way the delay enters the system. We consider a special class of nonlinear systems with delay obtained by taking a gradient dynamical system with a two-well "potential" function and replacing the argument of the right-hand side function with its delayed version. This choice of the system is motivated by the relative ease of its graphical interpretation, and by its relevance to a recent approach to use delay in finding the global minimum of a multi-well function. Here, the simplest type of such systems is explored, for which we hypothesise and verify the possibility to qualitatively predict the delay-induced effects, such as a chain of homoclinic bifurcations one by one eliminating local attractors and enabling the phase trajectory to spontaneously visit vicinities of all local minima. The key phenomenon here is delay-induced reorganisation of manifolds, which cease to serve as barriers between the local minima after homoclinic bifurcations. Despite the general scenario being quite universal in two-well potentials, the homoclinic bifurcation comes in various versions depending on the fine features of the potential. Our results are a pre-requisite for understanding general highly nonlinear multistable systems with delay. They also reveal the mechanisms behind the possible role of delay in optimisation.

Figures

Tables

Table 1. Table 1: Features of the minima and a maximum of the potential V 𝑉 V specified by ( IV.1 ), and the points of homoclinic bifurcations in ( 2 ), ( IV.1 ) as τ 𝜏 \tau grows. The Table provides the locations of the minima x 1 , 2 min superscript subscript 𝑥 1 2 min x_{1,2}^{\textrm{min}} and the maximum x max superscript 𝑥 max x^{\textrm{max}} , their depths V 1 , 2 min superscript subscript 𝑉 1 2 min V_{1,2}^{\textrm{min}} and V max superscript 𝑉 max V^{\textrm{max}} , Jacobians J 1 , 2 min superscript subscript 𝐽 1 2 min J_{1,2}^{\textrm{min}} and J max superscript 𝐽 max J^{\textrm{max}} , and the points of the first AH bifurcations τ 1 , 2 AH1 superscript subscript 𝜏 1 2 AH1 \tau_{1,2}^{\textrm{AH1}} and τ AH1 superscript 𝜏 AH1 \tau^{\textrm{AH1}} . For x max superscript 𝑥 max x^{\textrm{max}} , τ σ superscript 𝜏 𝜎 \tau^{\sigma} is the value of τ 𝜏 \tau at which its saddle quantity switches from negative to positive. All homoclinic bifurcations here are homoclinic loops of the saddle-focus fixed point at x max superscript 𝑥 max x^{\textrm{max}} , which are safe since all values of τ ​ J max 𝜏 superscript 𝐽 max \tau J^{\mathrm{max}} are below 1.975 1.975 1.975 .

Features of the potential $V$ in (IV.1)
$x_1^{\text{min}}$$=$$-2.77168$ $x_2^{\text{min}}$$=$$0.864388$ $x^{\text{max}}$$=$$-0.308654$ $V_1^{\text{min}}$$=$-1.7557 $V_2^{\text{min}}$$=$$-1.2417$ $V^{\text{max}}$$=$$-0.9207$ $J_1^{\text{min}}$$=$$-1.30015$ $J_2^{\text{min}}$$=$-2.19511 $J^{\text{max}}$$=$$0.490368$ ${\tau }_1^{\text{AH1}}$$=$$1.208$ ${\tau }_2^{\text{AH1}}$$=$$0.7155$ ${\tau }^{\text{AH1}}$$=$$9.61$ flattest, broadest, lowest ${\tau }^{\sigma }$$=$4.0276
Homoclinic bifurcation points in (2), (IV.1)
1st homoclinic (safe small loop of $x^{\text{max}}$): $\tau$$\approx$$1.53$, $\tau J^{\max }$$=$$0.75$
2nd homoclinic (safe large loop of $x^{\text{max}}$): $\tau$$\approx$$2.4307$, $\tau J^{\max }$$=$$1.19$
3rd homoclinic (safe small loop of $x^{\text{max}}$): $\tau$$\approx$$2.4499$, $\tau J^{\max }$$=$$1.201$

Equations

$$\dot{x}=-\nabla V(x),$$

$$\dot{x}=f(x_{\tau }),f(z)=-\frac{\mathrm{d}V(z)}{\mathrm{d}z},$$

$$\frac{\mathrm{d}x}{\mathrm{d}t}=g(x(t-1)),g(z)=-\frac{\mathrm{d}{V}_1(z)}{\mathrm{d}z}.$$

$$\frac{\mathrm{d}x}{\mathrm{d}t}=ah(x(t-1)),a>0,h(z)=-\frac{\mathrm{d}{V}_2(z)}{\mathrm{d}z}.$$

$$h(z)=\left(-\frac{1}{1+{\mathrm{e}}^{-\frac{z}{2}}}+{\mathrm{e}}^{-\frac{z}{2}}\right)\left(2-{\mathrm{e}}^{-\frac{z}{2}}\right)$$

$${\lambda }_k=\frac{W_k(J\tau )}{\tau }=\frac{J{W}_k(z)}{z},$$

$${\tau }_{\text{AH1}}=\frac{\pi }{2|J|},$$

$${\tau }_{\text{AH1}}=\frac{3\pi }{2J},$$

$$V(x)$$

$$f(x)$$

$$V(x)$$

$$f(x)$$

Full text

delay differential equation, homoclinic bifurcation, global bifurcation, predictability, optimisation

††preprint: AIP/123-QED

Can one predict, without resorting to numerical analysis, the behaviour of a nonlinear system with time delay described by a delay differential equation (DDE) as the delay is gradually increased? Generally not, as it is highly sensitive both to the form of nonlinearities, and to how the delay is introduced. Moreover, unlike in ordinary differential equations, their phase space is infinite-dimensional, and it is generally impossible to reverse the time, so nonlinear DDEs represent a challenge both for analytical, and for numerical treatment. However, here we construct a special form of highly nonlinear DDEs, in which one can qualitatively predict a sequence of bifurcations as the delay grows. Namely, we slightly modify the most basic setting for optimisation, when the parameter to be optimised decreases as the negative of the gradient of some multi-well “cost" function (called “potential energy landscape function" in physics problems), by delaying its argument. After overviewing some earlier rigorous results available for simpler DDEs, and interpreting them in terms of the potential function, we use these to predict phenomena in more complex DDEs with two-well potentials, which could be extended to multi-well potentials in the future. We hypothesise and verify some universality in the sequence of global homoclinic bifurcations, and also establish how different forms of these bifurcations are realised with different local features of the potential. Delay-induced global bifurcations can be the means to remove the barriers between the local minima of the cost function, and thus to allow the phase trajectory to approach all minima, similarly to what occurs in a famous optimisation method simulated annealing thanks to random forces. This effect seems promising for optimisation where the delay could replace random forces. Since the barriers are embodied in the manifolds of saddle points or saddle cycles, their reorganisation via homoclinic bifurcations is key to understanding how the barriers disappear as the delay grows. We explain rearrangement of manifolds to shed light on this mechanism.

I Introduction

Differential equations with time delay represent a special class of dynamical systems, which are routinely used to model the behaviour of both natural systems and artificial devices alongside with ordinary differential equations (ODEs). Delay equations have been introduced in 1940s as models of population dynamics [1, 2] and later of other biological phenomena [3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. In some models, the delay appears in the term(s) added to the components(s) of the original ODE [13, 14, 15]. Other models contain combinations of delayed and non-delayed terms from the outset [16, 2, 17].

Delay differential equations (DDEs) present a considerably greater challenge for the analysis than ODEs because their state is represented by a (vector-) function on an interval, rather than by a finite-dimensional vector, and hence the dimension of their phase space is infinitely large. Also, time reversal is generally not permitted, which complicates the detection of unstable objects. Moreover, the effects induced by delay greatly depend on a particular form of the DDE under study and on the exact way the delay is introduced. The latter makes it hardly possible to predict, before resorting to numerical analysis and actually observing the behaviour, the dynamics of even a scalar equation with a single delay \tau$$\geq[math], i.e. of \dot{x}$$=$$f(x,x_{\tau}) with $x,f$ $\in$ $\mathbb{R}$ and x_{\tau}$$=$$x(t-\tau) for an arbitrary $f$.

However, by extending the results from the qualitative theory of ODEs [18, 19, 20, 21] it has been possible to qualitatively predict certain phenomena in special cases. Namely, some predictions were made for DDEs reducible to the form \dot{x}$$=$$f(x_{\tau})$$-$$g(x), often with g(x)$$=$$\lambda x (\lambda$$\in$$\mathbb{R}) [22, 23, 17, 15, 24, 25], or to an even simpler form \dot{x}$$=$$f(x_{\tau}), for some special forms of $f$ [26]. With this, it is usually impossible to give accurate analytical predictions of the behaviour of general non-linear DDEs, and their studies heavily rely on numerical tools. Even more challenging in DDEs is the detection and interpretation of homoclinic and heteroclinic, i.e. global, bifurcations. The reason is that they are based on invariant manifolds, which are usually not confined to a small volume of the phase space (i.e. not localised) and can be one-, two-, many-, or infinite-dimensional. The tools available to date can handle only unstable invariant manifolds in DDEs, and examples considered involve one-dimensional [27] and two-dimensional [28, 29] manifolds.

Here, we explore a special class of nonlinear DDEs with bistability, for which we hypothesise and verify the possibility to predict qualitatively how the behaviour changes with the increase of delay. General bistable dynamical systems with delay are popular models in a range of areas, such as optical systems [30], laser systems [31], neural networks, atmospheric physics [32], and ice dynamics [33], so our work will provide an additional insight into typical bifurcations determining their behaviour.

The nonlinear DDE to explore is constructed by modifying a nonlinear ODE with arguably the most predictable behaviour, which is the gradient dynamical system of the form

[TABLE]

where {x}(t)$$\in$$\mathbb{R} is the state, $t$ is the time, \dot{{x}}$$=$$\frac{\textrm{d}{x}}{\textrm{d}t}, and V({x})$$:$$\ \mathbb{R}$$\to$$\mathbb{R} is the potential energy function at least twice continuously differentiable [34]. Also, $\nabla$ is the gradient operator, which for a scalar $x$ is equivalent to $\frac{\partial}{\partial x}$. This system models the behaviour of a particle in a potential energy landscape $V$ immersed in viscous fluid, so that the particle cannot oscillate. Assuming that $V$ is a multi-well function with several maxima separating them, the system converges to one of the local minima in a non-oscillatory manner, and the choice of the minimum depends on the initial conditions. System (1) represents the most basic setting for optimisation problem, which in practical applications is posed for {x}(t)$$\in$$\mathbb{R}^{N}, V({x})$$:$$\ \mathbb{R}^{N}$$\to$$\mathbb{R} with N$$\geq$$1, where $V$ is the multi-well “cost" function, and $x$ is the vector of parameters in need of optimisation [35, 36]. Solving this ODE can only deliver a local minimum, so to find the global one this setting is usually extended to enable the particle to overcome the barriers between the minima. Most popular extensions rely on incorporating random force that makes the particle visit various regions of the landscape, including the vicinity of the global minimum[37, 38, 39, 40].

We delay the argument of the right-hand side of (1) by some amount \tau$$\geq[math] and thus obtain a DDE to study,

[TABLE]

where $x$, $f$, $V$, $z$ $\in$ $\mathbb{R}$, x_{\tau}$$=$$x(t-\tau). This form of (2) allows us to make some rough qualitative predictions about its solutions based on the knowledge of general bifurcation theory and of general delay-induced effects in dynamical systems. Namely, in [41] it was hypothesised and numerically demonstrated that an increase of $\tau$ could eliminate one by one all local attractors via global homoclinic and possibly heteroclinic bifurcations, and give rise to a single large attractor embracing all local minima of $V$. Thus the phase trajectory could be forced to eventually approach all minima of $V$, including the global one, regardless of the initial conditions. This means that the delay could roughly mimic the effect from adding large noise to (1), as in global optimisation with simulated annealing [38]. However, here it would be achieved in a fully deterministic manner and through an entirely different mechanism.

Specifically, elimination of local attractors is only possible if their individual basins of attraction cease to be separated by boundaries formed by manifolds of saddle fixed points or of saddle cycles. In the context of optimisation, these manifolds serve as the barriers between the minima of $V$. After the homoclinic bifurcations these manifolds do not disappear, but reorganise in such a way, that they cease being the barriers. Thus, homoclinic bifurcations could effectively break down the barriers between the local minima. This would provide a mechanism, alternative to random noise, for overcoming barriers as required in optimisation. With this, understanding of the way the manifolds reorganise as the delay grows would provide the key to understanding if and how optimisation by delay could work.

As a starting point and a pre-requisite for understanding the phenomena in delay systems with general multi-well $V$ depending on one or many variables, here we focus on system (2) with two-well smooth landscape functions $V$, such that V(z)$$\rightarrow$$\infty as |z|$$\rightarrow$$\infty and $V$, $z$ $\in$ $\mathbb{R}$.

In Section II we overview mathematical theorems related to some simpler delay systems, interpret and illustrate them for the general reader, and put into context of a DDE (2) with energy $V$. In Section III we analyse stability of the fixed points of (2) and its relevance to homoclinic bifurcations. In Section IV, we combine rigorous theoretical and quantitative results overviewed in Sections II and III with the powerful apparatus of the qualitative theory of ODEs [18, 19, 20, 21], which has proved to be successfully applicable to DDEs as well, in order to predict, reveal and explain the rather intricate bifurcation phenomena induced by delay in (2) with a double-well $V$. Here, we illustrate these phenomena with specific examples of $V$ in (2), demonstrate homoclinic bifurcations of various types, and reveal distinctions and universality in the delay-induced behaviour of such systems. In Section V we put the delay-induced phenomena in the context of optimisation problem and demonstrate how optimisation could work with two local minima if one uses the delay. In Section VI we discuss the results obtained.

II Delay-induced behaviour in simple systems

In order to predict some phenomena that might be caused by delay in systems (2) with two- and multi-well landscapes $V$, it is important to know about the phenomena occurring in systems with $V$ having only one well. Note that multi-well functions $V$ can be obtained by gluing together the segments of single-well $V$’s and smoothing out all the joints. Thus, bifurcations involving objects localised within a single minimum could be expected in those with more minima, too. This Section overviews the results available for functions $V$ with a single minimum and up to one maximum.

Most theorems fomulated for (2) assume that \tau$$=$$1 and $f$ is amendable [42, 43, 44, 45, 46, 47, 48, 49, 50, 26], i.e. the equation reads

[TABLE]

However, introducing t$$=$$s\tau and x(t)$$=$$y(s) reduces (2) to \frac{\mathrm{d}y}{\mathrm{d}s}$$=$$\tau f(y(s-1)), i.e. to (3) with g(z)$$=$$\tau f(z) and the appropriate change in notations. Thus, the results valid for (3) can be easily adapted to (2) by using the fact that the increase of $\tau$ in (2) is equivalent to sharpening and deepening the wells of $V_{1}$ in (3), as illustrated in Fig. 1.

Note, that a typical global bifurcation in the DDEs being considered is associated with the formation of a homoclinic orbit, also called a homoclinic loop, starting and finishing at the same saddle fixed point. For a periodic orbit born from this loop, as the system approaches the bifurcation point, the period tends to become infinitely large. In ODEs there exists another global bifurcation with the same feature of the periodic orbit involved, called saddle-node homoclinic bifurcation and otherwise known as SNIPER, which was first described by Andronov in the 1940s (see 30 of [[51](#bib.bib51)]). A SNIPER bifurcation occurs when, as the control parameter changes monotonously, two fixed points (a saddle and a node) move towards each other while staying on the same limit cycle, and disappear in a saddle-node bifurcation. After this bifurcation, the whole of the limit cycle becomes the attractor. In the setting we consider, the only control parameter is the time delay\tau$, and in DDEs with constant delays the number or locations of the fixed points are not affected by the numerical value of the delay. Thus, the increase of the delay alone cannot lead to the movements or the disappearance of the fixed points, and the SNIPER bifurcation cannot occur in the situations we consider in this paper.

The idea of using delay for optimisation is inspired by the knowledge of typical behaviours of (3) with relatively simple $g$. Although the key relevant theorems are available from highly specialised literature, they have not been previously presented in a form accessible to a more general reader, and not interpreted in the context of the landscape, which we need before explaining how we combine them in order to predict the behaviour of (2) with complex-shaped $V$. In sub-sections below, three prominent special cases of $g$ in (3) are considered, namely, monotonically decreasing with a single zero-crossing, and making two zero-crossings with a single maximum or a single minimum. We illustrate these cases with specially constructed examples.

II.1 Existence of a periodic orbit

In [42, 43], Eq. (3) was considered with $g(z)$ continuously differentiable and monotonically decreasing on some open neighbourhood of z$$=[math], crossing zero at z$$=[math], and in addition satisfying zg(z)$$<[math] for all z$$\neq[math]. I.e. $g(z)$ should be strictly positive for negative $z$ and strictly negative for positive $z$. Therefore, the fixed point of (3) is at x$$=[math].

In addition, $g(z)$ should either be bounded from below for all $z$ (Fig. 2(b)), or satisfy |g(z)|$$\leq$$A if |z|$$\leq$$A, where $A$ is some positive constant (Fig. 2(d)). If |g^{\prime}(0)|$$>$$\frac{\pi}{2}, then (3) has a non-zero periodic solution. In [45] the stability theorem for this periodic solution is proved. From the viewpoint of the landscape, to enable a periodic solution, the single well of the twice differentiable $V_{1}$ should be sufficiently sharp at the bottom (Fig. 2(a), (c)).

Note, that the quoted theorems require that $g(z)$ crosses zero at z$$=[math]. However, obviously, the same results can be adapted to $g(z)$ crossing zero at any $z$ by making an appropriate change of variable. Thus, while the shape of $g$ contains all the necessary information for the predictions, its localisation on $z$-axis is not essential.

II.2 Existence of a homoclinic orbit to a fixed point

Equation (3) with a non-monotonic $g$ was studied in [52, 53] but rewritten as

[TABLE]

Here, $h(z)$ is required to cross zero at z$$=[math] from below to above, and at z$$=$$-|b|, b$$\in$$\mathbb{R}, from above to below. Parameter $a$ controls the sharpness of the two extrema of $V_{2}$. It has been shown that at certain $a$, in (4) there exists a stable periodic orbit around x_{2}$$=$$-|b| (red circle in Fig. 3). As $a$ increases, this orbit grows in size and, under some additional quantitative conditions on $h(z)$, at $a$ equal to some critical value $a^{*}$ can clash with the saddle fixed point at x_{1}$$=[math] (green circle in Fig. 3) and form a homoclinic loop. At a$$>$$a^{*} this orbit no longer exists.

This scenario is illustrated in Fig. 3 for a function $h(z)$, which has the required properties on the domain considered (see Fig. 3(a) and caption). With this $h$, Eq. (4) has two relevant fixed points: x_{2}$$=$$-2.5 (red circle) and x_{1}$$=[math] (green circle) corresponding to the local minimum and maximum of $V_{2}$, respectively (Fig. 3(a)–(b)). At small $a$, the point $x_{2}$ is stable, but it loses stability at a$$\approx$$1.396 at which $|ah^{\prime}(-2.5)|=\frac{\pi}{2}$. The point $x_{1}$ is unstable for any a$$>[math]. At a$$=$$2.24 there exists a stable periodic orbit around $x_{2}$ shown in Fig. 3(f) in projection on the plane $(x(t),x(t-1))$, and the respective solution $x(t)$ is given in (e). At a$$=$$a^{*}$$\approx$$2.3765 a homoclinic orbit is formed as the periodic orbit collides with $x_{1}$, as shown in (d). In (c) the respective homoclinic solution $x(t)$ is shown, which departs from x$$=[math] starting from some time instant in the past and approaches x$$=[math] as t$$\rightarrow$$\infty 111In ODEs, a homoclinic solution would be called doubly asymptotic to the fixed point (here x$$=[math]), which means that it approaches this point when $t$ tends not only to positive, but also to negative, infinity, given that time can be reversed. However, in DDEs time cannot be generally reversed, so the term doubly asymptotic might not be fully appropriate here..

One can interpret these events in the context of Shilnikov’s theorem about the birth of a periodic orbit from the breakdown of the homoclinic loop of a saddle-focus fixed point [55, 56]. For $h(z)$ in Fig. 3(b) for a range of values of a$$>[math], the fixed point x_{1}$$=[math] of (4) is a saddle focus with a single real positive eigenvalue $\lambda_{1}$, whereas all of its other eigenvalues are complex-conjugate with negative real parts (see Section III). According to Shilnikov’s theorem for ODEs and a relevant study of a special form of the DDE (3) [50, 52], the breakdown of a homoclinic loop of a saddle-focus gives birth to a periodic orbit, if at the instant of homoclinic bifurcation the loop is “safe". The latter means that the saddle quantity of this saddle fixed point is negative, i.e. \sigma$$=$$\lambda_{1}+Re(\lambda_{2,3})$$<[math], where $\lambda_{2,3}$ are eigenvalues with the negative real parts closest to zero.

From this viewpoint, in (4) the homoclinic loop shown in Fig. 3(d) is formed as the parameter $a$ decreases to its critical value $a^{*}$ from above. As $a$ is further decreased below $a^{*}$, from this homoclinic orbit a periodic orbit is born (Fig. 3(f)) provided the negativity of $\sigma$ of the fixed point x_{1}$$=[math]. If \sigma$$>[math], then in ODEs, according to Shilnikov, the breakdown of the homoclinic loop should produce chaos rather than periodic behaviour.

II.3 Existence of chaos from the homoclinic orbit to a saddle cycle

Equation (3) with a relatively simple non-monotonic $g$ can demonstrate a different sort of a homoclinic bifurcation arising from the tangency of a stable and an unstable manifolds of a saddle (hyperbolic) periodic orbit, which can lead to chaos. This global bifurcation was theoretically discovered for ODEs in [57, 58], and described in a more accessible manner in Sec. 7.2.1 of [20]. In [49, 59] the existence of such chaos has been proved for very specific DDEs, for which a solution could be constructed analytically. Also, some general properties of function $g$ of (3) have been outlined, which should lead to similar behaviour. Namely, $g(z)$ should cross zero twice: from below to above at smaller $z$, and from above to below at larger $z$. Here we give an example of a smooth function g(z)$$=$$ah(z) for (3), with $h(z)$ expressed as

[TABLE]

and illustrated in Fig. 4(a), which for a$$>[math] posseses the properties given in [49, 59]. The landscape $V_{2}$, satisfying $h(z)=-\frac{\mathrm{d}V_{2}(z)}{\mathrm{d}z}$, is given in Fig. 4(b). In this Section, we study (4) with $h$ given by (5), which is equivalent to (3) with g(z)$$=$$ah(z). We will follow the behaviour of (4), (5) as control parameter $a$ is varied.

There are two fixed points here: x_{1}$$=$$-1.38629 (maximum of $V_{2}$, green circle in Fig. 4), and x_{2}$$=$$0.962424 (minimum of $V_{2}$, red circle in Fig. 4), which are both saddle for the range of $a$ illustrated in Figs. 4 and 5.

The saddle cycle, whose manifolds eventually become tangent to each other at the homoclinic bifurcation, is born via Andronov-Hopf bifurcation from the saddle point $x_{1}$ at a_{AH}$$=$$\frac{1}{h^{\prime}(z)}\frac{3\pi}{2}$$\approx$$2.8274 (as explained in Section III), and exists for a$$>$$a_{AH}. This saddle cycle (not shown in Fig. 4) is located close to its parent $x_{1}$ (green circle).

Like in the example of Section II.2, here the homoclinic orbit arises from the tangency of the manifolds of the saddle cycle as the parameter $a$ decreases to its critical value $a^{*}\gtrapprox 3.913$ from above. The manifolds involved in the formation of this orbit are codimension-one (i.e. infinite-dimensional) stable, and two-dimensional unstable, and cannot be visualised with the numerical tools available to date. However, chaos just born from the manifolds ceasing to be tangent and instead intersecting transversally is visualised with blue line in Fig. 4(c) at a$$=$$3.913.

Generally, in Fig. 4(c)–(f) phase portraits are shown for (4), (5) for several values of $a$ as $a$ decreases. Attractors originating from this homoclinic bifurcation are shown by blue line, and attractors originating from the fixed point $x_{2}$ at the landscape minimum (red cirlce) are shown by orange line.

The description below can be compared with the bifurcation diagram in Fig. 5. At a$$=$$3.913 (Fig. 4(c)) and a$$=$$3.85 (Fig. 4(d)) chaos born from the homoclinic bifurcation is the only attractor. At a$$\in$$[3.77,3.84] two attractors coexist. Figure 4(e) shows at a$$=$$3.8325 chaos born from the homoclinic orbit to the saddle cycle together with the coexistent stable cycle born through Andronov-Hopf bifurcation from $x_{2}$. In Fig. 4(f) at a$$=$$3.773 two coexisting stable cycles are given: the one born from homoclinic chaos in the inverse cascade of period-doubling bifurcations as $a$ decreases (blue line), and another born from $x_{2}$ as $a$ increases (orange line).

III Eigenvalues of the fixed points

From Section II it becomes apparent that in (2) with a smooth $V$ having minima and maxima, one should expect an interplay between periodic solutions born from the minima, and attractors born from homoclinic bifurcations associated with the maxima directly or indirectly. To predict the behaviour of the solution near a minimum, or the one resulting from a homoclinic bifurcation related to a maximum, one needs to know the eigenvalues of the relevant fixed points. In Section S-I of Supplementary Note we performed the linear stability analysis of the fixed points of (2) and showed that their eigenvalues $\lambda_{k}$ can be expressed in terms of the Lambert function $W(z)$ with z$$\in$$\mathbb{R} and W$$\in$$\mathbb{C} as

[TABLE]

where J$$=$$f^{\prime}(x^{*}) and $x^{*}$ is the fixed point. Function $W(z)$ has countably many branches, as illustrated in Fig. S2, and $W_{k}$ is its $k$-th branch. Therefore, \frac{\lambda_{k}}{J}$$=$$\frac{W_{k}(z)}{z}$$=$$\frac{W_{k}(J\tau)}{(J\tau)}. Figures 6 and 7 show eigenvalues of the fixed points at the landscape minimum and maximum, respectively (compare with Fig. S2).

Note, that at a landscape minimum J$$<[math], so with \tau$$\geq[math], z$$=$$J\tau$$\leq[math]. However, in Fig. 6 for convenience we show \frac{\lambda_{k}}{|J|}$$=$$\frac{W_{k}(-\tilde{z})}{\tilde{z}} as a function of \tilde{z}$$=$$|J|\tau$$\geq[math]. For |J|\tau$$\in$$\big{(}0,\frac{1}{e}\big{)} the leading eigenvalue $\lambda_{1}$ (green line) is real and negative, whereas all other eigenvalues are complex with large negative real parts (this is well visible in Fig. 6(b), (d)). Therefore, the fixed point is effectively a stable node, and the solution converges to it without oscillations. At |J|\tau$$\in$$\big{[}\frac{1}{e},\frac{\pi}{2}\big{)} there is a pair of complex-conjugate leading eigenvalues with negative real parts (red lines in Fig. 6), so in the centre manifold of Andronov-Hopf (AH) bifurcation that occurs at |J|\tau$$=$$\frac{\pi}{2}, the fixed point is a stable focus and the solution converges to it in an oscillatory manner. The boundary between these two subtly different types of behaviour is \tau$$=$$\frac{1}{e}, which is highlighted by vertical green dashed line in Fig. 6, and the respective values of real and imaginary parts of the eigenvalues are highlighted by green filled circles.

In Section S-I of Supplementary Note we derive the first condition for the first AH bifurcation of the fixed point $x^{*}$ at the minimum of $V$ (Eq. (S11))

[TABLE]

and verify the second condition (Eq. (S18)). Equation (7) is consistent with the predictions for the existence of a periodic solution in (3) for a special form of $g(z)$ as discussed in Section II.1. However, this result is more general and applies to a smooth $f$ of any shape. Thus, Eq. (7) allows one to determine the value of $\tau$ at which it is possible for the stable cycle to be born from the fixed point at the landscape minimum. The first AH bifurcation, AH1, is highlighted by the vertical red dashed line in Fig. 6, and the respective values of real and imaginary parts of the eigenvalues are highlighted by red filled circles. At \tau_{\textrm{AH2}}$$=$$\frac{5\pi}{2|J|} the second AH bifurcation occurs, as a result of which a saddle cycle could be born around the minimum. The second AH bifurcation, AH2, is highlighted with the vertical dashed cyan line in Fig. 6, and the respective values of real and imaginary parts of the eigenvalues are highlighted by yellow filled circles.

For a fixed point at the maximum of $V$, J$$>[math], and Fig. 7 shows \frac{\lambda_{k}}{J}$$=$$\frac{W_{k}(J\tau)}{J\tau} as a function of J\tau$$\geq[math]. One can see that, as also confirmed by the analysis in Section S-I of Supplementary Note, at any \tau$$\geq[math] there is one real positive eigenvalue $\lambda_{1}$ (orange line), so this fixed point is always unstable for non-negative $\tau$. However, it is important to appreciate that for J\tau$$\in$$\big{(}0,\frac{3\pi}{2}\big{)}, this point is a saddle-focus with a one-dimensional unstable manifold and an infinite-dimensional (codimension-one) stable manifold. As shown in (S19) of Supplementary Note, at J\tau$$=$$\frac{3\pi}{2}, or at

[TABLE]

the first AH bifurcation occurs, which is highlighted by a vertical red dashed line in Fig. 7, and the respective values of eigenvalues are marked by red filled circles. For visualisation purposes, for J\tau$$\in$$\big{(}0,\frac{3\pi}{2}\big{)} we can mentally replace the stable infinite-dimensional manifold of the saddle-focus with a two-dimensional centre manifold of its first AH bifurcation. In this approximation, the given fixed point would represent a saddle-focus with a one-dimensional unstable manifold and a two-dimensional stable manifold, as explained in more detail in Section IV. If there is a well of $V$ nearby, the stable and unstable manifolds can form a homoclinic loop at some $\tau$ from inside \big{(}0,\frac{3\pi}{2J}\big{)}. With this, Shilnikov’s theorem for ODEs [55, 56] verified for a special form of (3) [50, 52], predicts that for a safe loop with \sigma$$<[math] (see II.2 for the definition of $\sigma$ and Shilnikov’s theorem), which exists for J\tau$$\in$$[0,1.975), the homoclinic loop breaks down to form a stable periodic orbit. If the loop is dangerous with \sigma$$>[math], the resultant regime should be chaotic at least in ODEs [55, 56], although to the best of our knowledge this was not verified for DDEs.

Note, that at $\tau_{\textrm{AH1}}$ a saddle periodic orbit is born from the saddle point at the maximum, which is an intersection of a stable and an unstable manifolds. At larger $\tau$ the manifolds of this saddle orbit can form a loop, whose breakdown can give birth to chaos, as discussed in Section II.3.

IV Delay-induced behaviour in systems with two-well potentials

A two-well landscape function $V$ in (2) can be constructed by gluing together segments of single-well landscapes considered in Section II and smoothing out the joints. This observation leads us to suggest that the phenomena discussed in Section II should also occur in different parts of the phase space of (2) with a two-well $V$. However, a two-well function has a different quality as compared to a single-well one, so it is reasonable to also expect new phenomena not covered in Section II.

It follows from Section II that homoclinic bifurcations play a central role in DDEs (3) with non-monotonic $g(z)$ of even very simple shapes, and result in the disappearance of local attractors when $g$ becomes sufficiently steep. The nature of the homoclinic bifurcations depends on the fine local features of $g$, but their occurrence seems inevitable as the steepness parameter grows. For (2), an equivalent of the steepness parameter is $\tau$, so the homoclinic bifurcations are expected to occur as $\tau$ grows.

In this Section we verify our prediction, that homoclinic bifurcations should occur in systems (2) with two-well potentials $V$. We reveal that the occurrence of the particular forms of such bifurcations depends on the individual features of the potential wells and a hump, such as the depth, width and sharpness, as well as on the relationships between them.

We illustrate our findings using two subtly different examples of $V$ which lead to homoclinic bifurcations of different kinds. For this purpose, we construct the functions $V$ in such a way, that by adjusting parameters we could control their local features. The first example is considered in this Section, and the second one in Supplementary Note, Section S–II. However, our additional studies with the forms of $V$ specified in Supplementary Note, Section S–III, which are not reported here due to the lack of space, suggest the generality and reproducibility of the phenomena discovered for the landscapes with the same qualitative features.

In relation to optimisation, a particularly important phenomenon at large delays is the disappearance of localised attractors and the existence of the global chaos, which enables the phase trajectory to spontaneously visit the neighbourhoods of all local minima. This effect is similar to the action of random noise on (1) within simulated annealing [38] and superficially might seem quite simple. However, the mechanisms leading to global chaos are quite intricate and involve considerable rearrangements of manifolds belonging to saddle fixed points or to saddle cycles while the system undergoes a chain of homoclinic bifurcations. It is thanks to these global bifurcations that the barriers separating the vicinities of local minima disappear. Therefore, understanding reorganisation of manifolds is crucial for understanding why and how the delay could help in optimisation. In this Section we explain what happens to manifolds by applying the knowledge about homoclinic bifurcations available from the qualitative theory of ODEs, and verify its predictions with numerical simulations of (2).

IV.1 Model

Here we consider bifurcations in (2) with a double-well $V$ and the respective f$$=$$-V^{\prime} specified as follows

[TABLE]

The functions $V(x)$, $f(x)$ and J(x)$$=$$f^{\prime}(x) are shown in Fig. 8(a) by blue, red and green lines, respectively. The constant $(-0.88)$ was added in $V$ in order to make the graphs better distinguishable from each other. The upper part of Table I summarises the local features of $V$ in (IV.1), including the positions of two minima $x_{1,2}^{\mathrm{min}}$ and one maximum $x^{\mathrm{max}}$, their depths $V_{1,2}^{\textrm{min}}$ and $V^{\textrm{max}}$, Jacobians $J_{1,2}^{\textrm{min}}$ and $J^{\textrm{max}}$, and the points of the first AH bifurcations $\tau_{1,2}^{\textrm{AH1}}$ and $\tau^{\textrm{AH1}}$. In Fig. 8(a), red circles indicate the positions of $x_{1,2}^{\mathrm{min}}$, which are stable fixed points of (2) at \tau$$=[math], and green circle shows the fixed point at the landscape maximum $x^{\mathrm{max}}$, which is unstable at \tau$$=[math] and saddle at any \tau$$>[math].

In the context of optimisation illustrated in Section V, we remark that in this $V$ the minimum $x_{1}^{\mathrm{min}}$ is the lowest of the two, and the respective well is the broadest and has the flattest bottom.

IV.2 Overview of bifurcations

The single control parameter in (2), (IV.1) is $\tau$, and the bifurcation diagram is shown in Fig. 8(b). Namely, vertical red lines show the locations of fixed points $x_{1}^{\mathrm{min}}$ and $x_{2}^{\mathrm{min}}$ for the range of $\tau$ inside which they remain stable. At \tau$$=$$\tau_{1,2}^{\textrm{AH1}} where the lines stop (values are given in Table I) AH bifurcations occur in agreement with (7). Green vertical line shows $x^{\mathrm{max}}$ in the whole range of $\tau$ considered.

For oscillatory attractors, red/black dots indicate local maxima/minima of $x(t)$, which constitute the projections of the Poincaré sections defined as \dot{x}$$=[math], \ddot{x}$$<[math]/\ddot{x}$$>[math], respectively. Note, that this diagram is different from the classical bifurcation diagrams, such as the one in Fig. 5, which usually show only one kind of the attractor extrema. We show both kinds of the extrema because we need to illustrate how the sizes of the limit cycles born from AH bifurcations grow with $\tau$, and how the attractors approach the saddle point $x^{\mathrm{max}}$. This representation allows us to register homoclinic bifurcations from the abrupt changes in the locations and/or amplitudes of the numerically found attractors under very small changes in $\tau$.

Altogether there are three homoclinic bifurcations, as indicated in Fig. 8(b) and in Table I. Whereas the 1st homoclinic bifurcation can be easily detected from Fig. 8(b), the details of the 2nd and the 3rd ones are hard to see in this scale. To demonstrate these convincingly, an enlarged segment of the bifurcation diagram in the relevant area is provided in Fig. 9.

Below we will state our predictions about the details of delay-induced homoclinic bifurcations, which can be made on the basis of rigorous results overviewed in Sections II.3 and III. These will be verified and illustrated with phase portraits in Fig. 10, which complement the bifurcation diagram in Fig. 8(b).

The point $x^{\mathrm{max}}$ is of a saddle-focus type for \tau$$>[math]. Within the range of $\tau$ shown in Fig. 8(b), it has a single real positive eigenvalue $\lambda_{1}$ corresponding to a one-dimensional unstable manifold, and countably many complex-conjugate eigenvalues with negative real parts corresponding to a stable manifold of codimension one. In other words, $x^{\mathrm{max}}$ does not undergo AH bifurcation and does not give birth to a saddle cycle. Therefore, we should not be expecting the situation described in Section II.3, and all homoclinic bifurcations here must consist in the closure of the manifolds of $x^{\mathrm{max}}$ itself to form a one-dimensional loop, following the scenario of Section II.2.

Based on Section II.2, we can predict that the formation from the loop of a new attractor, which would be localised near $x_{1,2}^{\textrm{min}}$, would occur as $\tau$ passes the bifurcation point from above to below, i.e. decreases. If we consider the sequence of events as $\tau$ is increased, then the localised attractor existing at $\tau$ below the bifurcation point would approach $x^{\mathrm{max}}$, collide with it while the manifolds close at the bifurcation, and cease to exist for $\tau$ above the bifurcation.

The predicted scenario is confirmed for the 1st and 3rd homoclinic bifurcations illustrated with phase portraits in Fig. 10 (a)–(d) and (i)–(l), respectively. As a result of these bifurcations, the localised attractors disappear one after another. However, the 2nd homoclinic bifurcation does not follow from Section II.2 and is of a different type. Namely, as $\tau$ grows, this bifurcation produces a large limit cycle embracing all fixed points (see Fig. 10(e)–(h)).

Figure 9(a)–(b) shows that in a small range \tau$$\in$$[2.4307,2.4499], i.e. between the 2nd and the 3rd homoclinics, the system demonstrates bistability as the large cycle coexists with the smaller one around $x_{1}^{\textrm{min}}$, until the latter disappears via the 3rd homoclinic. There is also a hysteresis meaning that by slowly and continuously increasing and decreasing $\tau$, one observes different attractors. Namely, Fig. 9(a) shows all local maxima of an attractor as $\tau$ increases (blue filled circles) and as $\tau$ decreases (orange empty circles). The lack of the full coincidence between the blue and the orange circles is an evidence of hysteresis. Also, a hysteresis loop is shown by green arrows in (b) where only the largest maxima of the attractors are given for clarity.

To predict whether the attractors colliding with $x^{\mathrm{max}}$ at homoclinic bifurcations are limit cycles or chaotic attractors, we need to establish whether the respective loops are safe or dangerous, respectively, based on the sign of the saddle quantity $\sigma$, as explained in Sections II.2 and III. Namely, as mentioned in Section III, for DDEs the quantity $\sigma$ is negative and the loop is safe when \tau J^{\mathrm{max}}$$<$$1.975, implying the collision with a stable limit cycle. For all three homoclinic bifurcations, Table I gives the values of $\tau$ together with $\tau J^{\mathrm{max}}$, the latter being smaller than $1.975$. Therefore, we expect safe loops in all three homoclinic bifurcations.

Phase portraits in Fig. 10(a)–(l) comply with this prediction and show limit cycles being destroyed by ((a)–(d) and (i)–(l)), or born from ((e)–(h)), the homoclinic bifurcations as $\tau$ increases.

IV.3 Role of manifolds in homoclinic bifurcations

The key components of homoclinic bifurcations are invariant manifolds, and it is their reconfiguration which induces most drastic changes in the observed system behaviour associated with the death or the birth of attractors. Here, at \tau$$>[math] different basins of attraction are separated by the stable manifolds of the saddle-focus $x^{\mathrm{max}}$. Therefore, in the context of optimisation, these stable manifolds form the barriers between the vicinities of the local minima of $V$. Therefore, to understand if and how local attractors in (2), (IV.1) can disappear as $\tau$ grows, one needs to understand how these manifolds are reconfigured.

Whereas the attractors of nonlinear DDEs can be detected numerically with relative ease, revealing and plotting manifolds of such systems is a challenge. Some methods to visualise two-dimensional manifolds in nonlinear ODEs have been introduced in [60, 61], and unstable manifolds in DDEs in [62, 29]. However, the phase space of a DDE is infinite-dimensional, implying that the stable manifold of $x^{\mathrm{max}}$ has dimension infinity and is thus highly challenging to detect and to depict. Moreover, revealing a stable manifold in an ODE would require reversing time, which is generally not possible in DDEs [63]. With this, to the best of our knowledge, there are no techniques available to date to numerically obtain stable manifolds of DDEs.

While being unable to reveal by direct numerical visualisations how the manifolds reorganise in DDEs, we can still explain this at a qualitative level by extending the results from the qualitative theory of ODEs [21], which are supported by numerical visualisations of such reorganisations during homoclinic bifurcations of a saddle-node [64] and a saddle-focus [65] in a three-dimensional system of ODEs. Our explanations and predictions are illustrated with Figs. 11, 12 and 13, which can be compared with the phase portraits in Fig. 10. Although the phase portraits can supply only a highly limited information, their agreement with the predicted effects would provide a reasonably acceptable level of verification for these.

The manifolds of a saddle-focus near a homoclinic loop are quite intricate in shape. To make our explanation clearer, we start from considering some fictional two- and three-dimensional systems, which we assume to demonstrate similar homoclinic bifurcations occurring to a saddle-node, whose manifolds are simpler than those of a saddle-focus (see Figs. 11(a)–(j) and 12).

Next, in order to schematically illustrate homoclinic bifurcations in a DDE, we will proceed by analogy with a centre manifold reduction in ODEs [66]. Namely, we will assume that the dynamics on the infinite-dimensional (codimension-one) stable manifold of the saddle-focus $x^{\mathrm{max}}$ of (2), (IV.1) can be approximated by the dynamics on the two-dimensional manifold associated with the leading pair of complex-conjugate eigenvalues of this point. Thus, we will explain reorganisation of manifolds of a saddle-focus in (2), (IV.1) by sketching manifolds, fixed points, limit cycles and homoclinic loops in a three-dimensional phase space in Fig. 11(k)–(o) and Fig. 13.

IV.4 First homoclinic bifurcation

In Fig. 11, panels (a)–(e) show snapshots of manifolds, fixed points and limit cycles at consecutive values of a certain control parameter of a fictional two-dimensional system undergoing the simplest orientable homoclinic bifurcation, which is roughly similar to the 1st homoclinic bifurcation in (2), (IV.1). Here, the saddle-node (green circle) is an equivalent of $x^{\mathrm{max}}$. Dashed lines show manifolds affected by the homoclinic bifurcation, and solid lines show manifolds unaffected by it. Other notations are given in figure caption.

An equivalent sequence of events in a three-dimensional system is illustrated in Fig. 11(f)–(j). The only difference from Fig. 11(a)–(e) is that the stable manifold of the saddle-node is two-dimensional here. Depiction of all manifolds in a 3D space helps to make a transition from a homoclinic loop of a saddle-node to that of a saddle-focus, because the latter can exist only in spaces of dimension three and higher.

Next, Fig. 11(k)–(o) schematically illustrates manifolds of the saddle-focus and other objects in the phase space of (2), (IV.1): (k)–(m) before, (n) at, and (o) after the 1st homoclinic bifurcation. The notations are specified in the figure caption.

The saddle-focus (green filled circle in Fig. 11(k)–(o)) is an intersection of the two-dimensional stable manifold (cyan surface, a segment is shown) and the one-dimensional unstable manifold (blue line) associated with the only real and positive eigenvalue of this point. The stable manifold is the boundary between the attractor basins of two attractors (red filled circles or red lines).

Below we compare Fig. 11(k)–(o) with the bifurcation diagram in Fig. 8(b) and with the phase portraits in Fig. 10(a)–(d). As $\tau$ grows from zero, the following events precede the 1st homoclinic bifurcation. At \tau$$\in$$[0,0.7155] there are two coexisting stable fixed points (Fig. 11(k)) at $x_{1}^{\mathrm{min}}$/$x_{2}^{\mathrm{min}}$ below/above the stable manifold of $x^{\mathrm{max}}$ (compare with Fig. 8(b)).

At \tau$$\approx$$0.7155, $x_{2}^{\mathrm{min}}$ undergoes AH bifurcation and the stable cycle is born from it. The situation just above this bifurcation is illustrated in Fig. 11(l), where red line above the stable manifold shows the newly born cycle. At \tau$$\approx$$1.208, $x_{1}^{\mathrm{min}}$ undergoes AH bifurcation and the stable cycle is born from it, which for a slightly larger $\tau$ is given by red line below the cyan surface in Fig. 11(m) (compare with Fig. 10(a)). With this, the first limit cycle (upper red line in 11(m) and turquoise in 10(a)) has grown in size and stretched towards the saddle focus.

At \tau$$\approx$$1.53, the 1st homoclinic bifurcation occurs, namely, the limit cycle around $x_{2}^{\mathrm{min}}$ collides with $x^{\mathrm{max}}$, the unstable manifold of the latter “sticks" to its own stable manifold and forms a homoclinic loop (magenta lines in Fig. 11(n) and in Fig. 10(b)–(c)). This bifurcation is identical to the one illustrated in Fig. 3 as one goes from (f) to (d).

As $\tau$ increases slightly beyond the value of the 1st homoclinic $1.53$, the local attractor on the right side of $x^{\mathrm{max}}$, i.e. near $x_{2}^{\mathrm{min}}$, ceases to exist. The system has a single attractor left, which is the limit cycle localised near $x_{1}^{\mathrm{min}}$ (compare red lines in Fig. 10(d) and Fig. 11(o)). If at $\tau$ just above $1.53$ the initial conditions are set on or near the just disappeared local attractor around $x_{2}^{\mathrm{min}}$, the phase trajectory swiftly leaves this region and converges to the limit cycle around $x_{1}^{\mathrm{min}}$ in what might feel like a “jump" as indicated in Fig. 8(b).

Unlike in the homoclinic bifurcation of a saddle-node illustrated in Fig. 11(a)–(j), the shape of the stable manifold of the saddle-focus is quite intricate. Namely, close to the homoclinic bifurcation, part of the stable manifold of $x^{\mathrm{max}}$ returns to the vicinity of $x^{\mathrm{max}}$ from above and takes the shape of an open helicoid with a finite number of turns (see upper part of Fig. 11(k)-(m)), with turns coming closer to each other as they come closer to the saddle-focus, as numerically demonstrated for a 3-dimensional dynamical system in [65].

At the instant of homoclinic bifurcation (n), the helicoid becomes closed and develops an infinite number of turns, such that the distance between the consecutive turns becomes smaller as the turns come closer to $x^{\mathrm{max}}$. At the same instant, the unstable manifold “sticks" to this closed helicoid and forms a loop (n). After the homoclinic bifurcation, the helicoid opens and possesses a finite number of turns again (o). As a result of homoclinic bifurcation, the one-dimensional unstable manifold (blue line) “permeates" through the stable one (cyan surface), and the latter no longer separates basins of attraction. A more detailed discussion of the attractor basins is given in Section IV.5.

Figure 11(k)-(o) schematically shows only a small portion of the stable manifold near the homoclinic loop. However, in reality this manifold has an even more complex shape. E.g. in (o) a part of this manifold should unwind from the vicinity of the unstable fixed point (empty circle in the upper part), but we cannot show this without overloading the figure. Also, this manifold extends well beyond the boundaries shown and has a similarly complex structure in the lower parts of panels in Fig. 11(k)-(o), which we again do not show for the sake of clarity. Thus, the complexity of the shape of the stable manifold of the saddle-focus prevents us from illustrating this bifurcation in full.

With this, reorganisation of the manifolds away from the homoclinic loop during this bifurcation would be qualitatively the same if the saddle-focus is replaced by a saddle node, whose manifolds have a simpler shape and are easier to depict (see Fig. 11(f)–(j)).

IV.5 Basins of attraction

In the context of optimisation, our main interest is in how homoclinic bifurcations amend the basins of attraction and lead to their merging. For a fictional two-dimensional system, Fig. 11(a)–(e) shows the basins of two coexisting attractors by different shades. Depiction of the basins is easy on the plane, but difficult in higher-dimensional spaces, so this figure can be used for a visual reference. Here, yellow shade shows the basin of an attractor at or near an equivalent of $x_{2}^{\mathrm{min}}$ (upper right part of the panel), whereas the basin of attractor at or near $x_{1}^{\mathrm{min}}$ (lower left part of the panel) is not shaded.

By analogy with similar situations in ODEs, one can hypothesise that the stable manifold should make several turns around the three fixed points, thus making both basins of attraction stripy. Although we cannot verify this hypothesis by building stable manifolds of the DDE, we can numerically find the basins themselves and reveal the stable manifolds as their boundaries, as was done in [67] for a four-dimensional system.

The basin of attraction in (2) is a set of all initial functions $\varphi(t)$, $t\in[-\tau,0]$, such that all solutions starting from these functions converge to the given attractor. It is a subset of an infinite-dimensional phase space and difficult to visualise in full. However, we utilise an approach of [68, 69] and consider only a subset of all possible initial conditions, namely, functions of a certain class parametrised by a finite number of parameters. In our case we choose the simplest class of constant initial functions, \varphi(t)$$=$$x$$=$$\mathrm{const}, $t\in[-\tau,0]$. Such functions can be directly compared with the positions of fixed points in (2).

Importantly, numerical simulation of (2), (IV.1), as well as of (2) with several similar potentials, including those specified in Sections S–II and S-III of the Supplementary Note, revealed that at positive $\tau$ from some initial conditions such systems systematically go to infinity. Thus, infinity represents an additional attractor with its own basin and boundaries, which could be revealed numerically.

The basins of attraction of (2), (IV.1) are given in Fig. 14 for a range of $\tau$ in the vicinity of the 1st homoclinic bifurcation. The striking feature of this plot is the large basin of attraction of infinity (blue shade), which expands with growing $\tau$. Below the bifurcation point \tau$$=$$1.53, in addition to infinity, there are two attractors on different sides of $x^{\mathrm{max}}$, and on the plot we do not distinguish between fixed points and limit cycles. Due to the 1st homoclinic bifurcation at \tau$$\approx$$1.53, the basin of attractor at or near $x_{2}^{\mathrm{min}}$ (red shade) collapses, and at 1.53$$<$$\tau$$<$$2.43 besides the infinity there is only one attractor at or near $x_{1}^{\mathrm{min}}$ (green shade).

Figure 14 reveals that the basins of attraction at least before the 1st homoclinic bifurcation are stripy, and therefore confirms the hypothesis that the stable manifold of $x^{\mathrm{max}}$ makes several turns around the three fixed points (compare with Fig. 11(b)–(c)).

IV.6 Second homoclinic bifurcation

As mentioned in Section IV.2, the 2nd homoclinic bifurcation cannot be predicted based on the theorems overviewed in Section II.2. The numerical simulations suggest that it consists in forming a large homoclinic loop embracing all fixed points (magenta line in Fig. 10(f)), and results in the birth of a large limit cycle (turquoise line in Fig. 10(h)).

The formation of the large homoclinic loop involves a different pair of manifolds of the saddle-focus $x^{\mathrm{max}}$, as compared to those forming the 1st small loop. As before, given the difficulty of making clear and easily interpretable sketches of intricately shaped manifolds of a saddle-focus, we initially illustrate the respective homoclinic loop in a two- (Fig. 12(a)–(c)) and a three-dimensional (Fig. 12(f)–(h)) fictional systems with saddle-nodes. Notations are the same as in Fig. 11, and in addition brown lines in (b)–(e) show typical phase trajectories.

Figure 12(a)–(c) explains how the basin of a newly born attractor (large cycle) is formed thanks to the rearrangement of manifolds. Namely, in (a) yellow shade shows the basin of the only attractor available just before the bifurcation, i.e. of the limit cycle around an equivalent of $x_{1}^{\mathrm{min}}$ (red line, compare with Fig. 10(d)–(e)). At the 2nd homoclinic bifurcation in (b), this basin becomes bounded. After the bifurcation in (c), the manifolds are rearranged to form the boundary of the new basin (non-shaded) of the large limit cycle born from the homoclinic loop (large red closed curve).

Figure 12(f)–(h) demonstrates the same sequence of events, only in a three-dimensional system with a two-dimensional stable manifold of the saddle-node. This illustration is an intermediate stage before depicting manifolds of a saddle-focus in the three-dimensional space while they undergo a similar reorganisation.

The respective manifolds of a saddle-focus of (2), (IV.1), in the centre manifold reduction and in the vicinity of the large homoclinic loop, are sketched in Fig. 13. Panels (a)–(b) illustrate the situation just before the 2nd (large) homoclinic loop is formed (compare with Fig. 12(f)), and (c)–(d) illustrate the instant of this large loop formation (compare with Fig. 12(g)).

The picture of the manifolds here is even more complex than near the small homoclinic loop sketched in Fig. 11(n) because of the spiralling of the stable manifold as it approaches the saddle focus from both sides. However, away from the saddle focus, the stable manifold is qualitatively similar to the one of a saddle-node.

Before the homoclinic loop, the stable manifold, whose different segments are shown by cyan surface in panels (a) and (b) of Fig. 13, wraps itself around all the three fixed points more than once, thus making more than one layer similarly to the manifold shown in Fig. 12(a) by dashed cyan line.

As a result, the one-dimensional unstable manifold of the saddle-focus (blue line in Fig. 13(a)–(b)) goes between the two layers of the stable manifold upwards from $x^{\mathrm{max}}$ and then towards the local attractor near $x_{1}^{\mathrm{min}}$ (red line in the lower left parts of the panels), just like in Fig. 12(a), (f). At the instant of the homoclinic bifurcation at \tau$$\approx$$2.4307, a segment of the stable manifold collides with the saddle-focus and captures the unstable manifold, which now forms a large safe homoclinic loop (magenta line in Fig. 13(d), compare with Figs. 10(f)–(g) and 12(b), (g)). Note, that at the instant of bifurcation, the helicoid-shaped portion of the stable manifold containing the homoclinic loop (upper part of Fig. 13(d)) makes an infinite number of turns as it approaches the saddle focus.

As $\tau$ grows, the loop disappears and gives birth to a large limit cycle embracing all three fixed points, which are similar to large red closed curves in Fig. 12(c), (h). The respective large cycle of (2), (IV.1) is given by turquoise line in Fig. 10(h).

IV.7 Third homoclinic bifurcation, chaos and infinity

The 3rd homoclinic bufurcation takes place at \tau$$=$$2.4499 (see Fig. 8(b)) and eliminates the small limit cycle to the left of $x^{\mathrm{max}}$, i.e. around $x_{1}^{\mathrm{min}}$, as illustrated with phase portraits in Fig. 10(i)–(l) and with sketches of a similar bifurcation in systems with a saddle-node in Figs. 12(c)–(e) and (h)–(j). This bifurcation is qualitatively the same as the 1st homoclinic. At \tau$$>$$2.4499 there are no more local attractors in the system.

At \tau$$\in$$(2.4499,3,33], the system has a single large attractor enclosing all fixed points. This attractor is initially a limit cycle of period one (Fig. 10(l)), but with increasing $\tau$ it undergoes a cascade of period-doubling bifurcations (Fig. 10(m)–(n)) and becomes chaotic ((Fig. 10(o)). Note, that the phase trajectory on the chaotic attractor visits the close vicinities of all three fixed points, which is roughly similar to the effect caused by adding large random noise to (1) in simulated annealing.

As $\tau$ exceeds $3.33$, the chaotic attractor disappears and the trajectory goes to infinity (Fig. 10(p)). We cannot specify the exact reason for this, and can only hypothesise that the manifold bounding the basin of attraction of chaos appears involved in some global bifurcation. Since we cannot visualise this manifold, we cannot verify our hypothesis. However, our studies of a considerable number of systems of the form (2) with various multi-well landscapes $V$ suggest the universality of this phenomenon.

Namely, it seems that in such systems, as $\tau$ becomes sufficiently large, a chaotic attractor enclosing all extrema of $V$ is initially born (although it might turn into a periodic one at even larger $\tau$ e.g. when periodic windows in chaos appear [70]). However, when $\tau$ is increased further, the large attractor is inevitably destroyed and the trajectory escapes to infinity from all initial conditions.

It is important to appreciate that the infinity is present as some kind of an attractor at all values of \tau$$>[math]. To make more specific predictions about encountering attractors at infinity, a rigorous investigation of the necessary and sufficient conditions for the existence of unbounded solutions in nonlinear systems of type (2) would be required.

IV.8 Different forms of homoclinics

We also studied the phenomena induced by the increase of $\tau$ in a slightly different system of the form (2) with a double-well landscape $V(x)$ specified by Eq. (S21) in Supplementary Note, in which at the instants of homoclinic bifurcations the saddle focus was below AH bifurcation, but had a positive saddle quantity. According to Shilnikov’s theorem for ODEs [55, 56], in such cases the homoclinic loops are expected to be dangerous, and their breakdown should lead to the formation of chaotic attractors. However, to the best of our knowledge, this result has not been verified for DDEs. Our numerical studies did not reveal any obvious differences between the sequence of bifurcations in the systems of the form (2) with double-well landscapes $V$ with dangerous or safe homoclinic loops, and the observed phenomena looked very similar when studied with the same numerical accuracy. Thus, whether the dangerous homoclinic loop gave birth to chaos or not, it did not affect the order of bifurcations, and the same key phenomena were observed.

In addition, Section S–II of the Supplementary Note presents a similar, albeit briefer, analysis of bifurcations in (2) with a two-well potential $V$, for which at the instant of the 1st homoclinic bifurcation the fixed point $x^{\mathrm{max}}$ is past AH bifurcation. There, the homoclinic orbit is formed by the manifolds of the saddle cycle born from $x^{\mathrm{max}}$, which become tangent to each other, as described in Section II.3. Nevertheless, the general sequence of bifurcations there is very similar to the one in (2), (IV.1).

V Relevance to optimisation

In the context of optimisation, the delay-induced bifurcations in (2) inform us about the mechanisms behind the removal of the barriers between the local minima of the two-well potential function, which could be potentially applicable to multi-well landscapes as well. Indeed, in the DDE (2) the barriers between the local minima of $V$ are formed by the stable codimension-one manifolds of the saddle fixed point at the maximum of $V$, or of the saddle cycle born from this point. These manifolds are the boundaries of the basins of attractors localised near the minima of $V$. Delay-induced homoclinic bifurcations do not make these fixed points of manifolds disappear, but rearrange manifolds in such a way that they cease to separate different basins and hence to be the barriers.

The absence of barriers between the local minima and the birth of the large chaotic attractor enables the phase trajectory to visit the vicinities of all minima. Thus, the delay acts by a rough analogy with the random noise in a famous optimisation technique simulated annealing [38].

One possible optimisation procedure using the delay $\tau$ as control parameter is illustrated in Fig. 8(c). To obtain this figure, Eqs. (2), (IV.1) were solved as $\tau$ was decreased from some positive value to zero. Namely, at t$$=[math], the system was launched from initial conditions belonging to the large chaotic attractor at \tau$$=$$3.3. Then $\tau$ was monotonously decreased in a step-wise manner from $3.3$ to zero with small steps of the size $0.01$. The duration of each step was $50$ time units, during which the phase trajectory approached the attractor closest to the last state from the previous step. The decrease of $\tau$ was almost adiabatic.

With this procedure, at \tau$$=[math] the system ended up at the minimum $x_{1}^{\mathrm{min}}$, at which the modulus of the Jacobian, $|J|$, is smaller than at $x_{2}^{\mathrm{min}}$ (see Fig. 8(a)), and which therefore undergoes AH bifurcation at a larger value of $\tau$ (see Fig. 8(b)). This minimum $x_{1}^{\mathrm{min}}$ is not only flatter than $x_{2}^{\mathrm{min}}$, but also broader, due to which the local attractor around $x_{1}^{\mathrm{min}}$ survives for larger values of $\tau$ than the one around $x_{2}^{\mathrm{min}}$.

Under the adiabatic decrease of $\tau$ within approximately $[1.53,3.3]$, i.e. above the 1st homoclinic, at every instantaneous value of $\tau$ the solution is close to the only attractor corresponding to this value of $\tau$. The latter is either the one embracing two minima at larger $\tau$, or the one localised near, or at, $x_{1}^{\mathrm{min}}$ at smaller $\tau$. When $\tau$ decreases below the 1st homoclinic at \tau$$\approx$$1.53, the local attractor around $x_{2}^{\mathrm{min}}$ appears. However, the solution is not affected by this event and remains in the vicinity of $x_{1}^{\mathrm{min}}$ until $\tau$ reaches zero. Since in this example $x_{1}^{\mathrm{min}}$ happens to be the lowest minimum, adiabatic decrease of $\tau$ has demonstrated optimisation as required.

In Supplementary Note, Section S–II B, optimisation with adiabatic decrease of $\tau$ is demonstrated for a subtly different example of a system with a two-well potential, in which all homoclinic bifurcations take a different form as compared to the ones in (2), (IV.1). To summarise, adiabatic decrease of $\tau$ delivers the flattest and the broadest minimum of the two.

It is clear that in a general two-well potential such a minimum would not necessarily be the lowest one, and adiabatic decrease of $\tau$ would not result in optimisation. To overcome this issue and to reduce the dependence on the shape of the potential, we can use the fact that larger delay induces chaotic, i.e. random-looking, behaviour forcing the phase trajectory to visit the vicinities of all minima. Then at some sufficiently large value of $\tau$, at some randomly chosen time $t$ the system can be found in one of these vicinities. If at this moment $\tau$ starts to decrease relatively fast, i.e. in a non-adiabatic manner, the phase trajectory would fail to approach any attractor existing at any fixed value of $\tau$. This way, at \tau$$=[math] the system could end up in the same well it was before $\tau$ started to decrease, to whose bottom it would then converge following standard gradient descent (1). Several repetitions of the same experiment with $\tau$ being increased from, and then decreased to, zero would eventually reveal all available minima, whose depths could be compared at the end, and the lowest one would be identified.

To illustrate this effect, consider (2) with a double-well $V$ and the respective f$$=$$-V^{\prime} specified as follows

[TABLE]

The functions $V(x)$, $f(x)$ and J(x)$$=$$f^{\prime}(x) are shown in Fig. 15(a) by blue, red and green lines, respectively. Red/green circles indicate the positions of the fixed points at the potential minima/maximum. Here, the flattest and broadest minimum is $x_{1}^{\textrm{min}}$, but the lowest minimum is $x_{2}^{\textrm{min}}$.

First, we launch this system at t$$=[math] from the initial conditions x(t)$$=$$2 for $t\in[-2.1,0]$, and \tau$$=$$2.1, and start to find the numerical solution $x(t)$ while decreasing $\tau$ at a relatively low rate. Namely, after $0.1$ time units, $\tau$ is abruptly decreased to $2.0$ and kept at this value for another $0.1$ time units while the system is being solved. The process continues while $\tau$ is being decreased in a step-wise manner to zero in steps of $0.1$. The resultant solution $x(t)$ is given in Fig. 15(b) by black symbols and demonstrates that with this relatively slow decrease of $\tau$ the system ends up near the flattest and broadest minimum $x_{1}^{\textrm{min}}$.

Next, we repeat the process of decreasing $\tau$ step-wise, but make the duration of every step ten times smaller than above, i.e. $0.01$. The resultant solution is shown in Fig. 15(b) by red symbols, and demonstrates that the system ends up in the well of the lowest minimum $x_{2}^{\textrm{min}}$. In Fig. 15(c) the same solutions are shown as in (b), but now against the current value of $\tau$ for the convenience of comparison. Thus, the fast decrease of $\tau$ delivers the lowest minimum, which is not the broadest or the flattest.

VI Discussion and Conclusion

The fact of the occurrence of homoclinic bifurcations induced by the increase of delay in (2) with a multi-well potential $V$ could be predicted based on the theorems overviewed in Section II. However, the specific forms of these bifurcations, their dependence on the features of $V$, and their ordering do not follow from these theorems.

For a general nonlinear DDE with an arbitrary nonlinearity and an arbitrary dependence on delay, it is usually impossible to predict the delay-induced changes in the behaviour before actually observing this behaviour using numerical analysis. However, we hypothesised and verified that for a special class of scalar nonlinear models (2), in which the right-hand side depends only on the delayed variable and represents the negative of the gradient of a two-well potential, it is possible to make some qualitative predictions of the phenomena induced by the increase of the delay.

Specifically, we tested and confirmed our initial hypothesis, that the increase of the delay $\tau$ should lead to a chain of homoclinic bifurcations leading to the disappearance of attractors localised around the minima of $V$, and to the eventual birth of a large attractor embracing both local minima, which forces the system to visit the vicinity of every minimum, including the global minimum, as time goes by.

The latter resembles the effect from a random term within simulated annealing [38], but is achieved in a purely deterministic manner. It also has some similarity with quantum annealing assumed in quantum computers [71], which is performed thanks to the ability of particles to tunnel the potential barriers between the local minima of energy function. In our case, thanks to the delay, the barrier between the two minima of the cost function disappears, thus enabling the system to freely wander between both minima. Indeed, the barrier between the minima is ultimately the manifold of the saddle fixed point at the maximum, which separates the basins of attractors localised near the minima. Although after the homoclinic bifurcation neither this manifold, nor the maximum disappear, the manifold stops being the barrier, and hence the barrier ceases to exist.

Another important observation is that, for the values of the delay exceeding some threshold, no attractors at finite locations survive in the system, and the phase trajectory tends to infinity from any initial conditions. This effect has been observed numerically for all examples considered, however, it would be useful to verify its validity analytically in the future work. In this context, it would be interesting to understand whether the death of the last attractor at a finite location is linked to reorganisation of manifolds or is caused by other reasons. In the absence of numerical methods for visualisations of the relevant manifolds in DDEs, this matter cannot be resolved by their direct computation at this stage. With this, we are not aware of any relevant theoretical results which could readily suggest a plausible explanation behind this phenomenon. Hopefully, this could be clarified with further development of numerical approaches and/or theory of bifurcations in DDEs.

Our results show that the realisation of the particular forms of homoclinic bifurcations depends not only on the parameters of the potential wells and a hump, such as the depth, width and sharpness, but also on the relationships between them. Namely, the local attractor around one of the minima may collide either with the maximum itself, or with the saddle cycle born from this maximum. Prior to the disappearance of the local attractor through a homoclinic bifurcation, bistability may either occur, or fail to occur.

However, the exact nature of homoclinic bifurcations does not seem to change the general sequence of events as $\tau$ grows, which has been confirmed also with a number of additional examples, including those provided in Supplementary Note. Thus, we are satisfied with our ability to qualitatively predict the events induced by the increase of time delay in a highly nonlinear system of a special form.

Our results are a pre-requisite to understanding and prediction of the phenomena in general nonlinear delay-differential equations, in which the nonlinear function in the right-hand side depends solely on the delayed variable. They will inform building delay systems with prescribed controllable properties in applications. For example, if one wishes to obtain a system demonstrating a desirable sequence of bifurcations with the increase of delay, one could use the knowledge obtained in order to design the landscape $V$ of (2) with the necessary properties. Namely, one could fine-tune the relative depths, widths and shapes of the individual wells and humps in $V$ in order to ensure a certain order and forms of homoclinic bifurcations.

A considerable motivation for this study has been the idea of [41] to use delay-induced bifurcations for optimisation. Here we demonstrated how optimisation could be achieved for some two-well cost functions if one uses the delay as the only control parameter. Further studies will be needed to explore the possibility to extend this approach to cost functions depending on many variables.

VII Supplementary material

Supplementary Note contains an overview of the analysis of local dynamics of DDEs around the fixed points (Section S-I), considers an example of a delay system (2) with a double-well potential slightly different from that specified by (IV.1), where a similar, but subtly different, chain of homoclinic bifurcations is observed (Section S-II), and gives three additional examples of double-well potential functions resulting in a similar chain of bifurcations as in examples considered in paper (Section S-III).

VIII Authors’ contributions

NBJ proposed and designed the research, did all numerical simulations described in the paper except for Fig. 14, interpreted the results, and wrote the paper. CJM performed analytical treatment of a delay differential equation in Section III and in Supplementary Note, did preliminary numerical simulations of a system not illustrated in the paper, which delivered rough initial results, calculated basins of attraction for Fig. 14, and edited the paper. Both authors contributed to the discussions.

IX Data availability

Data sharing is not applicable to this article as no new data were created or analysed in this study.

X Acknowledgements

The authors are grateful to Dmitri Tseluiko for helpful discussions about Shilnikov’s theorem, and for help in preparing Figures 11(k)-(o); and to Alexander Balanov for helpful discussions and feedback on the draft of this manuscript. C.J.M. was supported by EPSRC (UK) during his PhD studies in Loughborough University, grant EP/P504236/1.