Synchronization Behavior in a Ternary Phase Model

TL;DR

This paper investigates a ternary phase model with quenched disorder, revealing first-order synchronization transitions with hysteresis, contrasting with Kuramoto model behavior, supported by theoretical analysis of the infinite-oscillator limit.

Contribution

It introduces and analyzes a novel ternary coupling phase model with quenched disorder, providing insights into its synchronization transitions and theoretical predictions in the infinite-oscillator limit.

Findings

First-order transitions with hysteresis observed for Gaussian and uniform disorder.

Theoretical predictions derived from the infinite-oscillator limit.

Comparison shows stark differences from Kuramoto model behavior.

Abstract

Localized traveling-wave solutions to a nonlinear Schrodinger equation were recently shown to be a consequence of Fourier mode synchronization. The reduced dynamics describing mode interaction take the form of a phase model with novel ternary coupling. We analyze this model in the presence of quenched disorder and explore transitions to partial and complete synchronization. For both Gaussian and uniform disorder, first-order transitions with hysteresis are observed. These results are compared with the phenomenology of the Kuramoto model which exhibits starkly different behavior. An infinite-oscillator limit of the model is derived and solved to provide theoretical predictions for the observed transitions. Treatment of the nonlocal ternary coupling in this limit sheds some light on the model's novel structure.

Figures

Equations

$$\begin{array}{rl}{\dot{\Delta }}_j& =\frac{K}{2N}\sum _{l=j-N}^N\sin ({\Delta }_l+{\Delta }_{j-l}-{\Delta }_j),;j=1,2,\dots ,N\end{array}$$

$${\Delta }_j=\alpha \cdot j$$

$${\theta }_j=\beta .$$

$${\dot{\Delta }}_j={\omega }_j+\frac{K}{2N}\sum _{l=j-N}^N\sin ({\Delta }_l+{\Delta }_{j-l}-{\Delta }_j)$$

$$V=-\frac{K}{6N}\sum _{j=1}^N\sum _{l=j-N}^N\cos ({\Delta }_l+{\Delta }_{j-l}-{\Delta }_j).$$

$$Z_j=\frac{1}{2N+1-j}\sum _{l=j-N}^N{\text{e}}^{\text{i}({\Delta }_l+{\Delta }_{j-l})}.$$

$${\Delta }_j\to \alpha \cdot j+{\delta }_j+{\Omega }_j\cdot t.$$

$${\Omega }_j=\Omega \cdot j.$$

$$\sum _{j=1}^N{\dot{\Delta }}_j\cdot j=\sum _{j=1}^N{\omega }_j\cdot j$$

$$\sum _{j=1}^N(\Omega \cdot j)\cdot j=\sum _{j=1}^N{\omega }_j\cdot j$$

$$\Omega =\frac{6}{N(N+1)(2N+1)}\cdot \sum _{j=1}^N{\omega }_j\cdot j.$$

$${\Delta }_j^{\prime }={\Delta }_j-\Omega \cdot j\cdot t.$$

$$\Gamma =\sum _{j=1}^N{\Delta }_j^{\prime }\cdot j.$$

$${\omega }_j=\frac{c}{N^3}\cdot j^3$$

$$\begin{array}{rl}{\dot{\Delta }}_j& ={\omega }_j+\frac{2N+1-j}{2N}K{R}_j\sin ({\psi }_j-{\Delta }_j),\\ Z_j& =R_j{\text{e}}^{\text{i}{\psi }_j}.\end{array}$$

$$Z_j=\frac{1}{2N+1-j}\sum _{l=j-N}^N{\text{e}}^{\text{sgn}(l)\text{i}{\Delta }_{|l|}}{\text{e}}^{\text{sgn}(j-l)\text{i}{\Delta }_{|j-l|}}$$

$$\begin{array}{r}{Z}_j=\frac{1}{2N+1-j}\sum _{l=j-N}^N\underset{0}{\overset{2\pi }{\int }}\text{d}\Delta \underset{0}{\overset{2\pi }{\int }}\text{d}{\Delta }^{\prime }\underset{-\infty }{\overset{\infty }{\int }}\text{d}\omega \underset{-\infty }{\overset{\infty }{\int }}\text{d}{\omega }^{\prime }\\ \times P(\Delta ,{\Delta }^{\prime },\omega ,{\omega }^{\prime },|l|,|j-l|){\text{e}}^{\text{sgn}(l)\text{i}\Delta }{\text{e}}^{\text{sgn}(j-l)\text{i}{\Delta }^{\prime }}\end{array}$$

