Periodic orbits in the Ott-Antonsen manifold

In their seminal paper [Chaos 18, 037113 (2008)], E. Ott and T. M. Antonsen showed that large groups of phase oscillators driven by a certain type of common force display low dimensional long-term dynamics, which is described by a small number of ordinary differential equations. This fact was later used as a simplifying reduction technique in many studies of synchronization phenomena occurring in networks of coupled oscillators and in neural networks. Most of these studies focused mainly on partially synchronized states corresponding to the equilibrium-type dynamics in the so called Ott-Antonsen manifold. Going beyond this paradigm, here we propose a universal approach for the efficient analysis of partially synchronized states with non-equilibrium periodic collective dynamics. Our method is based on the observation that the Poincar\'e map of the complex Riccati equation, which describes the dynamics in the Ott-Antonsen manifold, coincides with the well-known M\"obius transformation. To illustrate the possibilities of our method, we use it to calculate a complete bifurcation diagram of travelling chimera states in a ring network of phase oscillators with asymmetric nonlocal coupling.


Introduction
Synchronization of rhythmic processes is a fundamental dynamical mechanism that underlies the functioning of many natural and man-made systems [1,2,3,4]. Circadian clocks [5] and metachronal waves in cilia carpets [6,7], Josephson junction arrays [8] and power grids [9,10] are just a few examples of this kind. Other situations are also known where synchronization can occur, but is undesirable. These include different neurological disorders such as schizophrenia, epilepsy, Alzheimer's and Parkinson's disease [11,12]. The variety of the above examples was the motivation for the development of mathematical theory of synchronization [13,14], culminating in the method of master stability function [19] and numerous case studies for the paradigmatic Kuramoto model [15,16,17,18]. Initially, the main focus of research was on full synchronization [20], when all components of a system behave identically, but later more complex forms of synchronous collective dynamics were identified and studied, including clustered synchronization [21,22], generalized synchronization [23], phase synchronization of chaos [24] and self-organized quasiperiodicity [25].
In heterogeneous systems consisting of many non-identical oscillators, synchronization usually occurs as a dynamical aggregation process controlled by the intensity of the interaction between these oscillators. The course of the process depends on the properties of the oscillators and the nature of their interaction [26,27,28,29,30]. Roughly speaking, for each such system there is a critical coupling strength, above which the asynchronous state becomes unstable, while one or more frequency-synchronized clusters are formed. As the coupling strength increases, the clusters increase in size and merge together, so that a fully synchronized state is eventually achieved. The simplest and most detailed description of this phenomenon can be given using the Kuramoto model [15], where the state of each oscillator is represented by one scalar variable, its phase, and the oscillators are all-to-all coupled. More specifically, if ω j ∈ R denotes the natural frequency of the jth oscillator and K > 0 is the coupling strength, then the phases of the oscillators evolve according to A remarkable feature of model (1) is that in the thermodynamic limit N → ∞ its dynamics can be carefully investigated using the McKean-Vlasov equation (also called the continuity equation) for a distribution ρ(θ, ω, t) that yields the probability to find an oscillator θ j (t) ≈ θ with the natural frequency ω j ≈ ω at time t. Moreover, if the natural frequencies ω j are drawn from a Lorentzian distribution, then the long-term dynamics of ρ(θ, ω, t) converges to a low-dimensional manifold parameterized by the limit value (as N → ∞) of the global order parameter This manifold was first described by Ott and Antonsen in [31,32], and since then, it has been used to study various types of collective dynamics in systems of phase oscillators [28,33], Winfree oscillators [36], theta neurons [34,35], and quadratic integrate-and-fire neurons [37,38]. It should be noted that so far the full power of the Ott-Antonsen method has been demonstrated mainly for statistically stationary states. For example, in the Kuramoto model (1) such are asynchronous and partially synchronized states with nearly constant magnitudes |z(t)| of the order parameter (2). On the other hand, there are a number of cases where the Ott-Antonsen method allowed to detect more complex collective dynamics of oscillators and neurons, which has so far been studied only superficially. In this regard, we can mention non-stationary partially synchronized states with periodically oscillating magnitudes |z(t)|, which occur in the model (1) with a bimodal distribution of natural frequencies ω j [39], as well as other examples of periodically modulated collective dynamics, which were briefly reported in [28,34]. Motivated by the latter examples, in this paper we want to propose a mathematical approach for the efficient analysis of periodic dynamics in the Ott-Antonsen manifold. For this, in Section 2 we introduce two auxiliary models: (i) a system of phase oscillators driven by a common periodic force, and (ii) a system of theta neurons driven by a common periodic input current. For both models we write complex differential equations that determine the corresponding collective dynamics in the Ott-Antonsen manifold. In Section 3 we consider these equations and show that their Poincaré maps coincide with the well-known Möbius transformation. Next, in Section 4 we explain how this fact can be used for fast calculation of periodic orbits in the Ott-Antonsen manifold. A practical application of the proposed semi-analytical method is described in Section 5. There, we consider a nonlinear integro-differential equation of the Ott-Antonsen type, which describes the long-term dynamics of a ring network of nonlocally coupled phase oscillators with Lorentzian-distributed natural frequencies. We focus on the travelling wave solutions of this equation, which represent so called travelling chimera states [40], and show how these solutions can be quickly calculated using a kind of self-consistency argument. Finally, in Section 6 we discuss other problems, which can be considered by the method of this paper.

