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

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.

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 multiwell "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 a) E-mail: N.B.Janson@lboro.ac.uk 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-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 2, 16,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 τ≥0, i.e. oḟ x= f (x, x τ ) with x, f ∈ R and x τ =x(t − τ) 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ẋ= f (x τ )−g(x), often with g(x)=λ x (λ ∈R) 15,17,22-25 , or to an even simpler forṁ x= f (x τ ), 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 forṁ where x(t)∈R is the state, t is the time,ẋ= dx dt , and V (x): R→R is the potential energy function at least twice continuously differentiable 34 . Also, ∇ is the gradient operator, which for a scalar x is equivalent to ∂ ∂ 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)∈R N , V (x): R N →R with N≥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 τ≥0 and thus obtain a DDE to study, where x, f , V , z ∈ R, x τ =x(t − τ). 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 τ 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)→∞ as |z|→∞ and V , z ∈ 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 delayinduced 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 -2 FIG. 1. Illustration of the relationship between equations (2) with arbitrary τ, and (3) with τ=1 and g(z)=τ f (z). For (2) with τ=5, 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 τ=1 and f is amendable 26,[42][43][44][45][46][47][48][49][50] , i.e. the equation reads However, introducing t=sτ and x(t)=y(s) reduces (2) to dy ds =τ f (y(s − 1)), i.e. to (3) with g(z)=τ 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 τ 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 ). 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 τ, 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.
A. 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=0, crossing zero at z=0, and in addition satisfying zg(z)<0 for all z =0. 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=0.
In addition, g(z) should either be bounded from below for all z (Fig. 2(b)), or satisfy |g(z)|≤A if |z|≤A, where A is some  positive constant (Fig. 2(d)). If |g (0)|> π 2 , then (3) has a nonzero 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=0. 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.
B. Existence of a homoclinic orbit to a fixed point Equation (3) with a non-monotonic g was studied in 52,53 but rewritten as Here, h(z) is required to cross zero at z=0 from below to above, and at z=−|b|, b∈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 =0 (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 =0 (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≈1.396 at which |ah (−2.5)| = π 2 . The point x 1 is unstable for any a>0. 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 * ≈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=0 starting from some time instant in the past and approaches x=0 as t→∞ 54 .
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>0, the fixed point x 1 =0 of (4) is a saddle focus with a single real positive eigenvalue λ 1 , whereas all of its other eigenvalues are complexconjugate 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. σ =λ 1 + Re(λ 2,3 )<0, where λ 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 σ of the fixed point x 1 =0. If σ >0, then in ODEs, according to Shilnikov, the breakdown of the homoclinic loop should produce chaos rather than periodic behaviour.
C. 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 man-  ner 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 and illustrated in Fig. 4(a), which for a>0 posseses the properties given in 49,59 . The landscape V 2 , satisfying h(z) = − dV 2 (z) dz , 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 FIG. 5. A segment of bifurcation diagram of (4) with h given by (5), which for every value of a shows the local minima of the solution x(t) (after transients are removed). Blue dots correspond to attractors originating from the homoclinic orbit to the saddle cycle located around the saddle fixed point x 1 at the maximum of V 2 (green line). Orange dots show attractor originating from the fixed point x 2 at the minimum of V 2 , which is outside the range of x of this figure. One can clearly see bistability for a∈[3.77, 3.85], which means that two attractors of different origins coexist. via Andronov-Hopf bifurcation from the saddle point x 1 at a AH = 1 h (z) 3π 2 ≈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 B, 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 * 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∈[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 perioddoubling 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 result-ing 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 λ k can be expressed in terms of the Lambert function W (z) with z∈R and W ∈C as where J= f (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, 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<0, so with τ≥0, z=Jτ≤0. However, in Fig. 6 for convenience we show z as a function ofz=|J|τ≥0. For |J|τ∈ 0, 1 e the leading eigenvalue λ 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|τ∈ 1 e , π 2 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|τ= π 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 τ= 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)) 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 A. 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 τ 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 τ AH2 = 5π 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 Jτ as a function of Jτ≥0. One can see that, as also confirmed by the analysis in Section S-I of Supple-mentary Note, at any τ≥0 there is one real positive eigenvalue λ 1 (orange line), so this fixed point is always unstable for non-negative τ. However, it is important to appreciate that for Jτ∈ 0, 3π 2 , this point is a saddle-focus with a onedimensional unstable manifold and an infinite-dimensional (codimension-one) stable manifold. As shown in (S19) of Supplementary Note, at Jτ= 3π 2 , or at 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τ∈ 0, 3π 2 we can mentally replace the stable infinite-dimensional manifold of the saddle-focus with a twodimensional 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 τ from inside 0, 3π 2J . 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 σ <0 (see II B for the definition of σ and Shilnikov's theorem), which exists for Jτ∈[0, 1.975), the homoclinic loop breaks down to form a stable periodic orbit. If the loop is dangerous with σ >0, 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 τ 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 τ the manifolds of this saddle orbit can form a loop, whose breakdown can give birth to chaos, as discussed in Section II C.

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 τ, so the homoclinic bifurcations are expected to occur as τ 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).

A. Model
Here we consider bifurcations in (2) with a double-well V and the respective f =−V specified as follows The functions V (x), f (x) and J(x)= f (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 (9), including the positions of two minima x min 1,2 and one maximum x max , their depths V min 1,2 and V max , Jacobians J min 1,2 and J max , and the points of the first AH bifurcations τ AH1 1,2 and τ AH1 . In Fig. 8(a), red circles indicate the positions of x min 1,2 , which are stable fixed points of (2) at τ=0, and green circle shows the fixed point at the land-  (9). Red/green circles show positions of fixed points at the minima/maximum of V . (b) Bifurcation diagram of (2), (9). Local minima/maxima of attractors are shown by black/red dots. Fixed points at the minima of V are shown by red vertical lines for τ at which they are stable, and the saddle-focus at the maximum of V is shown by green vertical line. As τ varies, in this system only safe homoclinic loops are formed by the manifolds of the saddle-focus fixed point x max . (c) Demonstration of optimisation. Local maxima (red dots) and minima (black dots) of solutions to (2), (9) are shown, which are obtained as τ slowly decreases from 3.3 to zero. The solution spontaneously settles down at the lowest, flattest and broadest minimum x min 1 . Compare with (a) and (b).
scape maximum x max , which is unstable at τ=0 and saddle at any τ>0.
In the context of optimisation illustrated in Section V, we remark that in this V the minimum x min 1 is the lowest of the two, and the respective well is the broadest and has the flattest bottom.  Fig. 8(b) illustrating 2nd and 3rd homoclinic bifurcations, bistability and hysteresis in (2), (9). Circles show all maxima of the attractors (limit cycles): as τ increases (blue filled), and as τ decreases (orange empty). Vertical green line shows x max . (b) The largest maxima of the attractors shown to demonstrate hysteresis more clearly: as τ increases (blue filled), and as τ decreases (orange empty). Hysteresis loop is marked by green arrows. The relevant small and large cycles are shown in Fig. 10(h) by red and turquoise lines, respectively.

B. Overview of bifurcations
The single control parameter in (2), (9) is τ, and the bifurcation diagram is shown in Fig. 8(b). Namely, vertical red lines show the locations of fixed points x min 1 and x min 2 for the range of τ inside which they remain stable. At τ=τ AH1 1,2 where the lines stop (values are given in Table I) AH bifurcations occur in agreement with (7). Green vertical line shows x max in the whole range of τ considered.
For oscillatory attractors, red/black dots indicate local maxima/minima of x(t), which constitute the projections of the Poincaré sections defined asẋ=0,ẍ<0/ẍ>0, 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 τ, and how the attractors approach the saddle point x 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 τ.
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 C 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 max is of a saddle-focus type for τ>0. Within the range of τ shown in Fig. 8(b), it has a single real positive eigenvalue λ 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 max does not un- Homoclinic bifurcation points in (2) (9), and the points of homoclinic bifurcations in (2), (9) as τ grows. The Table provides the locations of the minima x min 1,2 and the maximum x max , their depths V min 1,2 and V max , Jacobians J min 1,2 and J max , and the points of the first AH bifurcations τ AH1 1,2 and τ AH1 . For x max , τ σ is the value of τ 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 , which are safe since all values of τJ max are below 1.975. dergo AH bifurcation and does not give birth to a saddle cycle. Therefore, we should not be expecting the situation described in Section II C, and all homoclinic bifurcations here must consist in the closure of the manifolds of x max itself to form a one-dimensional loop, following the scenario of Section II B.
Based on Section II B, we can predict that the formation from the loop of a new attractor, which would be localised near x min 1,2 , would occur as τ passes the bifurcation point from above to below, i.e. decreases. If we consider the sequence of events as τ is increased, then the localised attractor existing at τ below the bifurcation point would approach x max , collide with it while the manifolds close at the bifurcation, and cease to exist for τ 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 B and is of a different type. Namely, as τ 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 τ∈[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 min 1 , until the latter disappears via the 3rd homoclinic. There is also a hysteresis meaning that by slowly and continuously increasing and decreasing τ, one observes different attractors. Namely, Fig. 9(a) shows all local maxima of an attractor as τ increases (blue filled circles) and as τ 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 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 σ , as explained in Sections II B and III. Namely, as mentioned in Section III, for DDEs the quantity σ is negative and the loop is safe when τJ max <1.975, implying the collision with a stable limit cycle. For all three homoclinic bifurcations, Table I gives the values of τ together with τJ max , the latter being smaller than 1.975. Therefore, we expect safe loops in all three homoclinic bifurcations.

C. 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 τ>0 different basins of attraction are separated by the stable manifolds of the saddle-focus x 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), (9) can disappear as τ 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 29,62 . However, the phase space of a DDE is infinitedimensional, implying that the stable manifold of x 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 threedimensional 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 max of (2), (9) 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), (9) by sketching manifolds, fixed points, limit cycles and homoclinic loops in a three-dimensional phase space in Fig. 11(k)-(o) and Fig. 13.

D. 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), (9). Here, the saddle-node (green circle) is an equivalent of x 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 saddlenode 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), 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 τ grows from zero, the following events precede the 1st homoclinic bifurcation. At τ∈[0, 0.7155] there are two coexisting stable fixed points (Fig. 11(k)) at x min 1 /x min 2 below/above the stable manifold of x max (compare with Fig. 8(b)). At τ≈0.7155, x min 2 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 τ≈1.208, x min 1 undergoes AH bifurcation and the stable cycle is born from it, which for a slightly larger τ 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 τ≈1.53, the 1st homoclinic bifurcation occurs, namely, the limit cycle around x min 2 collides with x 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 τ increases slightly beyond the value of the 1st homoclinic 1.53, the local attractor on the right side of x max , i.e. near x min 2 , ceases to exist. The system has a single attractor left, which is the limit cycle localised near x min 1 (compare red lines in Fig. 10(d) and Fig. 11(o)). If at τ just above 1.53 the initial conditions are set on or near the just disappeared local attractor around x min 2 , the phase trajectory swiftly leaves this Homoclinic bifurcation with a saddle−node in 2D with a saddle−node in 3D with a saddle−focus in 3D Homoclinic bifurcation Homoclinic bifurcation with a saddle−focus in 3D  with a saddle−focus in 3D 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 max returns to the vicinity of x 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 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 E. 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)).
12. Schematic illustration of homoclinic bifurcations in some fictional (a)-(e) two-dimensional and (f)-(j) three-dimensional systems, respectively, with a saddle-node, which are similar to the 2nd and 3rd homoclinic bifurcations of the saddle-focus in (2), (9) with safe loops forming as τ is increased from 2.425 to 2.46 as shown in Fig. 9. Notations are as in Fig. 11, and in addition in (b)-(e) brown line shows a typical phase trajectory. Away from the saddle-node, the stable manifold (cyan line or surface) behaves in a qualitatively the same manner as with a saddle-focus (compare (f)-(g) with Fig. 13). In (b) and (g) the 2nd homoclinic bifurcation is illustrated, at which a large safe loop (magenta line) is formed by the manifolds of the saddlenode (compare with Fig. 10(f) and Fig. 13(d)). This bifurcation leads to the birth of a large limit cycle (larger red closed curve in (c) and (h), compare with turquoise curve in Fig. 10(h)). In (d) and (i) the 3rd homoclinic bifurcation is illustrated, at which a small safe loop (magenta line) is formed by the manifolds of the saddlenode (compare with Fig. 10(j)). This bifurcation leads to the disappearance of the small limit cycle and leaves only one attractor being the large limit cycle (red line in (e) and (j), compare with Fig. 10(l)).

