## Abstract

A two-fold is a singular point on the discontinuity surface of a piecewise-smooth vector field, at which the vector field is tangent to the discontinuity surface on both sides. If an orbit passes through an invisible two-fold (also known as a Teixeira singularity) before settling to regular periodic motion, then the phase of that motion cannot be determined from initial conditions, and, in the presence of small noise, the asymptotic phase of a large number of sample solutions is highly random. In this paper, we show how the probability distribution of the asymptotic phase depends on the global nonlinear dynamics. We also show how the phase of a smooth oscillator can be randomized by applying a simple discontinuous control law that generates an invisible two-fold. We propose that such a control law can be used to desynchronize a collection of oscillators, and that this manner of phase randomization is fast compared with existing methods (which use fixed points as phase singularities), because there is no slowing of the dynamics near a two-fold.

## 1. Introduction

When an orbit of a system of ordinary differential equations approaches a stable periodic orbit, rather than the location of the orbit in phase space at a given time, a more useful quantity is often the asymptotic (or latent) phase of the orbit. The asymptotic phase distinguishes the long-term behaviour of different orbits, and sets of points with the same asymptotic phase are called *isochrons* [1]. Isochrons cannot intersect, but may accumulate at *phaseless sets* where the asymptotic phase is undefined. Near a point in such a set, termed a *phaseless point*, the asymptotic phase is highly sensitive to perturbations.

In smooth systems, phaseless points are unstable equilibria that can only be approached asymptotically (in backward time). In piecewise-smooth systems, we show that there exist phaseless points that can be intersected in finite time, without any slowing of the dynamics suffered near an equilibrium. The two-fold singularity, an enigma of piecewise-smooth dynamical systems theory, is one such phaseless point.

There are many reasons why it might be desirable to alter the phase of an existing oscillatory motion. One approach to tackling Parkinson’s disease, for instance, is to disturb the synchronized neural activity associated with physical tremors [2,3]. Desynchronization can be achieved with pulses [4], pulse trains [5], time-delayed control [6] or some other well-chosen feedback law [7]. Phase randomization for prototypical neuron models is achieved in [8] by briefly applying a control that transports the orbit to the close proximity of a phaseless point. The orbit subsequently returns to periodic motion but now has a different asymptotic phase. This phase is highly sensitive to the precise point at which the orbit is located when the control is removed, so small stochasticity in the system causes the resulting asymptotic phases of different neurons to be highly randomized.

Two-fold singularities of piecewise-smooth systems were first described by Filippov [9] and have garnered interest as points where both the forward and backward time uniqueness of a flow can break down in an otherwise deterministic flow [10–12]. In a vector field that is discontinuous on some hypersurface—the *discontinuity surface*—a two-fold is a point where the flow has quadratic tangencies (‘folds’) to both sides of the surface. Two-folds occur generically at isolated points in three-dimensional piecewise-smooth vector fields, and on (*n*−2)-dimensional manifolds in *n*-dimensional piecewise-smooth vector fields [12]. They have been identified in models of electronic circuits [13] and contact mechanics [14], and have deep connections to folded nodes of slow–fast systems [15]. The dynamics near a two-fold depends on whether the vector field on either side is curving towards or away from the discontinuity surface, and on the alignments of the two vector fields relative to each other [12].

In this paper, we illustrate the practical implications of a two-fold for phase randomization. As well as suggesting how two-folds might affect real systems, this suggests a use for them as control devices for the desynchronization of oscillators. To this end, we shall consider systems where a two-fold either occurs naturally in a piecewise-smooth system, or is introduced to a smooth system via a discontinuous control. We will simulate orbits in the presence of small noise as they pass close to a two-fold, and focus on how their phases are randomized by the presence of the singularity. The two-fold is not a zero of the vector field, so there is no slowing of the flow during phase randomization achieved in this way.

We shall first demonstrate that the two-fold singularity is a phaseless set, before deriving and simulating its effect on the distribution of the phase as a flow evolves through the singularity. In particular, we show that the distribution of the phase can be manipulated with higher-order terms in the vector field. We also explain how phase randomization can be achieved by using a discontinuous control action to generate a two-fold.

Our central mathematical result consists of a simple procedure by which the distribution of the asymptotic phase can be approximated. The results of this approximation are illustrated in figure 1 (to be produced in §4). To produce figure 1, we shall consider the normal form of the two-fold, add higher-order terms such that orbits emanating from the two-fold approach a stable periodic orbit, and add small noise. Figure 1*a* shows a histogram of the phase *ϕ*_{T}, computed relative to a reference time *T* and limiting to the asymptotic phase as *b* shows an analogous histogram using different higher-order terms, indicating that nonlinearity in the flow influences the distribution. The black curves show our theoretical approximations for the probability density functions of the phase derived in §4b.

