Fast phase randomization via two-folds

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.


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][11][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 1a shows a histogram of the phase φ T , computed relative to a reference time T and limiting to the asymptotic phase as T → ∞. This shows an apparently uniform distribution, hence the phase of solutions from the initial point is 'fully randomized'. Figure 1b 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.  Figure 1. Each panel shows a histogram of the phase φ T , of 10 4 sample solutions from the same initial condition, of the general stochastic system (4.1). The two panels correspond to different choices for the higher-order terms that affect the global dynamics of (4.1). One sample solution corresponding to each panel is shown in figure 7 (see caption of this figure for the parameter values used). The solid curves show our theoretical approximation for the probability density function of φ T , see §4b. (Online version in colour.) 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.

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 piecewisesmooth systemẊ and V − , V + ∈ R are parameters [16,17]. The two-fold is located at the origin, 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. x y z C + C -(0, y, gy) (0, µy, µ gy) Figure 2. A schematic of the two-fold of (2.1). The discontinuity surface x = 0 is made up of an attracting sliding region A (y, z > 0), a repelling sliding region R (y, z < 0), and two crossing sliding regions, C + (y < 0, z > 0) and C − (y > 0, z < 0). Parts of two orbits and their intersections with the discontinuity surface are shown. One orbit starts in R and has subsequent intersections in C + and C − . A second orbit starts and returns to an invariant line ζ ∈ C − ; see §2a. (Online version in colour.) The system is solved across the discontinuity by forming the Filippov system for (2.1). Using e i to denote the ith coordinate vector in R 3 , a Filippov solution of (2.1) is an absolutely continuous function ϕ(t) that satisfies for almost all values of 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ẋ = 0 gives the system that sliding motion obeys on A and R as Various details of the sliding flow can be derived from (2.4); see [16,17]. Throughout this paper, we assume (2.5) in that case a typical orbit of (2.3) and (2.4) will intersect the two-fold in forward time or backward time or both, and give the geometry shown in figure 2. With (2.5), the matrix in (2.4) has distinct negative eigenvalues. Of these, the eigenvalue that is smallest (in absolute value) is λ = 1 2 (V + + V − + [(V + − V − ) 2 + 4] 1/2 ). As orbits of (2.1) slide into the singularity inside 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),Ẋ = F L (X) andẊ = F R (X), have solutions and 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 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 − 2y when it is located at (0, −y, −2V + 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).

(Online version in colour.)
In the full three-dimensional space of (2.1), this eigenvector corresponds to a ray that emanates from the two-fold, as shown in figure 2. By evolving the flow forwards from ζ according to (2.7) and (2.8), we obtain a surface Λ that is given implicitly by Figure 3 shows a plot of Λ. 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 By performing this iteration backwards through repeated crossings, we can construct a Filippov solution that approaches the two-fold in backward time. This Filippov solution emanates from the two-fold and has infinitely many intersections with ζ 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 n ∈ Z. It follows that each ψ μ n a (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 The intersection of Λ 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 (x = r cos(θ), y = r sin(θ)) bear no relation to the behaviour of the system and yield complicated expressions that provide no insight. Instead we derive radial and angular 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 ∈ R. We let 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 z ∈ R, without again intersecting 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 We have ϕ L (τ L ; 0, y 0 , γ y 0 , 0) = (x, y, z), and so, by (2.7), We solve the first two of these equations simultaneously for y 0 and τ L . Eliminating y 0 yields to which an application of the quadratic formula gives Upon substituting (2.12) into (2.20), substantial simplification is possible and we obtain Moreover, y 0 is given by Second, choose any x > 0 and y ∈ R. We let 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 z ∈ R, without first re-intersecting 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 τ R = y − y 2 + 2x, (2.22) and 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. r(0, y) = y, θ (0, y) = 0, for all y > 0. (2.23) The positive 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θ to be inversely proportional to 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ṙ should be constant. In summary, we would like to define r and θ so that they satisfy (2.23) andṙ for some constants α, β > 0. The following choice of r and θ achieves this.

26)
where τ L and τ R are given by (2.21) and (2.22), and 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 θ = 2π ln((1 − 2α)μ)/ln(μ). Also by (2.27), we have 0 < α < 1 2 , with α → 0 as V − V + → 1 and α → 1 2 as V − V + → ∞ (assuming 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 for some C ∈ [0, 2π ). Therefore, the Filippov solution ψ a (t), defined in §2a, may be written in polar coordinates as 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.
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.

Deterministic global dynamics
The vector field near a generic invisible two-fold in a three-dimensional system can be transformed toẊ where F L and F R are the constituents of the normal form (2.2), and represent higher-order terms. As described in [16], there exists a surfaceΛ (that coincides with Λ in the limit |X| → 0) on which viable Filippov solutions to (3.1) evolve. This surface intersects x = 0 with y > 0 on a curvẽ ζ that matches ζ to first order. We letψ a (t) denote a viable Filippov solution to (3.1) withψ a (0) = 0 andψ a (a) ∈ζ .
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Λ approach Γ . 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φ(t; X 0 , t 0 ) be a Filippov solution to (3.1) that limits to Γ as t → ∞. The phase ofφ, relative to a reference time T that is assumed to be sufficiently large that the orbit lies close to Γ , is defined as Note that 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Λ. Figure 5 shows the intersection of five isochrons withΛ, produced using The surfaceΛ was computed by fitting a mesh to 1000 numerically computed forward orbits. The stable periodic orbit Γ 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 twofold 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ζ . First note thatψ a (a) ∈ζ by definition. Then, for any a > 0, let t = f (a) > a be the next time at whichψ a (t) ∈ζ . Figure 6 shows a plot of 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,ψ a (t) is located near Γ 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 Γ .

Stochastic dynamics and phase randomization
Here, we study a stochastic perturbation of (3.1) defined by the three-dimensional stochastic differential equation where 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; X 0 , 0), and a time 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