E. 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 min 2 (upper right part of the panel), whereas the basin of attractor at or near x min 1 (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 ϕ(t), t ∈ [−τ, 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 Importantly, numerical simulation of (2), (9), 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 τ 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), (9) are given in Fig. 14 for a range of τ 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 τ. Below the bifurcation point τ=1.53, in addition to infinity, there are two attractors on different sides of x max , and on the plot we do not distinguish between fixed points and limit cycles. Due to the 1st homoclinic bifurcation at τ≈1.53, the basin of attractor at or near x min 2 (red shade) collapses, and at 1.53<τ<2.43 besides the infinity there is only one attractor at or near x min 1 (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 max makes several turns around the three fixed points (compare with Fig. 11(b)-(c)). FIG. 14. Basins of attraction of (2), (9) with safe homoclinic loops, as a function of τ. Different shades mark points x such that, from initial conditions ϕ(t)=x, t∈[−τ, 0], the phase trajectory converges to the following attractors: infinity (blue), fixed point at, or limit cycle around, x min 1 (green), fixed point at, or limit cycle around, x min 2 (red). Above the 1st homoclinic bifurcation point τ=1.53 there is no local attractor at or near x min 2 . Compare the stripy structure of the basins below τ=1.53 here and in Fig. 11(a)-(c).

F. Second homoclinic bifurcation
As mentioned in Section IV B, the 2nd homoclinic bifurcation cannot be predicted based on the theorems overviewed in Section II B. 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 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 min 1 (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), (9), 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 max and then towards the local attractor near x min 1 (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 τ≈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 τ 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), (9) is given by turquoise line in Fig. 10(h).

G. Third homoclinic bifurcation, chaos and infinity
The 3rd homoclinic bufurcation takes place at τ=2.4499 (see Fig. 8(b)) and eliminates the small limit cycle to the left of x max , i.e. around x min 1 , 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 τ>2.4499 there are no more local attractors in the system.
At τ∈(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 τ 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 τ 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 τ 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 τ e.g. when periodic windows in chaos appear 70 ). However, when τ 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 τ>0. 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.

H. Different forms of homoclinics
We also studied the phenomena induced by the increase of τ in a slightly different system of the form (2) with a doublewell 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 max is past AH bifurcation. There, the homoclinic orbit is formed by the manifolds of the saddle cycle born from x max , which become tangent to each other, as described in Section II C. Nevertheless, the general sequence of bifurcations there is very similar to the one in (2), (9).

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 twowell 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 τ as control parameter is illustrated in Fig. 8(c). To obtain this figure, Eqs. (2), (9) were solved as τ was decreased from some positive value to zero. Namely, at t=0, the system was launched from initial conditions belonging to the large chaotic attractor at τ=3.3. Then τ 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 τ was almost adiabatic.
With this procedure, at τ=0 the system ended up at the minimum x min 1 , at which the modulus of the Jacobian, |J|, is smaller than at x min 2 (see Fig. 8(a)), and which therefore undergoes AH bifurcation at a larger value of τ (see Fig. 8(b)). This minimum x min 1 is not only flatter than x min 2 , but also broader, due to which the local attractor around x min 1 survives for larger values of τ than the one around x min 2 . Under the adiabatic decrease of τ within approximately [1.53, 3.3], i.e. above the 1st homoclinic, at every instantaneous value of τ the solution is close to the only attractor corresponding to this value of τ. The latter is either the one embracing two minima at larger τ, or the one localised near, or at, x min 1 at smaller τ. When τ decreases below the 1st homoclinic at τ≈1.53, the local attractor around x min 2 appears. However, the solution is not affected by this event and remains in the vicinity of x min 1 until τ reaches zero. Since in this example x min 1 happens to be the lowest minimum, adiabatic decrease of τ has demonstrated optimisation as required.
In Supplementary Note, Section S-II B, optimisation with adiabatic decrease of τ 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), (9). To summarise, adiabatic decrease of τ 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 τ 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 τ, at some randomly chosen time t the system can be found in one of these vicinities. If at this moment τ 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 τ. This way, at τ=0 the system could end up in the same well it was before τ started to decrease, to whose bottom it would then converge following standard gradient descent (1). Several repetitions of the same experiment with τ 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 specified as follows The functions V (x), f (x) and J(x)= f (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 min 1 , but the lowest minimum is x min 2 . First, we launch this system at t=0 from the initial conditions x(t)=2 for t ∈ [−2.1, 0], and τ=2.1, and start to find the numerical solution x(t) while decreasing τ at a relatively low rate. Namely, after 0.1 time units, τ 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 τ 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 τ the system ends up near the flattest and broadest minimum x min 1 . Next, we repeat the process of decreasing τ 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 min 2 . In Fig. 15(c) the same solutions are shown as in (b), but now against the current value of τ for the convenience of comparison. Thus, the fast decrease of τ 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 τ 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 τ 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 (9), 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 Here we show how the standard analytical tools can be applied to predict the local behaviour of the delay-differential equation (DDE) of the following forṁ where is some landscape function, and f (z) is a twice differentiable function. Specifically, by performing linear stability analysis, we will reveal at what values of time delay τ the fixed points at the minima of V (z) change from stable nodes to stable focuses and then undergo Andronov-Hopf bifurcation. Doing so will allow us to predict at what τ around the given fixed point a stable periodic oscillation could be observed. We also analyse the stability of the fixed points located at the maxima of V (z). Although these points remain unstable for any positive τ and thus unobservable in an experiment, perhaps counterintuitively, the parameters characterising their instability affect the observable behaviour of the system. For example, it is possible for the manifolds of an unstable fixed point of a saddle-focus type to close and to form a homoclinic loop. As as the parameter of the system changes, this loop breaks down and can give birth to a new attractor, which will extend well beyond the close vicinity of the relevant point. According to Shilnikov's theorem, this new attractor can be either periodic or chaotic, and its nature is determined by the local properties of the respective unstable point 55,56,S1,S2 .
We start from performing linear stability analysis of the fixed points x * of (S1), which satisfy Unlike a system of ordinary differential equations (ODEs), a DDE is an infinite-dimensional dynamical system because in order to solve such an equation, as an initial condition one needs to specify a function on an interval, namely, x(t)=ϕ(t) on t ∈ [−τ, 0], which contains infinitely many points. With this, if we wish to represent the dynamics of a DDE in the phase space by analogy with what is done for ODEs, we need to appreciate that the full state of a DDE at time t would a function on an interval [t − τ,t], and therefore the dimension of the respective phase space is infinity S3 . As a single state of a) E-mail: N.B.Janson@lboro.ac.uk the DDE, the fixed point can be viewed as a constant function x(t)=x * on an interval [t − τ,t], which stays the same at all times t. By analogy with the linear stability analysis of ODEs, we need to slightly perturb the system away from the fixed point. So we introduce where X(t) is the perturbation of the fixed point x * . Given thaṫ x =Ẋ, substitute the latter into (S1) to obtaiṅ where X τ (t)=X(t − τ). Because X is small, X τ 1, and we can approximate the above equation by expanding f in Taylor series and keeping the linear terms onlẏ where J= d f dz at z=x * . With account of (S2), we obtain the socalled linearized equation, which holds as long as X(t) stays smallẊ We assume that the linear DDE (S3) has exponential solutions, i.e. These two curves are guaranteed to intersect at some positive and real value of λ , which means that the local maxima of the landscape V (z) in (S1) are always unstable.  S2. Real (upper panel) and imaginary (lower panel) parts of Lambert Function W (z), assuming that z takes only real values. Red line shows the principal branch W 0 (z), magenta line shows W −1 (z), and blue line shows W 1 (z).

Im[W(z)]
With a positive J, Eq. (S4) always has a real positive root λ (see Fig. S1 illustrating the solution).
The solutions λ of (S4) can be expressed in terms of the Lambert function W (z), also called the product log function. Lambert function is a complex-valued function of a complex argument defined as an inverse of the function h(z)=ze z , where z, h ∈ C S5,S6 . The function W (z) is formed of solutions w of the equation so that for every value of z there is a countable infinity of the values of w. In applications to the stability of DDEs, only real arguments z of W (z) are used, which means that only the values of W (z) corresponding to the cross-section along Im[z]=0 are relevant to the problem. From now on, we will assume that z are real. The function W (z) with z∈R has countably many branches W k , k = 0, ±1, ±2, . . ., and W 0 is called the principal branch. In Fig. S2, real (upper panel) and imaginary (lower panel) parts of W (z) are shown as a function of z.
The solutions λ of (S4) are equal to where W k is the k-th branch of W (z). Note, that in (S6), for z > − 1 e the principal branch W 0 (z) takes real values only. In Fig. 7 eigenvalues λ k are illustrated as functions of (Jτ) for the fixed points at local maxima of the landscape V of (S1) with J>0. As confirmed by Fig. 7, such fixed points are always unstable because one of their eigenvalues is real and positive at all positive values of τ. Note, that at τ=0 these fixed points have only one positive real eigenvalue equal to J. Indeed, at τ=0 (3) becomes a first order ODE, in which fixed points have only one real eigenvalue equal to the derivative J of the righthand side. Fig. 6 shows eigenvalues as functions of (|J|τ) for the fixed points at local minima of V with J<0. One can see that for small τ these fixed points are stable, but as τ increases, they become unstable.
Next, we calculate the values of τ, at which Andronov-Hopf (AH) bifurcation occurs for the given fixed point. If a stable fixed point undergoes AH bifurcation, the point becomes unstable and possibly a stable periodic solution is born from it, thus giving rise to oscillatory behaviour of the system. The first condition of AH bifurcation is that the real part of W 0 (Jτ) must be equal to zero. Given J, we need to determine the respective value of τ at which (S4) has purely imaginary roots, i.e. λ =±β i with β ∈R. Substituting this into (S4) gives This can be written as Separating real and imaginary parts gives Equation (S7) has infinitely many solutions in the form of conjugate pairs, which we then substitute into (S8) to obtain τ = − π 2 + nπ J(−1) n , n = 0, 1, 2, . . . , (S10) which will have the same τ for both components of the conjugate pair. Assume that we consider the local minima of the landscape V (z), at which J<0. Also, in this work we assume that τ can only be positive, and to ensure that this is the case in (S10) n must be even. The first Andronov-Hopf bifurcation occurs for the smallest even value of n, therefore we take n=0 to obtain . (S11) This result is consistent with the earlier predictions, obtained with a very different approach, of the existence of periodic solutions in DDEs (4) for a special form of g(z) as discussed in Sec. IIA. At larger even values of n in (S10), more complex conjugate pairs of eigenvalues cross the imaginary axis from below.
The second criterion to be satisfied for an Andronov-Hopf bifurcation to occur is S4,S7 To verify this, differentiate both parts of (S4) with respect to τ to obtain (S13) Assume that λ =α+iω with α, ω ∈ R, and substitute into (S13). Separating real and imaginary parts gives Eq. (S12) implies dα dτ τ=τ AH > 0. (S15) At AH bifurcation α=0 and ω=β . By setting in (S9) n=0 and substituting τ from (S11) we obtain Substituting these values into (S14) gives (S17) Combining these equations to exclude dω dτ gives So we have verified both criteria for the occurrence of AH bifurcation, in which the fixed point loses stability. In the numerical simulations performed here, the bifurcation is always supercritical and leads to the birth of a stable limit cycle. Equation (S10) can also be used to find the values of τ at which the fixed points at the maxima of the landscape V (z) undergo AH bifurcation. In this case J≥0, therefore, to ensure the positivity of τ, n must be odd. Again choosing the smallest value of n, n=1, gives τ s AH = 3π 2J . (S19) Such fixed points would remain unstable for any τ. However, as a result of this bifurcation, a saddle periodic orbit could be born from the fixed point, together with its stable and unstable manifolds, which will affect the dynamics of the system. Further complex conjugate pairs will cross the imaginary axis for odd values of n in (S10). This theory can be extended to equations of the form (2) (or (S1)) with an n-dimensional vector x. The approach will be similar, however, an n×n Jacobian matrix will be involved, resulting in an order-n quasi-polynomial for the characteristic equation, which can be analysed to deduce the relevant bifurcations.

S-II. HOMOCLINIC ORBIT TO A SADDLE CYCLE IN DDE (S1) WITH TWO-WELL POTENTIAL
In an example below, we illustrate a chain of global homoclinic bifurcations in a DDE of the form Eq. (S1) (same as Eq. (2) in main paper) with a two-well potential V , in which all homoclinic bifurcations come in a different form as compared to the ones in Eqs. (2), (9). Namely, in (2), (9) all three global bifurcations occurring as the delay τ grows consist in the formation of homoclinic loops of the saddle-focus fixed point at the maximum of V . However, here some of these manifest themselves in the formation of homoclinic orbits to the saddle cycle. These homoclinic orbits arise as the manifolds of the saddle cycle, born from the saddle fixed point at the maximum, become tangent to each other. Although the latter bifurcations are somewhat more complex than the ones involving the fixed point, they can nevertheless be routinely expected to occur in the systems of the form (2).
With this, the locations, depths and shapes of both minima stay roughly the same as in (9) to make the results compara- (b) Bifurcation diagram of (2) (or (S1)), (S20). Local minima/maxima of attractors are shown by black/red dots. Fixed points at the minima of V are shown by red vertical lines for τ at which they are stable, and the saddle-focus at the maximum of V is shown by green vertical line. As τ varies, homoclinic orbits are formed by the manifolds of the saddle cycle around x max born at the point AH1. (c) Demonstration of optimisation. Local maxima (red dots) and minima (black dots) of solutions to (S1), (S20), which are obtained as τ slowly decreases from 2.8 to zero. The solution spontaneously settles down at the lowest, flattest and broadest minimum at x min 1 . Compare with (a) and (b).
ble. The functions V (x), f (x) and J(x)= f (x) are shown in Fig. S3(a) by blue, red and green lines, respectively. This way, in (S20) the Jacobian of the saddle-focus fixed point at the maximum x max ≈−0.0420246 is J(x max )≈4. 54 and is much larger than 0.49 in (9). As a result, its AH bifurcation occurs at a much smaller value of τ≈1.037783 as compared to 9.61 in (2), (9).  Fig. S3(b) illustrating bistability, hysteresis, and the 1 st homoclinic bifurcation in (2) (or (S1)), (S20) eliminating the localised attractor inside the right well of V (compare with Fig. 5). The version of the homoclinic bifurcation here manifests itself in the homoclinic orbit formed by the tangency of the manifolds of the saddle cycle around the potential maximum x max . Green line shows x max . Dots show local minima of attractors: orange dots show limit cycle born via the first AH bifurcation from x min 2 (outside the range of this figure), and blue dots show attractor developed from the homoclinic orbit of the saddle cycle via 1 st homoclinic bifurcation as τ decreases from 1.118. AH1 is the first AH bifurcation of x max , which gives rise to the saddle cycle together with its manifolds at τ≈1.118.
tors coexist: the only remaining local attractor around x min 1 (blue line in Fig. S5(i)) and the large attractor (cyan line in Fig. S5(i)).
At τ≈2.4968337 the 3rd homoclinic bifurcation occurs leading to the disappearance of the last of the local attractors. At higher τ, the only remaining attractor is the large chaos (cyan line in Fig. S5(l)), which survives until τ≈2.87 (blue line in Figs. S5(m)-(n)). Like in (2), (9), as τ reaches a certain threshold (here 2.876), even this large attractor vanishes and the phase trajectory goes to infinity (blue line in Fig. S5(o)).

B. Optimisation
One possible optimisation procedure is illustrated in Fig. S3(c). To obtain this figure, Eqs. (2) (or (S1)), (S20) were solved with the parameter τ being slowly decreased in a step-wise manner from the value 2.8 to zero, after the system was launched from initial conditions belonging to the chaotic attractor at τ=2.8. The procedure was the same as the one described in Section V of the main paper.
The system ended up at the minimum x min 1 , at which the modulus of the Jacobian, |J|, is smaller than at x min 2 (see Fig. S3(a)), and which therefore undergoes AH bifurcation at a larger value of τ (see Fig. S3(b)). This minimum x min 1 is not only flatter than x min 2 , but also broader, due to which the local attractor around x min 1 survives for larger values of τ than the one around x min 2 . For this reason, for the current shape of V adiabatic decrease of τ delivered the global minimum.
An example with the lowest minimum not being the flattest or belonging to the broadest well is considered in the main paper, Section V, where it is shown how such a minimum could FIG. S5. Phase portraits of (2) (or (S1)), (S20) at various τ illustrating bifurcation diagrams in Figs. S3 and S4 with τ given in panels. Notations are: fixed points x min 1,2 (red circles) and x max (green circle); filled/empty circles indicate points below/above AH bifurcations; attractors (blue, turquoise and orange lines). be found with the fast decrease of τ.

S-III. ADDITIONAL EXAMPLES WITH SIMILAR BEHAVIOUR
In the main paper we provided a detailed illustration of a chain of homoclinic bifurcations occurring in a system of the form (2) (or (S1)) with a double-well potential described by (9). In Section S-II we illustrated a similar chain with a slightly different form of homoclinic bifurcations in a subtly different system whose potential function is given by (S20). In addition, we observed qualitatively the same bifurcations in (2) with slightly different versions of double-well potentials V , which we do not illustrate here. Three of these examples are given below.