Building on these observations, we propose that a periodic orbit in a smooth system can, via the temporary application of a discontinuous control, be made to take an excursion through a two-fold and then return to periodic motion, but now with a randomized phase. We provide simulations of a simple model in which the effect is to desynchronize a collection of synchronized oscillators. As discussed in §6, two-folds have potential advantages over equilibria as phaseless points in the context of phase randomization. Although we provide a motivating toy example, calculating and optimizing realizable control actions (as in [8]) are beyond the scope of this paper.

In §2, we briefly set out the normal form equations of the invisible two-fold singularity and review a few pertinent details of the local dynamics. Here we include a novel definition of polar coordinates that fit naturally with the local dynamics (and should be of particular interest to dynamicists studying two-folds). We include higher-order terms with the normal form in §3, and extend these definitions before simulating a flow through the two-fold subjected to stochastic perturbations in §4. In §5, we give an example of an application to desynchronize smooth oscillators using a discontinuous control, and provide closing remarks in §6.

## 2. Deterministic dynamics of the normal form

There exist three generic kinds of two-fold singularity, depending on whether both folds are visible (curving away from the discontinuity surface) or invisible (curving towards the discontinuity surface), or one is visible and one is invisible. In each case, there exist local conditions under which the flow traverses the singularity in finite time. In this paper, we focus on the case where both folds are invisible, called an invisible two-fold or Teixeira singularity [10,16], in which the local flow passes through the singularity in finite time. The added interest of the invisible two-fold is that the local flow winds repeatedly around the singularity, giving oscillatory behaviour.

In three dimensions, *X*=(*x*,*y*,*z*), the normal form of the invisible two-fold is the piecewise-smooth system
*X*=(0,0,0). The discontinuity surface *x*=0 consists of attracting and repelling sliding regions, denoted *A* and *R* (where the flow is confined to the surface and so slides along it), and two crossing regions *C*^{±} (where the flow passes transversally though the surface), as illustrated in figure 2.

The system is solved across the discontinuity by forming the Filippov system for (2.1). Using *e*_{i} to denote the *i*th coordinate vector in *Filippov solution* of (2.1) is an absolutely continuous function *φ*(*t*) that satisfies
*t* [9]. We let *φ*(*t*;*X*_{0},*t*_{0}) denote a Filippov solution to (2.1) with initial condition *φ*(*t*_{0})=*X*_{0}. Each *φ*(*t*;*X*_{0},*t*_{0}) is a concatenation of several smooth orbit segments, including segments of evolution outside the discontinuity surface *x*=0, and segments of evolution on *x*=0 referred to as *sliding* motion. Solving the convex combination in (2.3) for *A* and *R* as
*A*, or slide out of the singularity in *R*, they do so tangent to the weak eigenvector of the matrix (‘weak’ meaning associated with the smallest eigenvalue) and so are tangent to the line *z*=(λ−*V*^{−})*y* on *x*=0. The time to reach or depart the singularity is finite. Specifically, the time is *t*=(*V*^{−}−λ−1)*y*/λ from an initial point on the line *z*=(λ−*V*^{−})*y* with *x*=0. Throughout the paper, we use the values

In the remainder of this section, we review crossing solutions to (2.1) in §2a, introduce polar coordinates for describing the flow as it spirals away from the two-folds in §2b, and define the notion of ‘viable’ Filippov solutions in the context of simulation in §2c. In later sections, we apply these concepts to more general piecewise-smooth systems.

### (a) Crossing dynamics and an unstable cone

The left and right half systems of (2.1), *X*_{0}=(*x*_{0},*y*_{0},*z*_{0}). In what follows, we use (2.7) and (2.8) to study Filippov solutions *φ*(*t*;0,*y*,*z*,*t*_{0}) of (2.1). As indicated in figure 2, orbits wind around regions *A* and *R*, repeatedly crossing *x*=0 in regions *C*^{±}. This can be understood with a return map that has been studied in considerable detail [16,18]. Here, we review a few details of this map that are needed for the ensuing analysis.

First consider the region *z*<0 on *x*=0, where the flow *φ*_{L} is directed away from the discontinuity surface. If *y*>0, then (0,*y*,*z*)∈*C*^{−}, so *φ*_{R} is directed towards the surface and hence *φ*(*t*;0,*y*,*z*,*t*_{0}) immediately enters *x*<0. If instead *y*<0, then (0,*y*,*z*)∈*R*, so both *φ*_{L} and *φ*_{R} are directed away from the surface and hence *φ*(*t*;0,*y*,*z*,*t*_{0}) is not uniquely determined (it may slide along *x*=0 following (2.4) or enter either *x*<0 or *x*>0 at any instant). For definiteness, when *y*≤0, we consider the Filippov solution *φ*(*t*;0,*y*,*z*,*t*_{0}) that immediately enters *x*<0. From (2.7), we find that *φ*(*t*;0,*y*,*z*,*t*_{0}) resides in *x*<0 until the later time *t*=*t*_{0}−2*z* when it is located at (0,*y*−2*V*^{−}*z*,−*z*).