$$P(\Delta ,{\Delta }^{\prime },\omega ,{\omega }^{\prime },|l|,|j-l|)=\rho (\Delta ,\omega ,|l|)\rho ({\Delta }^{\prime },{\omega }^{\prime },|j-l|)$$

$$\begin{array}{r}x=\frac{j-1}{N-1},\\ 0\le x\le 1.\end{array}$$

$$\frac{1}{2N+1-j}\sum _{l=j-N}^N\text{d}j=\frac{N-1}{2N-(N-1)x}\sum _{y=\frac{(N-1)x-N}{N-1}}^1\text{d}y.$$

$$\begin{array}{r}Z(x)=\frac{1}{2-x}\underset{x-1}{\overset{1}{\int }}\text{d}y\underset{0}{\overset{2\pi }{\int }}\text{d}\Delta \underset{0}{\overset{2\pi }{\int }}\text{d}{\Delta }^{\prime }\underset{-\infty }{\overset{\infty }{\int }}\text{d}\omega \underset{-\infty }{\overset{\infty }{\int }}\text{d}{\omega }^{\prime }\\ \times \rho (\Delta ,\omega ,|y|)\rho ({\Delta }^{\prime },{\omega }^{\prime },|x-y|){\text{e}}^{\text{sgn}(y)\text{i}\Delta }{\text{e}}^{\text{sgn}(x-y)\text{i}{\Delta }^{\prime }}.\end{array}$$

$$\begin{split}&\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial\Delta}\big{(}\rho\dot{\Delta}\big{)}=0,\\ \dot{\Delta}=\omega+&(1-x/2)KR(x)\sin(\psi(x)-\Delta),\\ &Z(x)=R(x)\textrm{e}^{\textrm{i}\psi(x)}\end{split}$$

$$\begin{array}{rl}& \underset{0}{\overset{2\pi }{\int }}\text{d}\Delta \underset{-\infty }{\overset{\infty }{\int }}\text{d}\omega \rho (\Delta ,\omega ,x)=1,\\ & \underset{0}{\overset{2\pi }{\int }}\text{d}\Delta \rho (\Delta ,\omega ,x)=g(\omega )\end{array}$$

$$R=\int\limits_{0}^{1}\textrm{d}y\Bigg{[}\frac{(2-y)KR}{2}\int\limits_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\textrm{d}\theta\;g\bigg{[}\frac{(2-y)KR}{2}\sin\theta\bigg{]}\cos^{2}\theta\Bigg{]}^{2},$$

$${\rho }_0=\frac{g(\omega )}{2\pi }.$$

$$\begin{array}{rl}\dot{\Delta }& =\omega -M(x)\sin \Delta ,\\ M(x)& =(1-x/2)KR\end{array}$$

$$\begin{split}&\rho(\Delta,\omega,x)=\delta\bigg{(}\Delta-\arcsin\frac{\omega}{M(x)}\bigg{)}g(\omega),\\ &\quad|\omega|\leq M(x).\end{split}$$

$$\begin{split}&\rho(\Delta,\omega,x)=\frac{a}{|\dot{\Delta}|}g(\omega)=\frac{a}{\big{|}\omega-M(x)\sin\Delta\big{|}}g(\omega),\\ &\quad|\omega|>M(x).\end{split}$$

