Border-collision bifurcations in a driven time-delay system
BBorder-collision bifurcations in a driven time-delay system
Pierce Ryan, Andrew Keane,
2, 1 and Andreas Amann School of Mathematical Sciences, University College Cork, Cork T12 XF62,Ireland Department of Mathematics, The University of Auckland, Private Bag 92019, Auckland 1142,New Zealand
We show that a simple piecewise-linear system with time delay and periodic forcing gives rise to a richbifurcation structure of torus bifurcations and Arnold tongues, as well as multistability across a significantportion of the parameter space. The simplicity of our model enables us to study the dynamical featuresanalytically. Specifically, these features are explained in terms of border-collision bifurcations of an associatedPoincaré map. Given that time delay and periodic forcing are common ingredients in mathematical models,this analysis provides widely applicable insight.
Both time-delay dynamical systems, and periodi-cally driven dynamical systems have been thor-oughly studied in the literature. This can beattributed to their great relevance to real-worldproblems. Time-delay systems arise naturally inphysical, biological or climate models due to finitepropagation speed; periodic drive is ubiquitous inengineering applications and is known to generatecomplex resonance phenomena. However, sys-tems that combine these two properties have re-ceived much less attention, despite being relevantin many real-world applications. In this paper,we study a simple piecewise-linear system withboth time delay and periodic forcing, which ex-hibits interesting dynamical features as a nontriv-ial consequence of this combination. These fea-tures include multistabilities, Arnold tongues andtorus bifurcations. Since the system is piecewiselinear and contains only two parameters, manyphenomena can be interpreted through an ana-lytically derived piecewise-smooth Poincaré mapand an analysis of the associated border-collisionbifurcations. The analysis explains the origin ofsimilar phenomena which has previously been ob-served numerically in more complicated relatedsystems.This article may be downloaded for per-sonal use only. Any other use requiresprior permission of the author and AIP Pub-lishing. This article appeared in Chaos30, 023121 (2020) and may be found athttps://aip.scitation.org/doi/10.1063/1.5119982.
I. INTRODUCTION
Periodically driven systems appear in many real-worldapplications. Examples include optical injection in lasersystems , vibration-driven energy harvesting devices ,injection-locked frequency dividers in electronics , or sea-sonal forcing in climate systems . They often give rise tointeresting resonance behaviour in damped oscillators and complex synchronization patterns in self-sustained oscillators .Similarly, time-delay systems also arise in many ex-perimental systems, for example in optics , electronics ,neuro-science , or climate systems , and also play animportant role in chaos control, for example, through theuse of time-delayed feedback control . From a mathe-matical point of view, time delay often leads to a for-mally infinite-dimensional phase space , which consid-erably complicates the analysis, but allows for a rich va-riety of phenomena.The combination of external forcing and time delay hasbeen studied, for example, in the context of the Duffingoscillator , the van der Pol oscillator and more re-cently in the context of climate systems . However,a general understanding of this class of dynamical sys-tems is not yet available. The objective of the currentpaper is, therefore, to study the fundamental features ofan elementary system with time delay and periodic forc-ing to obtain a broader insight into what phenomena areexpected to arise as a consequence of this combination.Let us consider a simple driven time-delay dynamicalsystem introduced by Ghil et al. as a model for a cli-mate phenomenon known as the El Niño Southern Os-cillation. The system of a real variable x ∈ R is definedby ˙ x ( t ) = − tanh [ κx ( t − τ )] + b sin (2 πt ) , (1) x ( t ) ∈ C ([ − τ, , (2)where b ≥ is the magnitude of the periodic forcing, τ > is the time delay of the delayed feedback, and κ > is the linear slope of the delayed feedback at the origin.A solution of the system is a trajectory x ( t ) ∈ R , t ∈ R .A consequence of the reliance of the delayed feedback ona continuous function x ( t ) over an interval [ − τ, is thatthe system has an infinite-dimensional phase-space. An-other key feature of this system is that it has the sym-metry x ( t ) → − x (cid:0) t + (cid:1) . This model has been stud-ied extensively by Keane et al. . Numerical analysisof this system demonstrated an extremely complex reso-nance structure . Furthermore, the autonomous system ( b = 0) has been studied analytically . In this casethe trivial solution x ≡ is only stable for τ < π κ . At τ = π κ , it becomes unstable, and a family of stable τ -periodic solutions is born. a r X i v : . [ n li n . C D ] F e b In order to analyse the phenomena seen in this modelfurther, let us consider a further simplification of the sys-tem by taking κ → ∞ . This has the effect of chang-ing the delayed feedback term from − tanh [ κx ( t − τ )] to − sgn [ x ( t − τ )] . We also apply the signum function tothe periodic forcing to obtain the dynamical system ˙ x ( t ) = − sgn [ x ( t − τ )] + b sgn (sin (2 πt )) , (3) x ( t ) ∈ C ([ − τ, , (4)where b ≥ and τ > . Critically, this simplifi-cation of the system preserves the symmetry x ( t ) →− x (cid:0) t + (cid:1) . The feedback term − sgn [ x ( t − τ )] takesvalues in { , , − } . However, while the feedback canin principle be , this occurs only under highly specificconditions which are not considered here. The forcingterm b sgn [sin (2 πt )] takes values in { b, , − b } . As theforcing is only at discrete times, we say that the forc-ing is positive for t mod 1 ∈ [0 , . , and negative for t mod 1 ∈ [0 . , . We now develop our method for solv-ing the system. II. NUMERICSA. Iterative map
In order to solve Eqs. (3,4), we note that ˙ x ( t ) canonly take discrete values in { b, − b, − b, − − b } ;therefore this continuous system can be modelled exactlyas a discrete time iterative map, or Poincaré map. Thestate of the system at time t is a tuple of variable length S ( t ) = ( x ; z , z , ..., z n − ) , (5)where x ∈ R is the position at time t and t − τ < z We now present sample solutions which demonstratesome of the characteristic features of this system. Fig. 1shows eight stable solutions found by simulating thesystem from different initial conditions and parameters ( b, τ ) .There are periodic and aperiodic solutions present inthe system. A periodic solution with period P is a tra-jectory that follows a cycle such that S ( t + P ) = S ( t ) .An aperiodic solution is considered to be a solution withinfinite period. We label solutions by the characteristicratio P : R , where R is the number of times the trajectorycrosses from x < to x > in one cycle. Note that the P : R notation differs from the p : q notation used in someprevious literature , where pq is the rotation number.We find that P : R is a useful measure of a solution, forreasons that will be made clear later.Fig. 1(a) is an example of the stable solution to theunforced system ( b = 0) . The change in the feedbackoccurs at a time t + τ after the trajectory passes through x = 0 at time t , resulting in a τ -periodic solution thatis stable for τ > . This is consistent with the analyticresults found for the unsimplified system. We consider τ to be the natural period of the feedback, as the solutionto the unforced system is τ -periodic. The characteristicratio of this solution is τ : 1 .Fig. 1(b) and Fig. 1(c) show a solution and a solution, respectively, that are stable for the sameparameters. This is an example of bistability, where thereare two stable solutions for the same parameters; thesolution to which the system converges depends on whichsolution’s basin of attraction the initial conditions are in.A second example of bistability is seen in Fig. 1(d) andFig. 1(e), which show a solution and a solution,respectively, that are stable for the same parameters.The solution in Fig. 1(f) is assumed to be aperiodic, asthe system does not converge to a periodic solution afterrunning a simulation up to T = 100000 . By comparingthe aperiodic solution to the solution in Fig. 1(g), wemay note that the aperiodic solution and the solu-tion have the same number of x = 0 crossings per periodof the forcing. However, the x = 0 crossings are evenlyspaced in the solution, but are not in the aperiodic x ( t ) b =0.00=0.60( a ) 4 :1101 x ( t ) b =1.35=0.95( b ) 4:2101 x ( t ) b =1.35=0.95( c ) 3:1101 x ( t ) b =2.30=1.17( d ) 1:1202 x ( t ) b =2.30=1.17( e ) 5:1202 x ( t ) b =4.10=1.51( f ) :202 x ( t ) b =4.10=1.45( g ) 1:10 1 2 3 4 5 6 7 8 9 10 t x ( t ) b =2.65=0.93( h ) 3:3 Figure 1. Sample solutions, excluding transients, obtainedfrom simulating the system from different initial conditionsand parameters ( b, τ ) shown in the bottom right corner of eachplot. The ratio of the period of the solution to the numberof crossings from x < to x > in one cycle is shown in thetop left corner of each plot. The vertical dotted lines indicatetimes when the forcing changes. solution. This may be an indicator that the aperiodicsolution is related to the solution in Fig. 1(g). Asimilar observation can be made for the solution inFig. 1(h) and the solution in Fig. 1(d). In later sec-tions of this paper, we will investigate these relationshipsin detail. C. Structure in the ( b, τ ) plane Having seen some interesting features of the system,we move on to understanding the overall dynamics inthe ( b, τ ) plane. Due to the bistability present in thesystem, simulating the system from arbitrary initial con- ditions across a ( b, τ ) mesh would produce an inconsistentpicture, as we have no prior knowledge of the basins ofattraction of bistable solutions. In order to circumventthis issue, for fixed τ , the solution is swept across a rangeof b by iteratively simulating the system, incrementing b slightly, then simulating again using the final state of theprevious simulation as the initial state of the next one.This allows a stable solution to be followed until it losesstability or ceases to exist, at which point the systemconverges to a nearby stable solution. By taking multi-ple sweeps in b for a range of fixed τ values, we obtain the P : R charts shown in Fig. 2. Fig. 2(a) shows the P : R chart under an upward sweep in b , from left to right.Fig. 2(b) shows the P : R chart under a downward sweepin b , from right to left. This figure demonstrates manystriking features of the system which will be explored inmore detail.First, let us focus our attention on the upward b sweepin Fig. 2(a). We observe that there are regions in the ( b, τ ) plane in which particular P : R solutions exist,such as the labelled , , and regionson the left side of the chart. These regions are sec-tions of Arnold tongues. An Arnold tongue is a regionof the ( b, τ ) plane, rooted on b = 0 , within which thefeedback and the forcing synchronise to produce a so-lution with period equal to a rational ratio of the forc-ing. In piecewise-linear systems, an Arnold tongue canhave shrinking points, at which the Arnold tongue haszero width. This phenomenon was first observed in thecircle map , and analysed in detail in the context ofpiecewise-linear continuous maps with single switchingmechanisms . This produces Arnold tongues thathave been compared to strings of sausages . We willrefer to an individual “sausage” as a tongue , and refer toa full “string” as an Arnold tongue . A notable featureof these Arnold tongues is that the characteristic ratio isdifferent in each tongue in the string. For example, the P = 7 Arnold tongue is rooted on b = 0 at τ = 1 . , andconsists of the tongue, the tongue, the tongue,and the unlabelled tongue. In each Arnold tongue,the leftmost tongue is a P : R tongue that is rooted on b = 0 at τ = P R . P is the same in every tongue in thechain, but R increases the further right the tongue is inthe string.Tongues with similar characteristic ratios tend to havesimilar shape, with some variation. For example, the and tongues have identical shape; the tongue isdifferent. The large P : 1 tongues noted earlier have iden-tical shape for P > . The Arnold tongue is unlikeany other Arnold tongue in shape. For large b , the forc-ing dominates the feedback, which results in the systemconverging to the solution exclusively. We refer tothe region in which b is large enough as the locked region,where the system is locked to the period of the forcing.The boundary of the locked region is unusual, being madeup of straight lines which meet at right angles. Branch-ing off from the horizontal lines of the boundary, thereare vertical stripes in which P : P solutions exist. We b ( a ) b ( b ) Figure 2. Colour maps showing the period of simulated solu-tions obtained by making 1124 sweeps of 1125 simulationsof duration T = 10000 in [0 , × [0 , . . The white ar-rows indicate the direction of the sweep. Dark blue indicates P = 1 solutions, shading to dark red for P = 999 solutions.White space indicates where the solution was aperiodic, orhad P > . The white text indicates the characteristic ra-tio of stable solutions found within the tongues. The dashedwhite line shows the border between the regions where the and solutions occur. will devote considerable attention to studying the solution later.Now we compare the upward b sweep in Fig. 2(a) to the downward b sweep in Fig. 2(b). The most significantdifference occurs around ( b, τ ) = (2 . , . . In the up-ward sweep, the system follows the solution to theedge of the tongue, and there is a complicated re-gion of smaller tongues and aperiodicity above the tongue. In the downward sweep, the system instead con-tinues to follow the solution into the region wheredifferent features occurred in the upward sweep. Thisagrees with the bistability of the and solutionsseen in Fig. 1(d,e). A less obvious difference is the bista-bility to the right of b = 1 ; note the apparent differencein shape of the P : 2 and P : 1 tongues near this line.This occurs because the P : 1 solutions overlap with the P : 2 solutions, in agreement with the bistability of the and solutions seen in Fig. 1(b,c). The absenceof P : 1 tongues for even P is notable. We numericallyobserve such solutions in simulations, but only for small b (cid:62) . and exactly τ = P . III. DYNAMICS Fig. 3 shows maximum charts in the ( b, τ ) plane over-layed with bifurcation curves. The maximum is takenas the maximum value of a solution over an interval oflength after a transient of length from the samesimulations that were used to generate Fig. 2. Some ofthe larger tongues seen in Fig. 2 can be seen withoutdifficulty in Fig. 3 due to a large jump in maxima atthe boundaries of the tongues. We note that the even P : R tongues are striped horizontally; this is most evi-dent in the and tongues. It appears that so-lutions with even P or R are not invariant under thesymmetry x ( t ) → − x (cid:0) t + (cid:1) ; rather there exists a pairof symmetry-related counterpart solutions, each with adifferent maximum value. The system converges to oneof these solutions depending on initial conditions, andremains at that solution until swept out of the tongue,resulting in stripes parallel to the direction of the sweep.This feature was also observed in the smooth system (1,2)by Keane, Krauskopf, and Postlethwaite .We will devote the rest of this section to deriving andcharacterising the bifurcation curves plotted in Fig. 3.In order to do so, we first need to establish a systematicmethod of analysing solutions. We apply this method todevelop Poincaré maps and border collision maps throughwhich we study the bifurcations present in this system. A. Symbolic representation We require a robust framework under which we cananalyse the dynamics of this system. We note that thecharacteristic ratio does not distinguish between the twodifferent forms of the solution shown in Fig. 1(d,g).Observing the order in which the feedback and forcingchange after the trajectory passes through x = 0 , we notethat the feedback changes before the forcing in Fig. 1(d), b ( a ) b ( b ) Figure 3. Maximum charts in the ( b, τ ) plane overlayedwith bifurcation curves. Blue indicates a low maximum, redindicates a high maximum. The black arrows indicate thedirection in which solutions were swept. The black text in-dicates the ratio of the period of the solution to the numberof Z symbols in the sequence for the stable solution withinthe associated tongue. BCSN bifurcations are shown in solidwhite, and T bifurcations are shown in solid black. The D ¯ D curve is shown in dashed white. and after the forcing in Fig. 1(g). In order to preciselycapture these differences, a more explicit labelling systemis required. We therefore adopt a symbolic representa-tion of solutions. Symbolic representations have previ- -2 -1 0 1 2 3 4 t x ( t ) Z D D HD D D Z D D HD D D Z Db =2.30=1.17 ( a ) t x ( t ) Z D D DH D D Z D D DH D D Z Db =2.30=1.33 ( b ) Figure 4. Derivation of the symbolic representations of the solution. ously been used to great effect in the study of iterativemaps .Let a solution be represented by a sequence of events ...X X ...X n ... where X i ∈ (cid:8) D, ¯ D, Z, ¯ Z, H, ¯ H (cid:9) . D de-notes a transition of the forcing from − b to b , Z de-notes a transition from x < to x > , and H de-notes a transition of the feedback from − to . A barover a symbol causes it to denote the opposite transi-tion; a symbol with two bars over it is the same asthe symbol unbarred. Fig. 4(a) shows the solu-tion seen Fig. 1(e), labelled with the events that occurin the trajectory. As this solution is periodic, the se-quence repeats, so we abbreviate the sequence of eventsrepresenting the solution to a minimal repeating sequence (cid:2) Z, ¯ D, D, ¯ H, ¯ D, D, ¯ D, ¯ Z, D, ¯ D, H, D, ¯ D, D (cid:3) . In general, a P : R solution is represented by a minimal repeating se-quence of events [ X , X , ..., X n ] containing: • P D events and P ¯ D events, • R Z events and R ¯ H events, • R ¯ Z events and R H events.Every cyclic permutation of a sequence represents thesame solution. For the sake of consistency, all sequencesbegins with a Z . We observe that the solutionshown in Fig. 4(a) is invariant under the symmetry x ( t ) → − x (cid:0) t + (cid:1) ; this causes the second half of the se-quence (cid:2) Z, ¯ D, D, ¯ H, ¯ D, D, ¯ D, ¯ Z, D, ¯ D, H, D, ¯ D, D (cid:3) to bethe same as the first half with all symbols barred. Weabbreviate (cid:2) Z, ¯ D, D, ¯ H, ¯ D, D, ¯ D, ¯ Z, D, ¯ D, H, D, ¯ D, D (cid:3) as (cid:2) Z, ¯ D, D, ¯ H, ¯ D, D, ¯ D (cid:3) − for convenience. In general, a P : R solution of the form (cid:2) X , X , ...X n , ¯ X , ¯ X , ..., ¯ X n (cid:3) canalso be represented by a half-sequence [ X , X , ...X n ] − .A sequence is considered legal within a subset of the ( b, τ ) plane if it represents a solution that exists withinthat subset. Each P : R solution is represented by a setof legal sequences, each one existing in a unique subsetof the ( b, τ ) plane. The union of these subsets is the re-gion in which the solution exists. The solution existswithin the tongue shown in Fig. 2(a). It is repre-sented by the half-sequences (cid:2) Z, ¯ D, D, ¯ H, ¯ D, D, ¯ D (cid:3) − for τ ≤ . and (cid:2) Z, ¯ D, D, ¯ D, ¯ H, D, ¯ D (cid:3) − for τ ≥ . , asshown in Fig. 4. At τ = 1 . , the increasing time delaybetween Z and ¯ H causes ¯ H to swap with ¯ D , changing thesequence. Determining whether a given sequence is legalwithin a subset of the ( b, τ ) plane is a nontrivial prob-lem. Appendix A presents a general method to determinewhether a sequence is legal for a given ( b, τ ) , which canalso be used to calculate a solution analytically from asequence for a given ( b, τ ) . This method was used to plotthe solutions in Figs. 4,5 and the bifurcation curves inFig. 3. B. Torus bifurcation of the solution We now apply our sequence representation in analysingthe solution; specifically, we determine what happensat the vertical black lines along the boundary of thelocked region in Fig. 3. Consider the set of half-sequencesthat represent the solution. Each half-sequencecontains one Z or ¯ Z , one H or ¯ H , and one D or ¯ D . Weneed only consider half-sequences starting with Z , as [ X , X , X ] − = [ X , X , X ] − = [ X , X , X ] − . Thereare only eight possibilities: [ Z, H, D ] − , [ Z, D, H ] − , (cid:2) Z, ¯ H, D (cid:3) − , (cid:2) Z, D, ¯ H (cid:3) − , (cid:2) Z, H, ¯ D (cid:3) − , (cid:2) Z, ¯ D, H (cid:3) − , (cid:2) Z, ¯ H, ¯ D (cid:3) − , (cid:2) Z, ¯ D, ¯ H (cid:3) − . [ Z, H, D ] − and [ Z, D, H ] − arenot legal, as they require Z to occur when ˙ x < , whichis impossible. Similarly, (cid:2) Z, ¯ H, D (cid:3) − and (cid:2) Z, D, ¯ H (cid:3) − arelegal only for b < , and (cid:2) Z, ¯ D, H (cid:3) − and (cid:2) Z, H, ¯ D (cid:3) − arelegal only for b > . The solutions shown in Fig. 1are (cid:2) Z, ¯ H, ¯ D (cid:3) − in Fig. 1(d) and (cid:2) Z, ¯ D, ¯ H (cid:3) − in Fig. 1(g).In the case b > , the sequence representing the solution is restricted by τ in the following way: (cid:2) Z, ¯ H, ¯ D, ¯ Z, H, D (cid:3) for τ mod 1 ∈ [0 , . , (cid:2) Z, ¯ D, ¯ H, ¯ Z, D, H (cid:3) for τ mod 1 ∈ [0 . , . , (cid:2) Z, H, ¯ D, ¯ Z, ¯ H, D (cid:3) for τ mod 1 ∈ [0 . , . , (cid:2) Z, ¯ D, H, ¯ Z, D, ¯ H (cid:3) for τ mod 1 ∈ [0 . , . (7)This can be explained by observing that for small τ , ¯ H must follow the Z that created it almost immediately.As τ increases, ¯ H drifts further away from Z in the se-quence, drifting past ¯ D at τ = 0 . . At τ = 0 . , ¯ H drifts past the subsequent ¯ Z . At τ = 0 . , ¯ H drifts pastthe subsequent D . At τ = 1 , the ¯ H created at Z driftspast the subsequent Z , and the pattern repeats. At thesame time, the same drift pattern occurs between ¯ Z and H . We now apply this information to analyse the solution. We construct a Poincaré map on the solution. Let t z be a time on the solution at which x = 0 ; thenthe state of the system at time t z is S ( t z ) = (0; z , z , ..., z n − , t z ) T . (8)As x = 0 at t = t z , we drop x and rewrite the state as S z = ( z , z , ..., z n − , t z ) T . (9)Let t ∗ z be the time on the solution at which x = 0 immediately after t z . Then the state of the system attime t ∗ z is S ∗ z = ( z , z , ..., t z , t ∗ z ) T . (10)As the solution is invariant under the symmetry x ( t ) → − x (cid:0) t + (cid:1) , it has the property S z = z z ...z n − t z = z − z − ...t z − t ∗ z − = S ∗ z − . (11)We define a Poincaré map P : S z → S ∗ z − so that the solution is a fixed point of P . Therefore, P is definedby P z z ...z n − t z = z − z − ...t z − t ∗ z − . (12)For the purpose of deriving P , we will assume w.l.o.g.that the trajectory is transitioning from x < to x > at t z ∈ [0 , . . Let t h ∈ [ t z , t ∗ z ) be the time at which thefeedback changes. First, we consider the case where τ < . ; then the zero element generated at t z is consumedat t h = t z + τ . Then P is one-dimensional as S z has onlyone variable, t z . We write the one-dimensional map as P ( t z ) = t ∗ z − . (13)If τ < . , then the solution is represented by the se-quence (cid:2) Z, ¯ H, ¯ D (cid:3) − ; therefore, the feedback is positivefor t ∈ [ t z , t h ) . We write an equation for the displace-ment of the trajectory between t z and t ∗ z as b (cid:18) − t z (cid:19) − b (cid:18) t ∗ z − (cid:19) +( t h − t z ) − ( t ∗ z − t h ) = 0 . (14)We substitute t h = t z + τ and solve for t ∗ z − to obtain P ( t z ) = − (cid:18) b − b + 1 (cid:19) t z + b + 2 τb + 1 − . (15)As (cid:12)(cid:12)(cid:12) b − b +1 (cid:12)(cid:12)(cid:12) < for b ≥ , P is stable for τ < .If τ ∈ [0 . , , then t h was generated, not at t z , butat the previous x = 0 crossing z . Then P is two-dimensional as S z has two variables, and the feedbackis negative for t ∈ [ t z , t h ) , where t h = z + τ . These con-ditions can be generalised for larger τ . The dimension ofthe system is n = (cid:100) τ (cid:101) , where (cid:100) τ (cid:101) is the smallest inte-ger greater than τ ; then the feedback for t ∈ [ t z , t h ) is ( − n − . We can then generalise Eq. (14) for arbitrary n and solve for t ∗ z − to obtain t ∗ z − 12 = − t z + b + 2 t h ( − n − b + ( − n − − . (16)For τ > , P can written as P ( S z ) = AS z + B, (17)where A is an n × n matrix A = · · · ... . . . . . . ...... ... . . . . . . 00 0 · · · − n − b +( − n − · · · · · · − , (18) B is given in Appendix B. Note that A only depends onthe parameter b and not on τ . This means that for fixed n , the value of of b at which the fixed point of P is bi-furcating, b bif , is constant. By solving the characteristicequation, we calculate b bif ( n ) = 1cos (cid:16) π ( n − n − (cid:17) − ( − n − , (19)where n = (cid:100) τ (cid:101) . Full calculations may be found in Ap-pendix B. The curve ( b bif ( (cid:100) τ (cid:101) ) , τ ) can be seen plottedas black vertical lines against maximum charts in Fig. 3 . If we examine the maximum chart where b is swept fromright to left, we see that the system ceases to convergeto the solution along this curve. We now know thatthis is because the fixed point of P , and hence the solution, loses stability at b = b bif . By calculating theeigenvalues of P explicitly for n = 2 and n = 3 , we findthat the loss of stability occurs because a pair of complexconjugate eigenvalues λ , cross | λ | = 1 . Therefore thefixed point of P loses stability due to a Neimark-Sacker(NS) bifurcation . As P is a Poincaré map on a periodicorbit, a NS bifurcation of the fixed point of P correspondsto a torus (T) bifurcation of the solution. However, P only shows the existence of the T bifurcation along thevertical sections of the boundary of the locked region,where n is constant and P is smooth. To fully under-stand the horizontal sections of the boundary, we mustlook to non-smooth bifurcation theory. IV. BORDER-COLLISION BIFURCATIONS The Arnold tongues seen in Fig. 2 have sharply de-fined boundaries, made up of curves and straight lines. We seek to determine what happens to solutions at theseboundaries and derive analytic expressions for the bound-aries using Poincaré maps.By simulating solutions near the boundaries of theArnold tongues, we observe that moving closer to theboundaries causes a D or ¯ D to move closer to x = 0 .An example of this is Fig. 4(a), where the D at t = 0 in the solution occurs for x > near x = 0 . Ifthis D crossed x = 0 and occurred at x < , this wouldsignificantly impact the feedback. There would be twoadditional x = 0 crossings in the trajectory, changingthe D to ¯ ZDZ and adding a H and a ¯ H elsewhere inthe sequence. Therefore, when we construct a Poincarémap to describe the dynamics of the system close to theboundary of such a tongue, the map must have a borderat x = 0 , such that the map is continuous across the bor-der but not differentiable at the border. Such maps, andthe associated border-collision bifurcations, have recentlyreceived systematic analysis in the literature . A. Border-collision saddle-node bifurcation of the and solutions We construct a Poincaré map B on the and solutions represented by (cid:2) Z, ¯ D, D, ¯ H, ¯ D, D, ¯ D (cid:3) − and (cid:2) Z, ¯ D, ¯ H, H, D, ¯ H, ¯ D, ¯ Z, D, Z, ¯ D (cid:3) − respectively, whichcan be seen in Fig. 5(a), such that our fixed point isthe state of the system at time t = 0 , at the position x D at which the D that will cross x = 0 occurs. As the ¯ H created at Z occurs before this D , B is one-dimensional.We divide B into B + and B − for x ≥ and x < respec-tively. The and solutions are invariant underthe symmetry x ( t ) = − x (cid:0) t + (cid:1) ; therefore we constructour map B : x D (0) → − x D (cid:0) (cid:1) . B + is a map on the so-lution represented by (cid:2) Z, ¯ D, D, ¯ H, ¯ D, D, ¯ D (cid:3) − . Followingthe blue curve in Fig. 5(b), we can derive B + ( x D ) = − (cid:18) b − b + 1 (cid:19) x D + 2 bb + 1 − b + 52 + 2 τ. (20) B − is a map on the solution represented by (cid:2) Z, ¯ D, ¯ H, H, D, ¯ H, ¯ D, ¯ Z, D, Z, ¯ D (cid:3) − , which is shown inFig. 5(b) in red. The feedback change due the trajec-tory dipping below x = 0 occurs at some t ∈ (1 , . hasduration − bxb − ; B − is otherwise identical to B + . Thus wederive B as B ( x D ) = − (cid:16) b − b +1 (cid:17) x D + bb +1 − b +52 + 2 τ x D ≥ (cid:16) bb − − b − b +1 (cid:17) x D + bb +1 − b +52 + 2 τ x D < (21)By setting x D = 0 and solving B for τ , we obtain thecurve τ = b +2 b +54( b +1) . This matches the lower right bound-ary of the tongue spanning from (1 , to (3 , . .For b ∈ [1 , and τ ∈ (cid:104) b +2 b +54( b +1) , . (cid:105) , B has two fixedpoints; a stable fixed point that exists for x D ≥ andan unstable fixed point that exists for x D ≤ . Thesetwo fixed points collide and vanish at the border x D = 0 at τ = b +2 b +54( b +1) in a border-collision saddle-node (BCSN)bifurcation . Hence the tongue is bounded by aBCSN bifurcation. Fig. 5(b,e,f) shows the BCSN bifur-cation of the and solutions.We generalise B for every P : 1 tongue for odd P ≥ as B P ( x D ) = − (cid:16) b − b +1 (cid:17) x D + bb +1 − b + | P − τ | x D ≥ (cid:16) bb − − b − b +1 (cid:17) x D + bb +1 − b + | P − τ | x D < (22)where | P − τ | is the term that causes the P : 1 tonguesto be symmetric across τ = P . By setting x D = 0 , wecalculate the boundaries of all P : 1 tongue for odd P ≥ and b ∈ [1 , as τ P ( b ) = P ± b − b b + 1) . (23)We find that this phenomenon persists throughout thesystem. The vast majority of tongues investigated arebounded by BCSN bifurcations occurring when a D crosses x = 0 . For b > , a stable P : R solution andan unstable P : R + 2 solution undergo a SN bifurcationwhen a D crosses x = 0 if the solution has the symmetry x ( t ) = − x ( t + P ) . Otherwise, the stable P : R solutionundergoes a SN bifurcation with an unstable P : R + 1 solution when a D crosses x = 0 .For b < , the mechanism by which BCSN bifur-cations occur is slightly different. Instead of a D crossing x = 0 and adding new symbols to the se-quence, D instead crosses x = 0 by swapping or-der with an existing Z . The stable solution (cid:2) Z, ¯ D, D, ¯ H, ¯ D, D, ¯ D (cid:3) − undergoes a BCSN bifurcationwith the unstable (cid:2) Z, D, ¯ D, D, ¯ H, ¯ D, D (cid:3) − at (cid:0) b, − b (cid:1) , b < , as shown in Fig. 5(a,c,d). While the unstable and solutions shown in Fig. 5 are two differ-ent solutions under our sequence representation, they areidentical at b = 1 , where a D passes through x = 0 with-out a bifurcation. A noteworthy feature of the BCSNbifurcations at the boundary of the tongue is thatthe pair of D events which collide at x = 0 is differentfor b < and b > , as seen in Figure Fig. 5(a,b,c,e).Generally, for b < , a stable P : 1 solution undergoes aBCSN bifurcation with an unstable P : 1 solution. Again,we produce a general form for the curves where a P : 1 solution undergoes a BCSN bifurcation for odd P ≥ and and b ∈ [0 , as τ P ( b ) = P ± b . (24)The solution is a special case. The stable (cid:2) Z, ¯ H, ¯ D (cid:3) − solution and the unstable (cid:2) Z, D, ¯ H (cid:3) − undergo a BCSNbifurcation along (cid:0) b, n + − b (cid:1) , b < , n ∈ Z + . Thestable (cid:2) Z, ¯ D, ¯ H (cid:3) − solution and the unstable (cid:2) Z, ¯ H, D (cid:3) − t x ( t ) b =0.50=1.17 ( a ) t x ( t ) b =2.10=1.17 ( b ) b x D ( c ) b = 0.5 ( d ) b x D ( e ) b = 2.1 ( f ) Figure 5. BCSN bifurcations at the boundary of the tongue. The stable solution is plotted in solid blue, while theunstable solution is plotted in dashed red. In ( a ) , the stable solution (cid:2) Z, ¯ D, D, ¯ H, ¯ D, D, ¯ D (cid:3) − and the unstable solution (cid:2) Z, D, ¯ D, D, ¯ H, ¯ D, D (cid:3) − are plotted close to (cid:0) b, − b (cid:1) .In ( b ) , the stable solution (cid:2) Z, ¯ D, D, ¯ H, ¯ D, D, ¯ D (cid:3) − and theunstable solution (cid:2) Z, ¯ D, ¯ H, H, D, ¯ H, ¯ D, ¯ Z, D, Z, ¯ D (cid:3) − areplotted close to (cid:16) b, b +2 b +54( b +1) (cid:17) . The D events associated withthe BCSN bifurcation for b < are plotted as triangles, whilethose associated with the BCSN bifurcation for b > areplotted as circles; ( c ) , ( d ) , ( e ) and ( f ) show their evolution as b and τ vary. undergo a BCSN bifurcation along (cid:0) b, n + b (cid:1) , b < , n ∈ Z + .A selection of these curves can be seen plotted in solidwhite in Fig. 3. There is a conspicuous outlier to thispattern. The solution does not undergo a BCSNbifurcation at the right boundary of the tongue for b > . We will now investigate why this is the case. B. Border-collision torus bifurcation of the solution Fig. 2 shows that the tongue has the same shapeas the other P : 1 tongues, rooted at τ = 0 . . Alongsome of the boundary, the solution changes to the solution without bifurcation. We test the observedsequences which represent the solution and find thatthe solution exists in the region [1 , × [0 . , outsideof the tongue. However, the solution is not alwaysstable; this can be seen in the top half of the regionin Fig. 2. In order to investigate this phenomenon, weconstruct a Poincaré map T by a similar method as forthe BCSN bifurcation. For τ > . , the solutionis represented by (cid:2) Z, ¯ D, D, ¯ H, ¯ D (cid:3) − from which we derive T + for x D ≥ . For x D < , the sequence changes to (cid:2) Z, ¯ H, H, ¯ D, ¯ Z, D, Z, ¯ H, ¯ D (cid:3) − , from which we derive T − .There is a significant difference however; the solutionis one-dimensional, but for τ > . , the solutionis two-dimensional. This occurs because the additional ¯ H and H events occur between the D that crosses x =0 and the Z present in both sequences. Therefore, wealso need to know the time t ¯ H at which the ¯ H presentin both sequences occurs. We derive T : ( x D , t ¯ H ) T → (cid:0) x ∗ D , t ∗ ¯ H − (cid:1) T as T (cid:18) x D t ¯ H (cid:19) = (cid:32) − − b +1 2 b +1 (cid:33) (cid:32) x D t ¯ H (cid:33) + C x D ≥ (cid:32) bb − − − b +1 2 b +1 (cid:33) (cid:32) x D t ¯ H (cid:33) + C x D < (25)where C = (cid:18) − b bb +1 + τ − (cid:19) . Note that T + has two rowswhich are multiples of each other; this is due to the solution being one-dimensional, and results in a eigen-value. By setting x D = 0 and solving T for τ , we ob-tain the curve τ = + b − b b +1) , which matches the upperright boundary of the tongue spanning from (1 , to (3 , . . T has a single fixed point which crosses x D = 0 at the boundary of the tongue. For x D > thefixed point is always stable. For x D < the fixed pointis stable for b > . , and loses stability when a pairof complex conjugate eigenvalues pass through | λ | = 1 ,resulting in a NS bifurcation at b = 2 . . | λ ± | < for b > . . Therefore, for τ > . the solution isstable for b > . . This agrees with the shape of theregion in which the solution is observed in Fig. 2.For b < . , we observe an interesting interactionbetween T + and T − . As T + has a eigenvalue, any point ( x D , t ¯ H ) for x D > is mapped directly to a nullcline on Figure 6. Border-collision Neimark-Sacker (BCNS) bifurca-tion of the fixed point of the solution. The upper twoplots show the BCNS bifurcation of the half-period map T .The lower plot is based on simulation using our iterative map,and shows the stable closed invariant curve expanding fromthe fixed point from x = 0 . The fixed point is plotted as ablue circle where it is stable, and as a red triangle where it isunstable. which the fixed point sits for x D > . The trajectorythen undergoes decaying oscillations to the fixed pointin the nullcline. However, if the fixed point is close to x D = 0 , T + can map the point into x D < , where thenullcline does not exist. In the absence of a stable fixedpoint at x D < , the trajectory starts to spiral out toinfinity. As this spiral must cross x D = 0 eventually, thetrajectory is caught by the nullcline again, causing an in-teresting half-spiral attraction. When the fixed point isat x D < , the trajectory converges to a stable attractorthat strikes the nullcline at multiple points. Fig. 6(a,b)show the interaction between T + and T − . Fig. 6(c) showsthe results of simulations in a ( b, τ ) sweep that crosses theboundary of the Arnold tongue. The simulated re-sults agree with those derived from T , and confirm thatthe closed invariant curve is born from the fixed pointas it crosses x D = 0 . We consider this to be a border-collision Neimark-Sacker (BCNS) bifurcation of T , simi-lar to that seen by Meiss . It corresponds to a border-collision torus (BCT) bifurcation of the solution inthe continuous system. The BCT bifurcation is super-critical, as a stable closed invariant curve expands fromthe fixed point when it becomes unstable. The BCT bi-furcation of the solution and the T bifurcation of the solution are plotted in black in Fig. 3.0 Figure 7. x-coordinates of H events from simulations. Thefour plots show the different cases of the solution at theborder. C. Border-collision torus bifurcation of the solution Now that we have studied the BCT bifurcation of the solution, we return to the solution. Previously wenoted that the Poincaré map P only proved the existenceof the torus bifurcation along the vertical boundaries ofthe locked region. We can now consider the horizontalboundaries in terms of border collisions. Recall that the solution is represented by four different sequences. At τ = n , n ∈ N , the sequence representing the solutionchanges when a H and a ¯ H cross x = 0 . This results ina change in the dimension of P . Let P n denote P when P contains an n × n matrix. When τ increases past τ = n , P changes from P n to P n+1 . Let us consider x H , theposition of the H event, as the fixed point of the map.Then x H = 0 when P changes from P n to P n+1 , resultingin a border collision at x H = 0 . Fig. 7 shows four casesthat occur in the border collision when we sweep acrossthe border τ = 1 . for fixed b .Fig. 7(a) shows the stable solution remaining sta-ble, as both P and P are stable at b = 6 . . Fig. 7(b)shows the stable solution becoming unstable, gen-erating a 13 : 13 solution in a supercritical BCT bifurca-tion. This solution exists in one of the vertical stripes wenoted in Fig. 2. Fig. 7(c) shows the stable solutionbecoming unstable, generating an aperiodic solution ina supercritical BCT bifurcation, showing that there aregaps between the vertical stripes. Fig. 7(d) shows thestable solution becoming unstable; however, the ape-riodic solution that the system converges to afterwards isnot generated at the border, indicating that the BCT bi- furcation is subcritical at this point. Therefore the BCTbifurcation of the solution must change criticality atsome point along τ = 1 . . The BCT bifurcation of the solution produces the horizontal black lines in Fig. 3,completing the boundary of the locked region.We examine the difference in maxima of nearby solu-tions on either side of τ = 1 . in Fig. 3(b). The point ( b, τ ) = (1 . , . stands out; the difference in maxima issudden to the left of that point, and gradual to the right.The gradual change in maxima occurs where the BCT bi-furcation is supercritical, and the abrupt change in max-ima occurs where the BCT bifurcation is subcritical, as il-lustrated by Fig. 7(b-d). We note that ( b, τ ) = (1 . , . forms one corner of a roughly triangular region contain-ing vertical P : P stripes seen in Fig. 2; every P : R solutionin this region is P : P . This leads us to the dashed whitecurve in Fig. 3. We refer to this curve as the D ¯ D curve.To the right of the D ¯ D curve, all D events occur at x < and all ¯ D events occur at x > , and so a legal sequencecannot contain a D adjacent to a ¯ D ; a P : R solution inthis region must, therefore, be P : P . To the left of the D ¯ D curve, a legal sequence representing a solution cancontain a D adjacent to a ¯ D ; we refer to such a solution asa D ¯ D solution. The point at which the BCT bifurcationchanges criticality occurs where the D ¯ D curve intersectsthe boundary of the locked region.As noted previously, when sweeping from right to left,the system converges to the solution until it becomesunstable at the T bifurcation. Now, we see that whensweeping from left to right, the system converges to D ¯ D solutions until they vanish at the D ¯ D curve; we notethat a and a solution undergo a BCSN bifur-cation at the D ¯ D curve. If the T bifurcation lies to theleft of the D ¯ D curve, then the solution remains sta-ble while the D ¯ D solutions exist, resulting in regions ofmultistability. We observe that this behaviour, togetherwith the change in criticality of the BCT bifurcation, isreminiscent of a Chenciner bifurcation, a co-dimension-2bifurcation that was observed in the smooth system byKeane and Krauskopf where it produced rich dynam-ics.If the T bifurcation lies to the right of the D ¯ D curve,this produces regions where vertical stripes can occur be-tween the BCT bifurcation of the solution and the D ¯ D curve. Note that the vertical stripes do not stretchall the way between the D ¯ D curve and the BCT bifur-cation. There is a second condition that a region mustsatisfy for the existence of vertical stripes, which is shownby the dashed black curve in Fig. 3. Vertical stripes onlyoccur where the order of D , ¯ D , H , ¯ H symbols in thestable solution is consistent with the solution. Thedifference between stable P : P solutions and the unstable solution they coexist with lies only in the order of ZH , ¯ Z and ¯ H symbols. This arises because these solutionsare generated when Z , H , ¯ Z and ¯ H symbols swapped or-der in the sequence representation of the solution at τ = n , n ∈ N .1 V. CONCLUSION We thoroughly analysed an elementary two-parametersystem which combines the effects of time-delayed feed-back and periodic forcing. In spite of its simplicity, itdemonstrates a complex structure of Arnold tongues withzero-width shrinking points and a high degree of multi-stability. Due to the system being piecewise-linear, weare able to solve the system analytically using an iterativemap. We investigate the existence and stability of solu-tions through the development of a symbolic representa-tion of solutions and the analysis of the subsequently de-veloped Poincaré and border-collision maps. This analy-sis reveals that the Arnold tongues are bounded by curvesof border-collision saddle-node bifurcations of periodicorbits. Additionally, we find curves of torus bifurcationsconnected to curves of border-collision torus bifurcations,and investigate changes in the criticality of these bifur-cations.Our analysis sheds new light onto previously ob-tained results in related smooth systems, particularlythe El Niño Southern Oscillation climate model stud-ied by Keane, Krauskopf, and Postlethwaite . Com-paring the numerically calculated bifurcation structurefound in that paper with the analytically calculated bi-furcation structure found here reveals that the prominentfeatures of the smooth system (1,2) are preserved in thenon-smooth limit of κ → ∞ . Indeed, the solutions inthe smooth system generally appear as “smoothed out”counterparts to the piecewise-linear solutions of the non-smooth system. There are some significant differences.The tongues in the smooth system are not connected byshrinking points, which is to be expected in the pres-ence of nonlinearity, according to Simpson and Meiss .Additionally, the smooth system contains period dou-bling bifurcations, which we do not observe in our sys-tem, likely due to the loss of nonlinearity in the limitof κ → ∞ . However, the analysis presented here doesprovide new insights into dynamics previously observednumerically .Both delayed feedback and periodic forcing are verycommon mathematical model ingredients and can befound in a variety of models used to study, for exam-ple, laser dynamics and chimera states . The currentwork reveals the phenomena which are a genuine conse-quence of this combination. Therefore, we expect thiswork to be of interest to a wide readership. VI. ACKNOWLEDGEMENTS This research was funded by the Irish Research Coun-cil and McAfee LLC through the Irish Research CouncilEmployment-Based Postgraduate Programme. We alsothank Dr. Sorcha Healy and the Applied Data Scienceteam at McAfee for their technical support. VII. APPENDIXA. Determining legality of a sequence Each event in a sequence has a time and a positionassociated with it. The times can be used to generate asystem of simultaneous equations: • Each D occurs at a time t D = n , n ∈ Z • Each ¯ D occurs at a time t ¯ D = n + , n ∈ Z • Each H occurs at a time t H = t ¯ Z + τ , where t ¯ Z isa time at which a ¯ Z occurs • Each ¯ H occurs at a time t ¯ H = t Z + τ , where t Z isa time at which a Z occurs • Each Z is connected to the previous ¯ Z by the equa-tion showing that the displacement of the trajec-tory between the two events is . • Each ¯ Z is connected to the previous Z by the equa-tion showing that the displacement of the trajec-tory between the two events is .As there is one equation for each symbol, we can solvethe system of equation to calculate the times. We thensolve for the positions, using the events to calculate theslope between the calculated times. Finally, we generatea set of inequalities for each event: • The time at which each event occurs must be con-sistent with the order of sequence. • Each event occurring between consecutive Z and ¯ Z events occurs at x > . • Each event occurring between consecutive ¯ Z and Z events occurs at x < .If the solved system of equations satisfies the set of in-equalities then the sequence is legal, for a given ( b, τ ) .These techniques allow us to use symbolic sequences tostudy solutions systematically. Solving for the times andpositions associated with events in a sequence, we canplot the solution represented by the sequence. Addition-ally, by changing the inequalities to equalities, we obtainthe bifurcation curves shown in Fig. 3. B. Calculation of the bifurcation curve of P For τ > , P can written as P ( S z ) = AS z + B, (26)2where B = − − − − ... − b +2 τ ( − n − b +( − n − − . (27)Due to the sparsity of A , the characteristic equation canbe readily calculated as ( − − λ )( − λ ) n − + ( − n − − n − b + ( − n − = 0 . (28)This simplifies to λ n + λ n − − − n − b + ( − n − = 0 . (29)By substituting λ = e iρ , we can solve for b to determine b = b bif for which the fixed point of P is bifurcating as e iρn + e iρ ( n − − − n − b bif + ( − n − = 0 . (30)Note that as − n − b bif +( − n − ∈ R , e iρn = e iρ ( n − . (31)So, e iρ (2 n − = 1 = e i πk , k ∈ Z , and so ρ = πk n − . Inorder to satisfy Eq. (31), k = n − . Substituting ρ backinto Eq. (30), cos (cid:18) πn ( n − n − (cid:19) +cos (cid:32) π ( n − n − (cid:33) − − n − b bif + ( − n − = 0 , (32)which simplifies under a sum-to-product cosine identityto cos (cid:18) π ( n − n − (cid:19) − b bif + ( − n − = 0 . (33)Therefore, b bif ( n ) = 1cos (cid:16) π ( n − n − (cid:17) − ( − n − , (34)where n = (cid:100) τ (cid:101) . REFERENCES S. Wieczorek, B. Krauskopf, T. B. Simpson, and D. Lenstra,“The dynamical complexity of optically injected semiconductorlasers,” Phys. Rep. , 1 (2005). S. P. Beeby, R. Torah, M. Tudor, P. Glynne-Jones, T. O’donnell,C. Saha, and S. Roy, “A micro electromagnetic generator forvibration energy harvesting,” J. Micromech. Microeng. , 1257(2007). S. Daneshgar, O. De Feo, and M. P. Kennedy, “Observationsconcerning the locking range in a complementary differential lcinjection-locked frequency divider–part I: Qualitative analysis,”IEEE Trans. Circuits Syst. I , 179 (2010). E. Tziperman, L. Stone, M. A. Cane, and H. Jarosh, “El Niñochaos: Overlapping of resonances between the seasonal cycle andthe Pacific ocean-atmosphere oscillator,” Science , 72 (1994). U. Parlitz and W. Lauterborn, “Superstructure in the bifurcationset of the Duffing equation,” Phys. Lett. A , 351 (1985). S. Boccaletti, A. N. Pisarchik, C. I. Del Genio, and A. Amann, Synchronization: from coupled systems to complex networks (Cambridge University Press, 2018). A. Marchionne, P. Ditlevsen, and S. Wieczorek, “Synchroni-sation vs. resonance: Isolated resonances in damped nonlinearoscillators,” Physica D , 8 (2018). T. Heil, I. Fischer, W. Elsässer, J. Mulet, and C. R. Mirasso,“Chaos synchronization and spontaneous symmetry-breaking insymmetrically delay-coupled semiconductor lasers,” Phys. Rev.Lett. , 795 (2001). L. Larger, B. Penkovsky, and Y. Maistrenko, “Virtual chimerastates for delayed-feedback systems,” Phys. Rev. Lett. ,054103 (2013). E. Schöll, G. Hiller, P. Hövel, and M. A. Dahlem, “Time-delayedfeedback in neurosystems,” Phil. Trans. Royal Soc. A , 1079(2009). J. Runge, V. Petoukhov, and J. Kurths, “Quantifying thestrength and delay of climatic interactions: The ambiguities ofcross correlation and a novel measure based on graphical mod-els,” J. Climate , 720 (2014). K. Pyragas, “Control of chaos via extended delay feedback,”Phys. Lett. A , 323 (1995). J. K. Hale and S. M. V. Lunel, Introduction to functional differ-ential equations (Springer Science & Business Media, 1993). H. Hu, E. H. Dowell, and L. N. Virgin, “Resonances of a harmon-ically forced Duffing oscillator with time delay state feedback,”Nonlin. Dynamics , 311 (1998). A. Maccari, “Vibration control for the primary resonance of thevan der Pol oscillator by a time delay state feedback,” Int. J.Non-Linear Mechanics , 123 (2003). M. Ghil, I. Zaliapin, and S. Thompson, “A delay differentialmodel of ENSO variability: parametric instability and the dis-tribution of extremes.” Nonlin. Processes Geophys. (2008). A. Keane, B. Krauskopf, and C. Postlethwaite, “Delayed Feed-back Versus Seasonal Forcing: Resonance Phenomena in an ElNiño Southern Oscillation Model,” SIAM J. Appl. Dyn. Syst. ,1229 (2015). A. Keane, B. Krauskopf, and C. Postlethwaite, “Investigatingirregular behavior in a model for the El Niño Southern Oscilla-tion with positive and negative delayed feedback,” SIAM J. Appl.Dyn. Syst. , 1656 (2016). R. D. Nussbaum, “Uniqueness and nonuniqueness for periodicsolutions of ˙ x ( t ) = − g ( x ( t − ,” J. Diff. Equ. , 25 (1979). S.-N. Chow and H.-O. Walther, “Characteristic multipliers andstability of symmetric periodic solutions of ˙ x ( t ) = − g ( x ( t − ,”Trans. American Math. Soc. , 127 (1988). Y. Wei-Ming and H. Bai-Lin, “How the Arnold tongues becomesausages in a piecewise linear circle map,” Comm. Theo. Phys. , 1 (1987). D. J. W. Simpson and J. D. Meiss, “Shrinking point bifurcationsof resonance tongues for piecewise-smooth, continuous maps,”Nonlinearity , 1123 (2009). D. J. W. Simpson, “The structure of mode-locking regions ofpiecewise-linear continuous maps: I. nearby mode-locking regionsand shrinking points,” Nonlinearity , 382 (2016). D. J. W. Simpson, “The structure of mode-locking regions ofpiecewise-linear continuous maps: II. skew sawtooth maps,” Non-linearity , 1905 (2018). N. Metropolis, M. Stein, and P. Stein, “On finite limit sets fortransformations on the unit interval,” J. Combinat. Theory, Se-ries A , 25 (1973). A. N. W. Hone, M. V. Irle, and G. W. Thurura, “On theNeimark-Sacker bifurcation in a discrete predator-prey system,”J. Biol. Dyn. , 594 (2010). M. di Bernardo, F. Garofalo, L. Iannelli, and F. Vasca, “Bifur-cations in piecewise-smooth feedback systems,” Int. J. of Control , 1243 (2002). S. Banerjee and C. Grebogi, “Border collision bifurcations intwo-dimensional piecewise smooth maps,” Phys. Rev. E , 4052(1999). A. Colombo, M. di Bernardo, S. J. Hogan, and M. R. Jeffrey,“Bifurcations of piecewise smooth flows: Perspectives, method-ologies and open problems,” Physica D , 1845 (2012). M. di Bernardo, C. Budd, A. Champneys, P. Kowalczyk,A. Nordmark, G. Tost, and P. Piiroinen, “Bifurcations in Nons-mooth Dynamical Systems,” SIAM Rev. , 629 (2008). A. Granados, M. Krupa, and F. Clément, “Border Collision Bi-furcations of Stroboscopic Maps in Periodically Driven SpikingModels,” SIAM J. Appl. Dyn. Syst. , 1387 (2014). D. Meiss, “Neimark-Sacker Bifurcations in Planar, Piecewise-Smooth, Continuous Maps,” SIAM J. Appl. Dyn. Syst. , 795(2008). H. E. Nusse, E. Ott, and J. A. Yorke, “Border-collision bifur-cations: An explanation for observed bifurcation phenomena,”Phys. Rev. E , 1073 (1994). U. Holmberg, Relay feedback of simple systems. , Ph.D. thesis,Lund Institute of Technology (1991). A. Colombo, M. di Bernardo, S. J. Hogan, and P. Kowalczyk,“Complex dynamics in a hysteretic relay feedback system withdelay,” Journal of Nonlinear Science , 85–108 (2007). A. Keane and B. Krauskopf, “Chenciner bubbles and torus break-up in a periodically forced delay differential equation,” Nonlin-earity , R165 (2018). D. J. W. Simpson and J. D. Meiss, “Resonance near border-collision bifurcations in piecewise-smooth, continuous maps,”Nonlinearity , 3091 (2010). T. Sorrentino, C. Quintero-Quiroz, A. Aragoneses, M. C. Torrent,and C. Masoller, “Effects of periodic forcing on the temporallycorrelated spikes of a semiconductor laser with feedback,” Opt.Expr. , 5571 (2015). V. Semenov, A. Zakharova, Y. Maistrenko, and E. Schöll,“Delayed-feedback chimera states: Forced multiclusters andstochastic resonance,” EPL115