(a) Sample solutions and phase randomization
A sample solution to (4.1) using (2.6) and (3.5) is shown in figure 7a. The initial point is X = (0, 1, 1), and we denote this solution asφ ε (t; 0, 1, 1, 0). The deterministic solutionφ(t; 0, 1, 1, 0) slides into the two-fold in a time t 0 = 3.0445 (to four decimal places). The perturbed solutionφ ε (t) initially follows a random path close toφ(t) near the discontinuity surface. Then while t ≈ t 0 , the perturbed solutionφ ε (t) is located near the two-fold. It is during this brief period that the order-ε noise has an order-1 influence on the location and phase ofφ ε (t) at later times. Nextφ ε (t) spirals outwards. The initial portion of this spiralling motion is close to the two-fold, and so the higher-order terms G L and G R have little influence. At later times,φ ε (t) is located relatively far from the two-fold.
Here G L and G R are important andφ ε (t) approaches the stable periodic orbit of (3.1), Γ .
We computed 10 4 sample solutions analogous to the solution shown in figure 7a (i.e. with the same initial point and parameter values). Figure 1a 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 7b shows a sample solution using (4.2) instead of (3.5). As in figure 7a, 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 1b 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 X 0 ∈ R 3 be such that the deterministic forward orbit of this pointφ(t; X 0 , 0) is located at the two-fold at some time t 0 > 0. For sample solutionsφ ε (t; X 0 , 0), we are only interested in the phase φ 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φ ε (t) will lie near the two-fold with high probability. Thus for the purposes of computing 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π ).
As sample solutionsφ ε (t) continue to evolve outwards from the two-fold, because there are no other singular points for solutions to encounter, the noise only has the effect of diffusing intersection times 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 Iterating 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 ln(a) then our procedure generates a uniform distribution for p T . Indeed, in figure 8a, f n is indistinguishable from an affine function on the given scale and so, in figure 1a, the distribution is approximately uniform. In figure 8b, f n has a notable nonlinearity, and, for this reason, the distribution in figure 1b is significantly non-uniform.

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. x y t y = a 3 t Figure 9. The dynamics of (5.2) for t 1 < t < t 2 (when the control is active) in the absence of noise using (5.5), t 1 = −5, and t 2 = 2.5. The evolution of x 2 + y 2 = 1-the periodic orbit of (5.1)-is indicated for t 1 < t < 0. As we evolve the set x 2 + y 2 = 1, it contracts to the discontinuity surface at t ≈ −3.7, then to the two-fold at t = 0. For 0 < t < t 2 , viable Filippov solutions emanating from the two-fold evolve onΛ. (Online version in colour.) To illustrate our method, we assume the oscillators are governed by the planar system (ẋ,ẏ) = F(x, y) = (x − y − x(x 2 + y 2 ), x + y − y(x 2 + y 2 )), (5.1) and apply a discontinuous control that creates a two-fold singularity. In principle, this can be achieved with any system exhibiting a stable periodic orbit, including excitable systems such as the FitzHugh-Nagumo equations [8].
Consider the two-dimensional stochastic differential equation where 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 where a 1 , a 2 , a 3 , a 4 ∈ R are control parameters, is sufficient for our purposes. An identification of more sophisticated functions providing superior performance is beyond the scope of this paper. By treating the time t as a third dynamic variable, i.e. withṫ = 1, the system is threedimensional and is piecewise-smooth when t 1 < t < t 2 with the discontinuity surface x = 0. The fold lines on x = 0 (whereẋ = 0) are 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 twofold 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  to be particularly effective. With smaller values of the 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Λ is shown for 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.

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 twofold, 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.