$$\begin{split}R=&I_{1}+I_{2},\\ I_{1}=&\frac{1}{2}\int\limits_{-1}^{1}\textrm{d}y\int\limits_{0}^{2\pi}\textrm{d}\Delta\int\limits_{0}^{2\pi}\textrm{d}\Delta^{\prime}\int\limits_{-M(|y|)}^{M(|y|)}\textrm{d}\omega\int\limits_{-M(|y|)}^{M(|y|)}\textrm{d}\omega^{\prime}g(\omega)g(\omega^{\prime})\\ &\times\textrm{e}^{\textrm{sgn}(y)\textrm{i}\Delta}\textrm{e}^{-\textrm{sgn}(y)\textrm{i}\Delta^{\prime}}\\ &\times\delta\bigg{(}\Delta-\arcsin\frac{\omega}{M(|y|)}\bigg{)}\delta\bigg{(}\Delta^{\prime}-\arcsin\frac{\omega^{\prime}}{M(|y|)}\bigg{)},\\ I_{2}=&\frac{1}{2}\int\limits_{-1}^{1}\textrm{d}y\int\limits_{0}^{2\pi}\textrm{d}\Delta\int\limits_{0}^{2\pi}\textrm{d}\Delta^{\prime}\int\limits_{|\omega|>M(|y|)}\textrm{d}\omega\int\limits_{|\omega^{\prime}|>M(|y|)}\textrm{d}\omega^{\prime}g(\omega)g(\omega^{\prime})\\ &\times\textrm{e}^{\textrm{sgn}(y)\textrm{i}\Delta}\textrm{e}^{-\textrm{sgn}(y)\textrm{i}\Delta^{\prime}}\\ &\times\frac{a}{\big{|}\omega-M(|y|)\sin\Delta\big{|}}\cdot\frac{a}{\big{|}\omega^{\prime}-M(|y|)\sin\Delta^{\prime}\big{|}}.\end{split}$$

$$\begin{split}R=&\frac{1}{2}\int\limits_{-1}^{1}\textrm{d}y\int\limits_{-M(|y|)}^{M(|y|)}\textrm{d}\omega\int\limits_{-M(|y|)}^{M(|y|)}\textrm{d}\omega^{\prime}g(\omega)g(\omega^{\prime})\\ &\times\exp\bigg{[}\textrm{i}\;\textrm{sgn}(y)\arcsin\frac{\omega}{M(|y|)}\bigg{]}\\ &\times\exp\bigg{[}-\textrm{i}\;\textrm{sgn}(y)\arcsin\frac{\omega^{\prime}}{M(|y|)}\bigg{]}.\end{split}$$

Full text

††preprint: AIP/123-QED

The damped, driven nonlinear Schrödinger equation (NLSE) with cubic nonlinearity is ubiquitous in nonlinear science. This equation models a broad range of physical systems, for instance charge-density waves, Josephson junctions, ferromagnets in microwave fields, quantum Hall ferromagnets, radio-frequency-driven plasmas, shear flows in liquid crystals, ocean and atmospheric waves, and spatial and temporal nonlinear waves in optical resonators with Kerr nonlinearity. These nonlinear wave phenomena can all exhibit sharply-peaked patterns of spatially and/or temporally localized pulses. In these systems, pulsation can be a signature of synchronization between many oscillating degrees of freedom. A previously derived modelTaheri et al. (2017a) explicitly relates the existence of such a pulse to the synchronization of its Fourier modes. Here, we analyze this model in the context of general coupled-oscillator systems and investigate its synchronization behavior.

I Introduction

Mode locking, i.e. alignment of laser mode frequencies on an equidistant spectral grid, has long been known to be a requirement for pulsation in actively and passively mode-locked lasers Mcduff and Harris (1967); Kuizenga and Siegman (1970); Sargent (1974). Variants of the nonlinear Schrödinger equation (NLSE) that admit sharply peaked soliton pulses have successfully been used to model many types of mode locked lasers Haus (1975); Haus, Fujimoto, and Ippen (1991); Haus (2000); Martinez, Fork, and Gordon (1984). The connection between mode locking (pulsation in the time domain) and the universal phenomenon of self-synchronization is evidenced by the alignment and synchronous oscillation of Fourier modes in the frequency domain, see, e.g., (Haken, 1985; Wen et al., 2016; Taheri et al., 2017b). Recently, this nexus was formally established through a phase model underlying spontaneous pulse formation in the damped and driven NLSE Taheri et al. (2017a), also called the Lugiato-Lefever equation (LLE) Lugiato and Lefever (1987).