Second, consider *y*<0 and assume *φ*(*t*;0,*y*,*z*,*t*_{0}) immediately enters *x*>0. From (2.8), and by applying arguments analogous to those in the previous paragraph, we find that *φ*(*t*;0,*y*,*z*,*t*_{0}) resides in *x*>0 until *t*=*t*_{0}−2*y* when it is located at (0,−*y*,−2*V*^{+}*y*+*z*).

From these observations, we can construct a return map capturing crossing dynamics. This is formulated in the following proposition (refer to [16] for a complete proof).

### Proposition 2.1

*Suppose V*^{−}, *V*^{+}<0, *V*^{−}*V*^{+}>1. *For any* *and z*<0, *let* (0,*y*′,*z*′) *denote the next intersection of φ*(*t*;0,*y*,*z*,*t*_{0}) *with the discontinuity surface at a point with y*′>0 (*assuming φ*(*t*) *immediately enters x*<0 *in the case y*<0). *If z*<2*V*^{+}*y*/(4*V*^{−}*V*^{+}−1), *then z*′<0 *and*
*and the evolution time of the orbit from* (*y*,*z*) *to* (*y*′,*z*′) *is*
*whereas if z*>2*V*^{+}*y*/(4*V*^{−}*V*^{+}−1), *then z*′>0.

The return map (2.9) is linear, area-preserving, and has a saddle-type fixed point at (*y*,*z*)=(0,0). The unstable multiplier (i.e. the eigenvalue of the matrix in (2.9) with modulus greater than 1) is
*γ*), where

By evolving the flow forwards from *ζ* according to (2.7) and (2.8), we obtain a surface *Λ* that is given implicitly by
*Λ*. The surface *Λ* is non-differentiable at *x*=0 and consists of part of a hyperbolic paraboloid on each side of *x*=0.

For any initial point on the ray *ζ*, the next intersection of *φ*(*t*;0,*y*,*γy*,*t*_{0}) with *C*^{−} is (0,*μy*,*μγy*), shown in figure 2. By substituting *z*=*γy* into (2.10), we find that this iteration of (2.9) corresponds to an evolution time of
*ζ* in finite time. Such an orbit is shown in figure 3.

Given any *a*>0, let *ψ*_{a}(*t*) (defined for *t*≥0) denote the unique Filippov solution that is located at the two-fold at *t*=0, and located on *ζ* at *t*=*a*. That is, *ψ*_{a}(0)=(0,0,0) and *ψ*_{a}(*a*)∈*ζ*. By using (2.16) and the classical formula for the sum of a geometric series, we determine the *y*-value of *ψ*_{a}(*a*) to be

After *t*=*a*, the Filippov solution *ψ*_{a}(*t*) next intersects *ζ* at *t*=*μa*. Consequently, *ψ*_{a}(*μ*^{n}*a*)∈*ζ*, for all *ψ*_{μna}(*t*) is the same orbit as *ψ*_{a}(*t*), and so we can characterize all the *ψ*_{a}(*t*) by restricting our attention to values of *a* in a fundamental domain *a**≤*a*<*μa**, for some *a**>0. Since *Λ* is solely composed of Filippov solutions *ψ*_{a}(*t*), we can write
*Λ* with *C*^{−} is *ζ*. Since *ζ* corresponds to an unstable eigenvector of the crossing map, we refer to *Λ* as an unstable manifold of the two-fold.

### (b) Polar coordinates