Two models and two Ott-Antonsen equations
The simplest and therefore the most popular mathematical models used in the study of synchronization phenomena include various types of networks consisting of phase oscillators and theta-neurons. Their collective dynamics is often investigated by means of the self-consistency method. Roughly speaking, one assumes that oscillators or neurons are influenced by a given external force and calculates the corresponding dynamics of each network's node. From this, the effective force of interaction between nodes is estimated and the obtained value is compared with the initially assumed value of the force. The resulting match relation is called the self-consistency equation and its analysis is often much easier than the analysis of the original network dynamics. Below, we describe two auxiliary models related to the application of the self-consistency method to networks of phase oscillators and networks of theta-neurons.
Phase oscillators. Using definition (2), the Kuramoto model (1) can be rewritten in the form as if each oscillator is driven by the global order parameter z(t) multiplied by K. Generalizing this setting, we can write another model: a population of phase oscillators driven by an arbitrary complex-valued force W (t) where θ j (t) ∈ R is still the phase of the jth oscillator and ω j ∈ R is its natural frequency. If the natural frequencies ω j are drawn randomly and independently from a Lorentzian distribution g(ω) = γ π 1 (ω − ω 0 ) 2 + γ 2 with ω 0 ∈ R and γ > 0, then using the Ott-Antonsen method [31,32] it can be shown [41,Sec. 3.1] that in the thermodynamics limit N → ∞ the long-term dynamics of system (3) is completely described by a scalar complex equation We call this equation the Ott-Antonsen equation corresponding to model (3). Theta neurons. The theta neuron is an excitable dynamical unit that is described by the normal form of a saddle-node-on-an-invariant-circle (SNIC) bifurcation [42,34]. The external driving usually acts on each neuron in the form of a real input current J(t) so that dθ j dt = 1 − cos θ j + (1 + cos θ j )(η j + J(t)), j = 1, . . . , N, where η j ∈ R is the excitability parameter of the jth neuron. In the case when η j are chosen from a Lorentzian distribution one can also apply the Ott-Antonsen theory and obtain the mean-field equation [41,Sec. 3 By analogy with (4), we call Eq. (6) the Ott-Antonsen equation corresponding to model (5).

Complex Riccati equation
It is easy to see that both Eq. (4) and Eq. (6) are particular cases of the more general complex Riccati equation Indeed, Eq. (7) coincides with Eq. (4), if we assume Similarly, Eq. (6) is obtained from Eq. (7) in the case Recall that in the Ott-Antonsen method, the function z(t), which solves Eq. (4) or Eq. (6), determines the limit value (as N → ∞) of the global order parameter (2), so by definition it must satisfy the inequality |z(t)| ≤ 1 for all t ≥ 0. Therefore, in the rest of the section we consider only such complex Riccati equations, which guarantee this property. For a more accurate statement, let us denote by D = {z ∈ C : |z| < 1} the open unit disc in the complex plane, by D = {z ∈ C : |z| ≤ 1} the closure of D, and by S = {z ∈ C : |z| = 1} the boundary of D. Then the next proposition provides a sufficient condition for Eq. (7) to be consistent with the Ott-Antonsen method.

Proposition 3.1 Suppose
Re(c 0 (t)z + c 1 (t) + c 2 (t)z) ≤ 0 for all z ∈ S and t ≥ 0, then the closed unit disc D is an invariant set of Eq. (7). In other words, if z(0) ∈ D, then the corresponding solution z(t) of Eq. (7) lies in D for all t > 0.
Now we focus on the periodic complex Riccati equation, i.e. Eq. (7) where c 0 (t), c 1 (t) and c 2 (t) are 2π-periodic continuous complex-valued functions. Our goal is to analyze all possible 2π-periodic solutions of this equation and determine the conditions that ensure the existence of a periodic solution that lies entirely in the unit disc D.
Let U (t, z 0 ) denote the solution of the initial value problem for Eq. (7) with the initial condition z(0) = z 0 . Then the mapping z 0 ∈ C → U (2π, z 0 ) ∈ C is called the Poincaré map of Eq. (7). For the periodic complex Riccati equation, it is known [43,44] that its Poincaré map is a one-to-one map of the extended complex planeĈ = C ∪ {∞} (also called the Riemann sphere) onto itself, which corresponds to the Möbius transformation with complex coefficients a, b, c and d such that ad − bc = 0. This correspondence provides a useful mathematical tool for analyzing the periodic solutions of Eq. (7). Indeed, every periodic solution of Eq. (7) corresponds to a fixed point of the transformation M(z). Since the Möbius transformation in general has two different fixed points (or one fixed point of multiplicity two in the degenerate case), the same can be said about the number of periodic solutions of Eq. (7). More detailed information about the position and stability of these fixed points can be obtained from the geometric properties of the Möbius transformation [45]. For this, one needs to consider the behaviour of map trajectories, i.e. complex sequences z n+1 = M(z n ), n = 0, 1, 2, . . . , with different initial conditions z 0 ∈ C. Note that every Möbius transformation is invertible. For instance, the inverse of (11) reads Therefore, the above map trajectory can be extended not only for increasing indices n but also for decreasing ones, namely Qualitative difference in the behaviour of map trajectories allows one to identify four main types of Möbius transformations: parabolic, elliptic, hyperbolic and loxodromic, see Fig. 1. Parabolic transform is the only type of Möbius transformation with one degenerate fixed point of multiplicity two. In this case, every map trajectory converges to this fixed point for both n → −∞ and n → +∞, e.g. Fig. 1(a). In contrast, all non-parabolic transforms have two different fixed points. More specifically, for an elliptic transform each map trajectory lies on a circle around one of the fixed points, e.g. Fig. 1(b). In the resonant case, the trajectory visits only a finite number of points of the circle, otherwise it fills the circle densely. For hyperbolic and loxodromic transforms, each map trajectory lies on a smooth curve that connects two fixed points, see Fig. 1(c) and (d) respectively. The trajectory converges to one of the fixed points for n → +∞ and to the other for n → −∞. The main difference between hyperbolic and loxodromic cases originates from the fact that in the former case the curves on which the different trajectories lie are circular arcs, while in the latter case they are logarithmic spirals. The location of the map trajectories in each of the above cases determines the stability of the fixed points of the corresponding Möbius transformation. For example, we can see that the degenerate fixed point of a parabolic transform is always unstable in the sense of Lyapunov. On the other hand, both fixed points of an elliptic transform are stable in the sense of Lyapunov, but not asymptotically stable. Finally, every hyperbolic or loxodromic transform has one asymptotically stable and one unstable fixed point. These facts are important for our further consideration, because they allow us to characterize the stability of periodic solutions of Eq. (7).
This means that every solution z(t) of Eq. (7) with z(0) ∈ S is trapped in the unit disc D, so that |z(t)| < 1 for all t > 0. Hence, the Möbius transformation M(z) representing the Poincaré map of such Eq. (7) satisfies |M(z)| < 1 for all z ∈ S, and therefore M(D) ⊂ D. Due to [46,Theorem 4.5], the latter is a sufficient condition for M(z) to be a contraction on D.
Then, Theorem 1.1 and Theorem 3.6 from [46] imply that M(z) is a hyperbolic or loxodromic Möbius transformation 1 , that its stable fixed point lies in D, and that its unstable fixed point lies outside D. Obviously, if we take the stable fixed point of M(z) as the initial condition z(0) for Eq. (7), we obtain a 2π-periodic solution z(t) such that |z(t)| < 1 for all 0 ≤ t ≤ 2π.
Remark 3.4 Note that the strict inequality in (12) is essential for the statement of Proposition 3.3. For example, if instead of (12) we use the non-strict inequality from Proposition 3.1, we cannot guarantee that the Poincaré map of Eq. (7) is represented by a hyperbolic or loxodromic Möbius transformation only. The corresponding counter-examples can be found in [47].

Remark 3.5 A statement similar to Proposition 3.3 can also be formulated in the case
Then, the Poincaré map of Eq. (7) is again described by a hyperbolic or loxodromic Möbius transformation. But, the stable fixed point of this map lies outside D, while the unstable fixed point lies in D. To see this, it is enough to use the substitution z(t) =z(−t) in Eq. (7) and consider the resulting equation as an equation with respect toz(t).
Therefore, the strict inequality (12) holds for all z ∈ S, except z = −1. To overcome this difficulty, we propose a modified version of Proposition 3.3.
Proof: We only need to show that M(S) ⊂ D, then we can repeat the remaining arguments of the proof of Proposition 3.3. For this, we consider a general solution z(t) of Eq. (7) with z(0) ∈ S. In this case, Proposition 3.1 ensures that z(t) ∈ D for all t > 0. Could it be that z(2π) ∈ S? The answer is no, because d|z| 2 dt = 2Re c 0 (t)z + c 1 (t) + c 2 (t)z < 0 for all z ∈ S\z * and t = 2π, and because of the assumption that the solution of Eq. (7) with z(2π) = z * satisfies z(0) / ∈ S. This ends the proof.

Solution of periodic complex Riccati equation
The complex Riccati equation (7) is a nonlinear differential equation that cannot be solved analytically. Therefore, its periodic solutions can usually be found only by numerical methods, such as the shooting method or the collocation method. However, without a good initial guess, each of these methods may involve a large number of iterations and thus require a long computational time to ensure the desired accuracy of the result. A more efficient way to calculate the periodic solutions of Eq. (7) can be proposed using the relation between the Möbius transformation and the Poincaré map of Eq. (7).
Suppose that the assumption of Proposition 3.3 is satisfied. Then, by choosing three distinct points z k ∈ D, k = 1, 2, 3, and solving Eq. (7) with different initial conditions z(0) = z k , we obtain three complex functions U k (t). Due to Proposition 3.1 each of these functions is bounded and satisfies |U k (t)| ≤ 1. Let us denote w k = U k (2π), then where M(z) is a Möbius transformation of the form (11) that represents the Poincaré map of Eq. (7). It is well-known [45] that the information contained in the relations (13) is sufficient to determine the coefficients a, b, c and d in the expression of M(z). The corresponding explicit formulas read Once the transformation M(z) is defined, its fixed points can be found by solving the equation The latter is obviously equivalent to the quadratic equation and has two roots Now both periodic solutions of Eq. (7) can be obtained by integrating this equation with the initial conditions z(0) = z − and z(0) = z + . Importantly, Proposition 3.3 ensures that one of the fixed points z − or z + lies in the open unit disc D and the corresponding periodic solution satisfies |z(t)| < 1.

Remark 4.1
In some cases, the map trajectories of the Möbius transformation M(z) converge to its fixed point in D so fast that the numbers w 1 , w 2 and w 3 lie extremely close to each other. Then, the above formulas for the coefficients a, b, c and d cannot be applied, due to the vanishing values of the determinants. In this case, the fixed point of M(z) that lies in D can simply be approximated by the mean (w 1 + w 2 + w 3 )/3.

Application to travelling chimera states
In [40], the author of this paper considered moving coherence-incoherence patterns, called travelling chimera states, in a ring network of nonlocally coupled phase oscillators. In the continuum limit of infinitely many oscillators, the long-term dynamics of such a network is described by an integro-differential equation [28,48] where z(x, t) is the unknown complex-valued function satisfying the periodic boundary condition z(x + 2π, t) = z(x, t), γ > 0 is a parameter analogous to the width of the Lorentzian distribution in Section 2, and is a convolution integral operator with a non-constant real 2π-periodic kernel G(x). Moreover, each travelling chimera state corresponds in Eq. (14) to a travelling wave solution of the form where a(x) ∈ D is the wave profile, s ∈ R is the wave speed and Ω ∈ R is its complex-phase velocity. The analysis of travelling waves (15) carried out in [40] revealed unexpected oscillatory properties of their profiles a(x) and a non-monotonic dependence of the speed s on the system parameters. Below, we extend this analysis, using the advantages provided by the semi-analytic description of the periodic solutions of the complex Riccati equation (7).

Self-consistency equation
In this section, we derive a self-consistency equation for travelling waves (15), which helps to significantly speed up their numerical calculation. Inserting ansatz (15) into Eq. (14) we obtain where a(x) is the unknown complex-valued function and a is its usual derivative. Next, subtracting iΩa in the both sides and dividing the resulting equation by (−s) = 0, we get The latter equation can be written in the form where This is the complex Riccati equation (7) with c 0 (x) = w(x), c 1 (x) = ζ and c 2 (x) = −w(x).
Obviously, for γ > 0 and s = 0 we have therefore for every w ∈ C per ([0, 2π]; C) and ζ ∈ {z ∈ C : Re z = 0} there exists a unique 2π-periodic solution of Eq. (17) that lies entirely in the unit disc D, see Proposition 3.3 and Remark 3.5. (Note that although for s > 0 the above solution a(x) is unstable with respect to Eq. (17), the corresponding travelling wave (15) may be stable or unstable with respect to Eq. (14).) In the following, we denote the obtained solution operator for Eq. (17) by U(w, ζ). More precisely, this operator determines a mapping (Note that in the above definition we have |a(t)| < 1 for all t ∈ [0, 2π]!) If a = U(w, ζ), then to agree with the formulas (18) we must have Remark 5.1 Every travelling wave solution of Eq. (14) is not uniquely determined. Given a wave profile a 0 (x), one obtains infinitely many other solutions of Eq. (14) by the formula a(x) = a 0 (x−ξ)e iφ with ξ, φ ∈ R. Moreover, it is easy to verify that the same symmetry property is inherited by Eq. (19). In other words, if w 0 (x) solves Eq. (19), then w(x) = w 0 (x − ξ)e iφ with ξ, φ ∈ R is a solution of Eq. (19) too.
Eq. (19) is the self-consistency equation to be solved with respect to w ∈ C per ([−π, π]; C), s ∈ R\{0} and Ω ∈ R. A unique solution of Eq. (19) is fixed by two supplementary pinning conditions

Self-consistency equation for a trigonometric coupling function
Now we consider travelling wave solutions of Eq. (14) in the case of a trigonometric coupling function with two real coefficients A and B. The integral operator G with such a kernel G(x) is a finite-rank operator. This follows from its representation formula where denotes the usual L 2 -scalar product, and are three basis functions that span the image of the operator G. Formula (23) implies that every solution of Eq. (19) with the above trigonometric coupling function G(x) has the form w(x) =ŵ 0 +ŵ 1 e ix +ŵ 2 e −ix (24) whereŵ 0 ,ŵ 1 ,ŵ 2 ∈ C. Moreover, if the pinning conditions (20) and (21) are satisfied, thenŵ 0 andŵ 1 must be real. Note that due to Remark 5.1, we can always achieve this by using two continuous symmetries of the solutions of Eq. (19). Indeed, suppose that w(x) is a solution of Eq. (19) in the most general form (24). Using the complex-phase-shift symmetry, we can obtain another solution with real first coefficient. Then, using the translation symmetry we can get a solution of the form which has real first and second coefficients and thus satisfies (20) and (21).
In the rest of this section we show how the system (19)- (21), which determines travelling wave solutions of Eq. (14), can be reduced to a finite-dimensional nonlinear system. For this, we look for solutions of Eq. (19) in the following form where p, q ∈ R andŵ * ∈ C. Such ansatz ensures that the pinning conditions (20) and (21) are satisfied automatically. Then, from formula (23) and from the linear independence of functions ψ k (x) it follows that Eq. (19) is equivalent to a system of three scalar complex equations: This system can be solved with respect to six unknowns p, q, Reŵ * , Imŵ * , s and Ω, using a standard Newton's method. Given a system's solution, the corresponding wave profile can be calculated by the formula a(x) = U w(x), γ + iΩ s .

Results
Using system (25)- (27), we carried out a detailed analysis of travelling wave solutions of Eq. (14) in the case of trigonometric coupling (22) with B = 0. We kept the parameters A = 0.9 and α = π/2 − 0.1 fixed and used a pseudo-arclength continuation to follow three branches of travelling waves for γ = 0.005, 0.01 and 0.02. For each travelling wave we calculated its twist, which is defined as the net number of multiples of 2π through which the argument of the complex wave profile a(x) decreases as the spatial domain is traversed once. Moreover, the stability of travelling waves was determined by the linearization of Eq. (14) about the found solution, which led to the consideration of the eigenvalue problem [40] λv where a(x), s and Ω are the wave profile, speed and complex-phase velocity of the reference travelling wave, λ and (v + (x), v − (x)) T are the eigenvalue and eigenfunction, and The above integral system was discretized on a uniform grid of 1000 points and solved as a matrix eigenvalue problem. Note that in contrast to the method based on the Lyapunov-Schmidt reduction of Eq. (16) and used in [40] to calculate the travelling wave solutions of Eq. (14), the numerical method proposed in this paper allowed us to perform the same calculations much faster. Therefore, we could calculate complete solution branches, also in the case of values γ smaller and larger than the value considered in [40]. Next we describe the obtained results. For γ = 0.01, the branch of travelling waves starts out as twist-0 for B ≈ 0, see Fig. 2. It remains twist-0 for small speeds s, but then it becomes twist-1 for s ∈ (0.0665, 0.0683), twist-2 for s ∈ (0.0683, 0.0749), and twist-3 for s ∈ (0.0749, 0.0827), as explained in Fig. 3. For larger values of s, the inverse order of transformations is observed: from twist-3 into twist-2 and then into twist-1. It is noteworthy that stable travelling waves occur with relatively low speeds only, while all waves with s > 0.0762 are unstable. Moreover, the endpoint of the branch, i.e. the wave with the largest speed s max ≈ 0.42, corresponds to a splay state z(x, t) = a max e −i(x−smaxt) e iΩmaxt = a max e −ix e i(Ωmax+smax)t (30) with a max ≈ 0.82 and Ω max ≈ −0.37. Remark that the exponent in the right-hand side of formula (30) is our motivation to show the sum Ω + s instead of the complex-phase velocity Ω in Figs. 2, 3 and later.  Fig. 2. The data are identical to those shown in Fig. 4 in [40]. But here they were computed with the new continuation method based on Eq. (19). Other notations the same as in Fig. 2.
The situation that the branch of travelling waves occurs as a "bridge" between a standing wave (s = 0) and some splay state is also observed for larger and for smaller values of γ, see Figs. 4 and 5. However, the smaller is γ the more oscillatory is the corresponding s-versus-B diagram. Moreover, the maximal twist found on the branch of travelling waves typically grows for decreasing γ. For example, for γ = 0.02 we find only twist-0, twist-1 and twist-2 waves. In contrast, for γ = 0.005 we find travelling waves with twists up to 5.

Discussion
Let us summarize the main results of this paper. In Section 3 we have outlined a wide class of periodic complex Riccati equations, whose Poincaré maps are represented by hyperbolic or loxodromic Möbius transformations, which map the closed unit disc D into itself. We have shown that this class includes several types of the Ott-Antonsen equation used in the analysis of the dynamics of phase oscillator networks and networks of theta neurons. In Section 4 we explained how the properties of Möbius transformations can be exploited to calculate periodic solutions of the complex Riccati equation by solving at most four intial value problems for this equation. Finally, a practical application of the latter fact was demonstrated in Section 5, where we derived the self-consistency equation for travelling chimera states arising in a ring network of phase oscillators with asymmetric nonlocal coupling, and used it to calculate several complete branches of such states. It is likely that the semi-analytical method of this paper can be generalized to study more complex coherence-incoherence patterns in large networks of phase oscillators, including breathing, pulsating and alternating chimera states [28,33,49], as well as moving chimera states on two-and three-dimensional oscillator lattices [50,51,52]. Another class of potential applications is concerned with Proposition 3.6, which can be useful in the study of moving and oscillatory bump states in theta neuron networks [53,54]. Moreover, recalling one of the recent applications of the Ott-Antonsen method to networks of quadratic integrate-and-fire neurons [37,38,55], we expect that our method can also be adapted for such systems as well. We plan to report on these issues in future work.