Interestingly, this phase model is of a novel type within the realm of self-synchronization phenomena: the phases interact in triples rather than pairs. The ternary coupling is rooted in the physics of the nonlinear interaction; specifically, the conservation of energy and momentum in nonlinear four-wave mixing in optical Kerr resonators. In this work, we treat the ternary phase model as a stand-alone model worthy of independent investigation from a nonlinear dynamical perspective. In particular, we explore the transitions to partial and complete synchronization, and compare the basic phenomenology to that of the Kuramoto model, the archetype of pairwise-interacting oscillator populations.

The paper is organized as follows. In Section II we introduce the ternary phase model and briefly review its relationship to the NLSE. In Section III, we present the results of numerical simulations showing the transition behavior for different choices of disorder, and compare this with the corresponding behavior of the Kuramoto model. In Section IV we construct a “thermodynamic limit” of the model and derive an expression we then use to calculate the order parameter quantifying synchronization. We end with a discussion of open questions.

II Background

In this paper, we explore a disordered version of the following coupled oscillator equation Taheri et al. (2017a)

[TABLE]

where $\Delta_{j}$ is the phase of the $j^{th}$ oscillator, and $K$ is the coupling constant. By definition, $\Delta_{-j}=-\Delta_{j}$. The equation bears some resemblance to the Kuramoto model Kuramoto (1975). We briefly review how Eq.(1) comes about; full technical details can be found in Ref. (Taheri et al., 2017a).

Equation 1 governs the slow dynamics of the LLE in the strongly pumped limit. As first noted by Wen et al, in this limit, the dynamical evolution of the system proceeds in three successive steps, which can be captured by a systematic expansion in inverse powers of the pump amplitudeWen et al. (2016). On the fastest time scale, the pumped mode steady-state is established. On the intermediate time scale, mode pairs (symmetrically located about the pump frequency) grow via parametric amplification. Each mode pair evolves independently, reaching its steady state amplitude, and simultaneously a partial ordering of its phases such that $\phi_{j}+\phi_{-j}-2\phi_{0}=0$. (Indices are chosen such that the pump index is zero.) It is during this intermediate time scale that a dense frequency comb is established; however, localized/pulsing dynamics emerges only on a still slower time scale, described by Eq.(1), during which the anti-symmetric parts of the phase pairs $\Delta_{j}=\phi_{j}-\phi_{-j}$ interact with each other. In the context of the LLE, $\Delta_{j}$ specifies a phase difference between modes. In this paper, we consider the set to be absolute phases of an oscillator population. The three-phase coupling of (1) is the direct consequence of the cubic nonlinearity in the underlying wave equation. We note that phase models with three- and higher-fold couplings are not unheard of, and in fact arise generically in certain contexts Ashwin and Rodrigues (2016).

Our first observation is that Eq.(1) admits a family of perfectly synchronized states, given by

[TABLE]

for any constant $\alpha$. In the context of the micro-resonator problem, these states correspond to a spatially localized pulse propagating around the micro-resonator ring, and $\alpha$ determines the location of the optical pulse in the co-moving frame. This definition of synchronization is analogous to the situation in the Kuramoto model, in which the oscillators converge to an arbitrary phase

[TABLE]

For the ternary model, the oscillators will have identical phases if $\alpha=0$, though there is nothing unique about that particular value.

It has been shown Taheri et al. (2017a) that this state is dynamically stable except for a single neutral eigenvalue corresponding to a shift in $\alpha$.

In what follows, we consider a generalization of Eq.(1), where we include some intrinsic quenched disorder:

[TABLE]

with the disorder constants $\omega_{j}$ drawn from a fixed, zero-mean distribution $g(\omega)$. We expect such imperfections on physical grounds. In the optical resonator context, the $\omega_{j}$ arise from cavity imperfections and/or high-order dispersion effects.

III Synchronization Behavior

We now consider the collective behavior of the ternary phase model, in particular the tendency to evolve into ordered states. Our approach is inspired by what is known about the Kuramoto model. To begin, we note that in the case of zero disorder ($\omega_{j}=0$), Eq.(4) is a gradient system $\dot{\Delta}_{j}=-\partial V/\partial\Delta_{j}$ with potential given by

[TABLE]