Since orbits of (2.1) on *Λ* spiral out from the two-fold (figure 3), it is natural to introduce polar coordinates. The surface *Λ* projects uniquely onto the (*x*,*y*)-plane (and the (*x*,*z*)-plane, but we use the former). The usual polar coordinates (*r*(*x*,*y*) and *θ*(*x*,*y*) (which we interpret as ‘polar coordinates’) in a way that the restriction of (2.1) to *Λ* takes a particularly simple form.

To obtain a suitable notion of the phase *θ*, we begin by finding the times *τ*_{L}>0 and *τ*_{R}<0 that an orbit takes to travel from *ζ* to a given point in *x*<0 and in *x*>0, respectively.

First, choose any *x*<0 and *y*_{0}>0 be such that by evolving (2.1) forwards from (0,*y*_{0},*γy*_{0}) we arrive at (*x*,*y*,*z*)∈*Λ*, for some *x*=0, and let *τ*_{L}>0 be the corresponding evolution time. This is illustrated in figure 3. Note that the value of *z* is given from (2.14) by *x*=−*Ξ*(*y*,*z*)/*V*^{−}.

We have *φ*_{L}(*τ*_{L};0,*y*_{0},*γy*_{0},0)=(*x*,*y*,*z*), and so, by (2.7),
*y*_{0} and *τ*_{L}. Eliminating *y*_{0} yields *y*_{0} is given by *y*_{0}=*y*−*V*^{−}*τ*_{L} with (2.21).

Second, choose any *x*>0 and *y*_{0}>0 be such that by evolving (2.1) backwards from (0,*y*_{0},*γy*_{0}) we arrive at (*x*,*y*,*z*)∈*Λ*, for some *x*=0, and let *τ*_{R}<0 be the corresponding evolution time. Here *z* is specified by *x*=*Ξ*(*y*,*z*)/*V*^{+}, and, in the same manner as above, we obtain
*y*_{0}=*y*−*τ*_{R}.

Using these expressions, we now build a useful set of polar coordinates for the projection of *Λ* onto the (*x*,*y*)-plane. We let the positive *y*-axis correspond to *θ*=0 and *r*=*y*, i.e.
*y*-axis corresponds to points on *ζ*. Naturally, we want to define *θ* such that a change of 2*π* corresponds to one complete revolution from *ζ* back to itself. By (2.16), forward evolution from (0,*y*,*γy*) returns to *ζ* in a time proportional to *y*. This suggests that we want *y* (or the distance from the two-fold, *r*). Moreover, the orbit returns to *ζ* at the point (0,*μy*,*μγy*). That is, the value of *r* increases from *y* to *μy* and so changes by an amount proportional to *y*, suggesting *r* and *θ* so that they satisfy (2.23) and
*α*,*β*>0. The following choice of *r* and *θ* achieves this.

### Proposition 2.2

*Suppose V*^{−}, *V*^{+}<0 *and V*^{−}*V*^{+}>1. *Let*
*and*
*where τ*_{L} *and τ*_{R} *are given by* (* 2.21*)

*and*(

*2.22*),

*and*

*where μ is given by*(

*2.11*).

*Then*(

*2.25*)

*and*(

*2.26*)

*define a continuous bijection from*

*to*

*that satisfies*(

*2.23*).

*Moreover, under this coordinate change the restriction of*(

*2.1*)

*to Λ is given by*(

*2.24*)

*with*(

*2.27*).

A constructive proof of proposition 2.2 using (2.7) and (2.8) is given in appendix A. Figure 4 illustrates contours of *r*(*x*,*y*) and *θ*(*x*,*y*) as defined in (2.25) and (2.26).

Given (2.26), the negative *y*-axis corresponds to *α*→0 as *V*^{−}*V*^{+}→1 and *V*^{−}/*V*^{+} is fixed as the limit is taken). Graphically, we have found also that 0<*β*<*π*.

The right-hand side of (2.24) is defined everywhere except at *r*=0 (the two-fold) and can be solved explicitly. Taking *r*→0 as *t*→*t*_{0} as an initial condition, for all *t*>*t*_{0} the solution to (2.24) is given by
*C*∈[0,2*π*). Therefore, the Filippov solution *ψ*_{a}(*t*), defined in §2a, may be written in polar coordinates as

### (c) Viable Filippov solutions

The following result shows that orbits emanating from the two-fold either dwell in the repelling sliding region *R* for some time, or evolve on *Λ* for all later times.

### Proposition 2.3

*Let φ*(*t*)=*φ*(*t*;0,0,0,0) *be a Filippov solution of* (*2.1*) *located at the two-fold at t*=0. *Then either*: (*i*) *φ*(*t*)∈*Λ* *for all t*>0 *or* (*ii*) *φ*(*t*)∈*R* *for some t*>0 *and φ*(*t*)∉*Λ for all t*>0.

### Proof.

First, suppose *φ*(*t*)∈*R* for some *t*>0. Since Filippov solutions to (2.1) can only enter the repelling region *R* by passing through the two-fold (illustrated in figure 2), then either *φ*(*t*)∈*R* for all *t*>0, in which case clearly *φ*(*t*)∉*Λ* for all *t*>0, or *φ*(*t*)∈*R* on time interval (0,*b*], and at *t*=*b*, the Filippov solution *φ*(*t*) is ejected into either *x*<0 or *x*>0. Suppose without loss of generality that *φ*(*t*) enters *x*<0. Then, subsequent crossing motion (for all *t*≥0, it can be shown [12]) is described by the return map (2.9) starting from *φ*(*b*). Since *φ*(*b*)∉*ζ*, iterations under (2.9) do not lie on *ζ*, hence *φ*(*t*)∉*Λ* for all *t*>0.

Second, suppose *φ*(*t*)∉*R* for all *t*>0. By solving backwards in time, the only way that solutions can reach the two-fold without intersecting *R* is through intersections with *ζ* (illustrated in figure 3). This is because orbits cannot reach the two-fold by evolving backwards in time on the attracting sliding region. Nor can solutions reach the two-fold by evolving purely in either the left half-space (*x*≤0) or the right half-space (*x*≥0) because such evolution follows (2.7) and (2.8). The only remaining possibility is to reach the two-fold via crossing motion. In view of the saddle-type nature of (*y*,*z*)=(0,0) for (2.9), this can only be achieved along the unstable eigenvector. Therefore *φ*(*t*)∈*Λ* for all *t*>0. ▪

Now, we introduce the notion of a ‘viable’ Filippov solution to represent an orbit that is robust for the purpose of forward time numerical simulation. Any orbit that evolves from a point in *R*, simulated using discretization or by modelling the switch using hysteresis, time-delay or noise, is likely to be immediately ejected into either *x*<0 or *x*>0. In view of proposition 2.3, any simulated orbit that arrives at the two-fold is likely to subsequently evolve on *Λ* (or at least very near *Λ*), as observed in [16]. Indeed, perturbations of a planar two-fold by hysteresis, time-delay and noise studied in [19] confirm that perturbed systems follow Filippov solutions provided they do not lie on a repelling sliding region.

### Definition 2.1

A Filippov solution is said to *viable* if it does not intersect (or travel along) a repelling sliding region in forward time.

For the system (2.1), *φ*(*t*;*X*_{0},*t*_{0}) is viable if and only if, for all *t*>*t*_{0}, *φ*(*t*;*X*_{0},*t*_{0})∉*R*. The idea is that viable Filippov solutions are the most relevant in real systems or in simulations. The last result of this section follows immediately from proposition 2.3 and provides us with a complete characterization of viable Filippov solutions to (2.1) that pass through the two-fold.

### Corollary 2.1

*Let φ*(*t*)=*φ*(*t*;0,0,0,0) *be a Filippov solution of* (*2.1*) *that is located at the two-fold at t*=0. *Then φ*(*t*) *is viable if and only if φ*(*t*)=*ψ*_{a}(*t*), *for some a.*

## 3. Deterministic global dynamics

The vector field near a generic invisible two-fold in a three-dimensional system can be transformed to
*F*_{L} and *F*_{R} are the constituents of the normal form (2.2), and

As described in [16], there exists a surface *Λ* in the limit |*X*|→0) on which viable Filippov solutions to (3.1) evolve. This surface intersects *x*=0 with *y*>0 on a curve *ζ* to first order. We let

Next, in §3a, we define the asymptotic phase for orbits of (3.1), study the associated isochrons in §3b, and find times of return to the discontinuity surface in §3c.

### (a) The asymptotic phase of an orbit

Let us suppose that (3.1) has a stable periodic orbit *Γ* of period *τ* such that orbits emanating from the two-fold on *Γ*. Examples are given below and also in [16]. In the examples considered below, *Γ* has exactly two intersections with the discontinuity surface, one with *y*>0 and one with *y*<0. We use the intersection point with *y*>0 as a reference point with zero phase.

Let *Γ* as *phase* of *T* that is assumed to be sufficiently large that the orbit lies close to *Γ*, is defined as
*s*_{T}≤*T* is the previous time at which the orbit lies on *x*=0 with *y*>0. The ‘*mod* 2*π*’ ensures *ϕ*_{T}=[0,2*π*), but is almost redundant, because the orbit is close to *Γ* and so *T*−*s*_{T} is unlikely to be greater than *τ*. Assuming the forward orbit of *X*_{0} is unique and converges to *Γ*, the *asymptotic phase* of *X*_{0} is defined as

### (b) Isochrons

An *isochron* is a set of points with the same asymptotic phase [1]. Isochrons were introduced by Winfree in [20] and shown to be (*n*−1)-dimensional *C*^{k} manifolds for *n*-dimensional *C*^{k} vector fields by Guckenheimer in [21]. Isochrons of the three-dimensional piecewise-smooth system (3.1) are therefore two-dimensional piecewise-smooth manifolds. Since we are only concerned with orbits that spiral out from the two-fold, it suffices to consider the intersection of the isochrons of (3.1) with *Γ* forms the boundary of this surface. The isochrons were computed by interpolating between computed points on the forward orbits and correspond to *ϕ*=2*πk*/5 for *k*=0,…,4. This method is described by Winfree in [1]. More sophisticated methods for computing isochrons are discussed in [22,23].

The isochrons accumulate at the two-fold where the asymptotic phase is undefined. The two-fold is a phaseless point and the basin of attraction of the two-fold, which is a three-dimensional region, is a phaseless set.

### (c) Times of return to the discontinuity surface

Here, we define a return time function *f* for the curve *a*>0, let *t*=*f*(*a*)>*a* be the next time at which *f* using (2.6) and (3.5) and was computed by numerical simulation.

For small *a*, *f*(*a*)≈*μa*, because near the two-fold the higher-order terms *G*_{L} and *G*_{R} have little effect and the return time is *μa* for the normal form (2.1), as stated in §2a. For large *a*, *Γ* and so *f*(*a*)≈*a*+*τ*.

In §4, we use *f* to extrapolate the probability density function for *ϕ*_{T} from points near the two-fold to points near *Γ*.

## 4. Stochastic dynamics and phase randomization

Here, we study a stochastic perturbation of (3.1) defined by the three-dimensional stochastic differential equation
*W*(*t*) is a standard three-dimensional Brownian motion, *D* is a 3×3 matrix of constants, and 0≤*ε*≪1 represents the noise amplitude. Given a sample solution to (4.1), *T*, we define *s*_{T} and *ϕ*_{T} by (3.4) in the same manner as for the deterministic system (3.1).

Here we show that the asymptotic phase is highly randomized for sample solutions to (4.1) that pass close to the two-fold before approaching a stable periodic orbit. We provide numerical evidence for this in §4a, then derive an approximation to the phase distribution in §4b. To illustrate the results, we use (2.6) and two choices for *G*_{L} and *G*_{R}, namely (3.5) and
*V*^{−} and *V*^{+} and different choices for *G*_{L} and *G*_{R}, we have observed similar results.

### (a) Sample solutions and phase randomization

A sample solution to (4.1) using (2.6) and (3.5) is shown in figure 7*a*. The initial point is *X*=(0,1,1), and we denote this solution as *t*_{0}=3.0445 (to four decimal places). The perturbed solution *t*≈*t*_{0}, the perturbed solution *ε* noise has an order-1 influence on the location and phase of *G*_{L} and *G*_{R} have little influence. At later times, *G*_{L} and *G*_{R} are important and *Γ*.

We computed 10^{4} sample solutions analogous to the solution shown in figure 7*a* (i.e. with the same initial point and parameter values). Figure 1*a* shows a histogram of the phase *ϕ*_{T} of these solutions. Here we used *T*=15, because transient dynamics appears to have decayed by this time. We note that the phase is roughly uniformly distributed on [0,2*π*).

Figure 7*b* shows a sample solution using (4.2) instead of (3.5). As in figure 7*a*, the solution approaches the two-fold (the deterministic sliding time is *t*_{0}=2.9763, to four decimal places), then spirals out towards a stable periodic orbit *Γ*. Figure 1*b* shows a histogram of the phase *ϕ*_{T} of 10^{4} sample solutions, using *T*=40. Again, the phase is highly random, but, in this case, the phase distribution is not well approximated by a uniform distribution. A phase of *ϕ*_{T}≈3*π*/2 appears to be about twice as likely as a phase of *ϕ*_{T}≈*π*/2. In further simulations (not shown) with different values of *ε* and *T*, we observed similar phase distributions to those shown in figure 1.

In §4a, we combine the polar coordinate description of the local dynamics (given in §2b) with the global return time function *f* (given in §3c) in order to construct approximations to the phase distributions of figure 1.

### (b) The probability distribution of the asymptotic phase

Let *t*_{0}>0. For sample solutions *ϕ*_{T} at large values of *T* (in order to adequately approximate the asymptotic phase *ϕ*) but our analysis requires considering all *T*>*t*_{0}.

For any *T*>*t*_{0}, let *p*_{T}(*s*_{T}) denote the probability density function for the value of *s*_{T} (the previous time of intersection with the discontinuity surface at a point with *y*>0). We begin by explaining that for intermediate values of *T*−*t*_{0}, specifically *ε*≪*T*−*t*_{0}≪1, it is suitable to assume that *s*_{T}−*t*_{0} has a reciprocal probability distribution, that is,

Since the noise amplitude *ε* is small, if *T*−*t*_{0}≪1 then an arbitrary sample solution *s*_{T}, we can ignore the higher-order terms *G*_{L} and *G*_{R}. Also if *ε*≪*T*−*t*_{0}, then with high probability an arbitrary sample solution will lie sufficiently far from the two-fold that for the purposes of computing *s*_{T} we can ignore the noise. Therefore, here it is suitable to restrict our attention to the normal form (2.1), and we work with this system in polar coordinates (2.24).

Any solution to (2.24) that limits to the two-fold as *t*→*t*_{0} (with *t*>*t*_{0}) is given by (2.28) for some *C*∈[0,2*π*). The system (2.24) is independent of *θ*, which implies we should treat *C* as a random variable with a uniform distribution on [0,2*π*).

Given *C*, we can calculate *s*_{T} by using (2.28) to solve *θ*(*t*)=0 for *s*_{T}. We write *C* is uniformly distributed implies that *s*_{T}−*t*_{0} has a reciprocal probability distribution as given by (4.3).

As sample solutions *s*_{T} by a small amount. Therefore, in order to approximate *p*_{T} for larger values of *T*>*t*_{0}, we can again ignore the noise (as long as *T* is not so large that the contribution of the accumulated small diffusive effects is significant). In this approximation, *p*_{T} can be expressed in terms of *p*_{f−1(T)} by iterating the density under *f*, specifically
*n* times gives

We can therefore approximate *p*_{T} for a large value of *T* by using (4.7) with a value of *n* that is sufficiently large that *f*^{−n}(*T*)−*t*_{0} is small, and so the reciprocal probability density function (4.3) can be used for *p*_{f−n(T)}. This is the manner by which we obtained the two approximations shown in figure 1, using *n*=10 in both cases. The iterative procedure (4.7) was performed numerically, because analytic expressions for *f* appear to be unavailable.

Figure 8 shows plots of *f*^{n} for both (3.5) and (4.2). We note that if *f*^{n}(*a*) is an affine function of *p*_{T}. Indeed, in figure 8*a*, *f*^{n} is indistinguishable from an affine function on the given scale and so, in figure 1*a*, the distribution is approximately uniform. In figure 8*b*, *f*^{n} has a notable nonlinearity, and, for this reason, the distribution in figure 1*b* is significantly non-uniform.

## 5. Desynchronizing a collection of smooth oscillators

Now suppose that we wish to break the synchrony of a collection of coordinated oscillators. In [8], this is achieved by applying a control that pushes the state of the oscillators towards an unstable equilibrium, where small random perturbations efficiently randomize the phase of the oscillators. Upon removal of the control action, each oscillator returns to the original regular periodic motion but now the oscillators are desynchronized. Here we propose a ‘fast’ method to achieve this, using a two-fold as the phase singularity rather than an equilibrium; in that case there is no slowing of the dynamics as the oscillators return to periodic motion. The example we give is intended to illustrate this concept in a simple way, with practical applications left to future work.

To illustrate our method, we assume the oscillators are governed by the planar system

Consider the two-dimensional stochastic differential equation
*W*(*t*) is a standard two-dimensional Brownian motion, and 0≤*ε*<1 is the noise amplitude. The function *H* is the Heaviside function and *t*_{1} and *t*_{2} are the start and end times of the control action *c*(*t*). In order to generate a generic two-fold, *c*(*t*) must be discontinuous and time-dependent. Here we show that the piecewise-linear form

By treating the time *t* as a third dynamic variable, i.e. with *t*_{1}<*t*<*t*_{2} with the discontinuity surface *x*=0. The fold lines on *x*=0 (where *y*=*a*_{1}*t* and *y*=*a*_{3}*t*, so the origin (*x*,*y*,*t*)=(0,0,0) is a two-fold. Via elementary calculations, we find that we need *a*_{2}<*a*_{1}, *a*_{3}<*a*_{4} and *a*_{1}≠*a*_{3} in order for the two-fold to be generic and for both folds to be invisible, as in (2.1). We also require *a*_{1}<*a*_{3} in order to have *V*^{−},*V*^{+}<0 and *V*^{−}*V*^{+}>1, as in (2.5). In summary, we require

After some experimentation, we found the values
*a*_{i}, larger negative values of *t*_{1} are required in order for the control to direct orbits into the two-fold; in that case the phase randomization is not achieved as quickly. With larger values of the *a*_{i}, we observed distributions of the asymptotic phase that were substantially non-uniform; in that case the control does not ‘fully randomize’ the phase.

We require *t*_{1}<0 and |*t*_{1}| to be sufficiently large that all points on the periodic orbit of (5.1) are directed into the two-fold from *t*=*t*_{1} to *t*=0. The size of the set of all points (*x*,*y*) that are directed into the two-fold increases with |*t*_{1}|, and with (5.5) this set contains the periodic orbit for all *t*_{1}≤−4.1, approximately. For this reason, we use *t*_{1}=−5. The evolution of the periodic orbit from *t*=*t*_{1} to *t*=0 is shown in figure 9.

Finally, we would like the value of *t*_{2}>0, at which time the control is turned off, to be small so that the phase randomization is achieved relatively quickly. However, with smaller values of *t*_{2} solutions lie nearer to (*x*,*y*)=(0,0) at *t*=*t*_{2}. Since this point is an equilibrium of (5.1), smaller values of *t*_{2} cause solutions to dwell near (0,0) for longer periods of time, and so the phase randomization is not affected as quickly. Here we use *t*_{2}=2.5 as it proves to be a suitably intermediate value. In figure 9, the unstable manifold *t*=0 to *t*=*t*_{2}.

Figure 10 illustrates the phase randomization using (5.5), *t*_{1}=−5 and *t*_{2}=2.5. We see that the randomization appears to be highly effective. A quantitative comparison of our phase randomization with other methods remains for future work.

## 6. Discussion

In this paper, we have studied orbits that pass through an invisible two-fold (Teixeira singularity) then limit to a stable periodic orbit. We considered the presence of small noise in order to remove the ambiguity of evolution through the two-fold. We found that the phase of the limiting periodic motion is highly randomized and that the probability distribution of the phase depends crucially on the nature of the nonlinear dynamics experienced by orbits between the two-fold and the periodic orbit. By using polar coordinates to describe the motion of orbits as they initially depart from the two-fold, and a one-dimensional return map to describe the global dynamics approaching the periodic orbit, we constructed an approximation to the probability density function of the asymptotic phase that matches well to the results of numerical simulations (figure 1).

In §5, we showed how a simple discontinuous control law can generate a two-fold in a smooth system. We propose that in this fashion the two-fold can be used to desynchronize a collection of oscillators and that this may have the following advantages to desynchronization using an unstable equilibrium.

(i) With an unstable equilibrium, the control action must direct the state of each oscillator close to the equilibrium, and the effectiveness of the randomization correlates with the accuracy of this action. With a two-fold, however, the control does not need to be as precise, because it only needs to direct orbits to the basin of attraction of the two-fold.

(ii) With an unstable equilibrium, the effect of the randomization is due to inherent stochasticity in the system and any artificial randomness in the control law. With a two-fold, the randomization is inherent in the ambiguity of forward evolution through the two-fold.

(iii) Orbits dwell near unstable equilibria and so can take a relatively long time to return to regular periodic motion. With a two-fold there is no slowing of the dynamics, hence we refer to our desynchronization as ‘fast’ phase randomization. Indeed, in figure 10, we see that solutions return quickly to approximately regular periodic motion after the control is turned off. Whether this truly constitutes fast desynchronization requires a study of the full time for which the control action is required. This and a study of practical discontinuous control actions remain for future work.

Teixeira’s paper [10] inspired a legacy of intrigue around a specific case of the two-fold singularity, first studied more generally in [9]. Interest has only increased as its role as a determinacy-breaking singularity has become more clear [11,16]. This paper begins to reveal the practical side of these insights, including how two-folds may manifest in physical systems, and how they might be put to use as tools for control. Little is still known about where two-folds appear naturally in physical systems. Similarly, little is known about what applications they may have in control systems, but a role as a phase randomizer is suggested here. Perhaps the obvious next step is to design implementable control circuits to investigate the practical obstacles and opportunities that they present.

## Data accessibility

The figures in this paper are reproducible as we have provided all the necessary equations and parameter values in the text and figure captions.

## Authors' contributions

D.J.W.S. was responsible for calculations regarding polar coordinates and stochastic dynamics and performed the numerical simulations. M.R.J. was responsible for calculations regarding the control example. Both authors contributed roughly equally to the ideas and the writing of the paper.

## Competing interests

We have no competing interests.

## Funding

M.R.J.’s research is supported by EPSRC Fellowship grant no. EP/J001317/2.

## Appendix A. Proof of proposition 2.2

### Proof.

Since *ψ*_{a}(*a*)∈*ζ*, the *y* and *r*-values of *ψ*_{a}(*a*) are the same; see (2.23). The *y*-value of *ψ*_{a}(*a*) is given by (2.17). By (2.29), the *r*-value of *ψ*_{a}(*a*) is equal to *αa*. By matching these two expressions, we conclude that *α* is given by (2.27).

At the beginning of §2b, we showed that, given *x*<0 and *τ*_{L}>0 is given by (2.21) and *y*_{0}=*y*−*V*^{−}*τ*_{L}, then evolving (2.1) forwards from (0,*y*_{0},*γy*_{0}) for a time *τ*_{L} takes us to the point (*x*,*y*,*z*), where *x*=−*Ξ*(*y*,*z*)/*V*^{−} (2.15). Moreover, this defines a bijection from the region *x*<0, *y*_{0}>0, 0<*τ*_{L}<−2*γy*_{0}.

By (2.24), evolving (*r*,*θ*)=(*y*_{0},0) for a time *τ*_{L} takes us to the point
*y*_{0}>0, 0<*τ*_{L}<−2*γy*_{0} to *r*>0, *β*>0). By substituting (2.21) and *y*_{0}=*y*−*V*^{−}*τ*_{L} into (6.1) and simplifying, we arrive at (2.25) and (2.26) for *x*<0.

Similarly, given *x*>0 and *τ*_{R}<0 be given by (2.22) and *y*_{0}=*y*−*τ*_{R}. This is a bijection to *y*_{0}>0, −2*y*_{0}<*τ*_{R}<0, and evolving (2.1) backwards from (0,*y*_{0},*γy*_{0}) for a time *τ*_{R} takes us to (*x*,*y*,*z*), where *x*=*Ξ*(*y*,*z*)/*V*^{+} (2.15).

In polar coordinates evolving (*r*,*θ*)=(*y*_{0},0) for a time *τ*_{R} takes us to
*r*>0, *y*_{0}=*y*−*τ*_{R} into (6.2) leads to (2.25) and (2.26) for *x*>0.

Finally, by matching the limiting left and right values for *θ* on the negative *y*-axis as given by (2.26), we obtain
*μ*=(1−2*αγ*)/(1−2*α*) and solving (6.3) for *β*, we obtain the expression for *β* given in (2.27).

This completes the proof, because the last calculation provides continuity of *θ*(*x*,*y*), and by construction *r* and *θ* satisfy the ordinary differential equations (2.24). ▪

- Received November 9, 2015.
- Accepted January 14, 2016.

- © 2016 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.