This structure implies the existence of locally stable perfectly synchronized states Gushchin, Mallada, and Tang (2015). Extensive numerical simulations suggest that, for almost all initial conditions, the disorder-free system evolves to a perfectly synchronized state. (In simulations of Eq.(1), a state satisfying condition (2) may not be immediately obvious, owing to the angular nature of the variables. A simple check can be performed by plotting the difference of adjacent $\Delta_{j}$ modulo $2\pi$.) In view of this, we expect the dynamics of Eq.(4) to exhibit a competition between the synchronizing influence of the coupling and the de-synchronizing influence of the disorder. The central issue is whether, and under what circumstances, the system exhibits an order-disorder transition, and if so what is the nature of the transition(s).

In order to investigate this question quantitatively, we introduce the following set of complex order parameters:

[TABLE]

$j=1,2,\dots,N$. At steady state, $|Z_{j}|$ is observed to be independent of index $j$, up to fluctuations of order $N^{-1/2}$. For a perfectly synchronized state (2), $|Z_{j}|$ is exactly one, while for randomly distributed $\Delta_{j}$ it tends to zero. In practice, it is useful to average $|Z_{j}|$ over all oscillators and consider $\overline{|Z|}$ to be the overall order parameter.

Our simulations show that, for long times, trajectories of Eq.(4) evolve linearly in time for all $\Delta_{j}$, i.e.

[TABLE]

The $\delta_{j}$ reflect deviation from (2) due to the added disorder. From Eq.(4), this relationship can only be satisfied if

[TABLE]

Taking a weighted sum of the $N$ equations yields, since the coupling terms cancel,

[TABLE]

so that

[TABLE]

and so

[TABLE]

For convenience, we can consider the co-rotating variables

[TABLE]

From (7) and (8), we can easily see that the transformed system $\Delta^{\prime}_{j}$ evolves to a stationary state. Incidentally, this transformation also ensures the existence of a constant of motion

[TABLE]

Unless otherwise noted, the simulations described below are performed in this frame by generating a set of $\omega_{j}$, calculating $\Omega$, and subtracting according to Eq.(12).

To investigate order-disorder transitions in our system, we performed simulations using three different distributions for the $\omega_{j}$: uniform random, Gaussian random, and a systematic deterministic disorder to be described later.

The simulations were carried out as follows. We started by generating a random set of $\omega_{j}$ with unit variance. Multiplying this set by a factor $\sigma$ gives a standard deviation $\sigma$ while retaining comparable statistical properties. Starting with $\sigma=0$ and random initial conditions, we numerically integrated Eq. (4) until the system reached its stationary synchronized state. We then increased $\sigma$ by a small amount and continued integrating, until a new steady state was reached (identified by an unchanging averaged order parameter $\overline{|Z|}$). Iterating this process, the system eventually reaches the unsynchronized regime at $\sigma_{1}$. At some distance into the unsynchronized regime we reversed the process by gradually reducing $\sigma$. Eventually, at $\sigma_{2}$, the system begins to synchronize again.

Figures 1a) and 1b)

show typical results for random Gaussian and uniform disorder, respectively. The results are qualitatively similar. For small, increasing $\sigma$, the order parameter gradually decreases. At the critical value $\sigma_{1}$, there is a first-order transition beyond which the order parameter becomes negligible. Then, decreasing $\sigma$ induces another first-order transition at $\sigma_{2}$. The two transitions occur at different values of $\sigma$, resulting in a hysteresis loop. As expected, the fluctuations grow near the transition points.

This behavior is different than what one sees in the Kuramoto model. For comparison, results for the latter are shown in Fig. 1c) (Gaussian disorder) and Fig. 1d) (uniform disorder). For uniformly distributed disorder, the Kuramoto model shows a first-order transition with no hysteresis Pazó (2005). In the case of Gaussian distributed disorder, it shows a second-order transition.

In addition to random disorder, we also considered disorder of the form

[TABLE]

with $c$ the effective “width” of the distribution. This choice is motivated by higher-order dispersion effects that arise in the micro-resonator problem Taheri et al. (2017a); Wen et al. (2016). The behavior in this case can be understood heuristically as follows. As $c$ is increased, the oscillators with large $j$, i.e. those far from the pump mode, acquire an intrinsic frequency much larger than the coupling constant and thus become completely unlocked from the rest of the population. The lower $j$ oscillators still have relatively low disorder, and can synchronize. The result is that the population is synchronized up to a critical $j$ above which the oscillators are effectively free-running. This is illustrated in Figure 2a)

by plotting $\dot{\Delta}_{j}$ and $\omega_{j}$ versus $j$ in the non-co-rotating frame. The locked oscillators fall on a straight line reflecting the synchronization condition (7-8) while the unlocked oscillators, beginning at $j\approx 200$, rotate at their intrinsic frequencies. While the locked oscillators can be thought of as being frequency locked, this is really a result of their phase locking; for the ternary model the two are necessarily connected through the synchronization condition.

Carrying out simulations using the same protocol as before, we find that the fraction of locked oscillators (those with $\dot{\Delta}_{j}$ falling on the straight line) decreases continuously to zero with no observed transition or hysteresis (Figure 2b). This suggests that cubic dispersion merely dilutes the total synchronization; for random noise, the first-order transitions enforce a critical threshold above which no synchronization is possible.

IV Thermodynamic Limit

In order to move beyond simulations, we formulated a thermodynamic limit corresponding to an infinite number of oscillators on a lattice with uniform spatial density. We show that, in this limit, there are no fluctuations and the collective state can be uniquely identified.

Unlike the all-to-all coupling of oscillators in the Kuramoto model, the ternary coupling in (4) is index-dependent. We therefore consider the oscillators to lie on a discrete lattice corresponding to $j=1,2,\dots,N$. Inspired by Ref. (Gupta, 2017), we consider a spatial continuum of oscillators which can be described by an explicitly position-dependent probability density.

First, we note that the equation of motion (4) can be written in terms of the order parameters (6) as

[TABLE]

In this form, the individual oscillators appear decoupled which allows for a self-consistent calculation of the order parameters. The order parameters $Z_{j}$ are written explicitly as

[TABLE]

reflecting the antisymmetry condition $\Delta_{-j}=-\Delta_{j}$. The $Z_{j}$ can be expressed in terms of a probability distribution rather than an explicit sum as

[TABLE]

where $P(\Delta,\Delta^{\prime},\omega,\omega^{\prime},x,x^{\prime})$ is the joint probability density of the oscillator at $x$ having phase $\Delta$ and disorder $\omega$ while the oscillator at $x^{\prime}$ has phase $\Delta^{\prime}$ and disorder $\omega^{\prime}$. In the spirit of Eq.(15), we consider the oscillators to be independent so that the probability density $P$ can be approximated as

[TABLE]

with $\rho$ the probability density for a single oscillator parameterized by position.

Because the coupling depends on oscillators at two different positions, the density $\rho$ appears twice. This is one of the major structural novelties of the model and significantly complicates stability analysis of the unsynchronized state.

In order to treat the infinite oscillator limit, in which case the probability densities vary continuously in space, we introduce a spatial coordinate $x$ such that

[TABLE]

From (19), $\textrm{d}j=(N-1)\textrm{d}x$, and the summation over $l$ transforms according to

[TABLE]

Taking $N\rightarrow\infty$, the sum becomes an integral and position $x$ becomes continuous, leading to

[TABLE]

The full thermodynamic system, completely described by $\rho$, is governed by the continuity equation

[TABLE]

and the constraints

[TABLE]

where $g(\omega)$ is the probability density for each oscillator’s disorder.

Having established the evolution equations for the thermodynamic limit, we may calculate the order parameter $R=\overline{|Z|}$ as a functional of the disorder distribution $g(\omega)$. The details of the calculation are relegated to the appendix, but the final result is

[TABLE]

valid for even, unimodal $g(\omega)$.

One sees immediately that $R=0$ is always a solution. To find non-zero solutions, Eq.(24) can be numerically solved for any given $g(\omega)$. The resulting predictions for Gaussian and uniform disorder are shown in Figure 3

alongside values from direct simulation for different $N$. As $N$ is increased, agreement along the upper branch improves considerably. Since $|Z|$ is non-negative, convergence along the lower incoherent branch is much slower. The middle solution branch is presumably unstable, as is the case for a saturating subcritical pitchfork bifurcation Strogatz (1994a).

Curiously, the unsynchronized state is stable for a wider range of $\sigma$ as $N$ is increased, and may even be stable for all values as $N\rightarrow\infty$. This would be consistent with the theoretically predicted curve in which the middle branch and $R=0$ branch meet exactly at the origin. Bistability of the synchronized and unsynchronized states would have important implications for the micro-resonator problem since even a perfect device could fail to achieve mode locking.

While our previous analysis hasn’t addressed stability, we note that linear stability of even the trivial $R=0$ state is complicated by the fact that the order parameter is quadratic in $\rho$. This can be understood from the following argument without a complete stability calculation. The incoherent $R=0$ state is represented by the uniform density

[TABLE]

Evaluation of the order parameter (21) using a perturbed solution of the form $\rho_{0}+\epsilon\delta\rho$ yields zero to linear order in $\epsilon$; one of the two factors of $\rho$ must be $\rho_{0}$, causing its associated exponential in $\Delta$ to vanish upon integration. However, the order parameter is responsible for coupling the oscillators as shown in (22). When it vanishes, $K$ disappears from the equations and so cannot be responsible for a bifurcation. As a result, a nonlinear stability analysis is required to describe any synchronizing transition from the incoherent state.

V Conclusion and Open Questions

In this paper, we investigated the synchronization behavior of a new phase model. We found that its transitions are characteristically different from those found in the classic Kuramoto model. We developed a thermodynamic, continuum limit which we used to calculate the order parameter for different disorder distributions. Although the predicted values for the order parameter agree well with simulations, the onset of instability for the unsynchronized solution is currently unexplained. With the enormous utility of the Ott-Antonsen ansatzOtt and Antonsen (2008) for the continuum limit of the Kuramoto model, it is natural to wonder if a similar dimensional reduction can be performed for the ternary model.

In addition to its novel dynamical properties, investigating the ternary phase model provides insight into the robustness of the synchronization process in optical microresonator cavities. Additional effects could be explored such as finite-size fluctuations and the inclusion of non-identical coupling constants $K_{j}$. The latter has been investigated in the absence of disorder Taheri et al. (2017a).

Perhaps the greatest outstanding challenge in studying the ternary phase model is the issue of stability. The difficulty of even linear stability analysis for the Kuramoto model is well knownStrogatz (2000). For the ternary model, linear analysis fails and a strictly nonlinear approach is apparently needed.

Appendix A Calculation of the Order Parameter

Two key observations from the discrete formulation allow us to proceed in calculating the thermodynamic limit of the order parameter for a given disorder distribution. First, the steady-state phase of the order parameter is linear in space, i.e. $\arg(Z_{j})=\alpha j$. Since $\alpha$ is an arbitrary constant, it can be set to zero. Second, $|Z_{j}|$ is observed to be independent of $j$ in steady state. We can therefore choose to evaluate $Z(x)$ arbitrarily at $x=0$ and set its phase to zero. Equation (15) then simplifies to

[TABLE]

where $R=|Z(0)|$.

In the spirit of Kuramoto’s solution to the Kuramoto modelStrogatz (1994b); Kuramoto (1984), we seek stationary solutions of (22) assuming a constant value of $R$. Oscillators with small disorder will have their phase entrained to the equilibria of (26) such that

[TABLE]

The remaining oscillators are unsynchronized, and stationarity requires

[TABLE]

with $a$ determined by normalization.

Combining Eqs. (21), (27), and (28) yields a self-consistency equation for the order parameter:

[TABLE]

For even, unimodal $g(\omega)$, (28) admits a symmetry in $\Delta$ and $\omega$ that causes the integrals over $\omega$ to vanish and so the unlocked oscillators do not contribute to the order parameter (see for example Ref.Strogatz, 1994b). Evaluating the $\Delta$ integrals, (29) reduces to

[TABLE]

Expanding the complex exponentials, the imaginary terms are linear in $\omega$ and $\omega^{\prime}$ respectively and thus vanish upon integration with even $g$. As a result, both integrands are even in $y$ and (30) reduces to

[TABLE]

Taking a change of variables

[TABLE]

we end up with the indicated result

[TABLE]