---
abstract: 'In this paper we study the synchronisation of three identical oscillators, i.e., clocks, hanging from the same hard support. We consider the case where each clock interacts with the other two clocks. The synchronisation is attained through the exchange of small impacts between each pair of oscillators. The fundamental result of this article is that the final locked state is at phase difference of $\frac{2\pi}{3}$ from successive clocks (clockwise or counter-clockwise). Moreover, the locked states attract a set whose closure is the global set of initial conditions. The methodology of our analysis consists in the construction a model, which is a non-linear discrete dynamical system, i.e. a non-linear difference equation. The results are extendable to any set of three oscillators under mutual symmetric interaction, despite the particular models of the oscillators.'
author:
- 'Emma D’Aniello and Henrique M. Oliveira'
- Dipartimento di Matematica e Fisica
- Università degli Studi della Campania Luigi Vanvitelli
- 'Viale Lincoln n. 5 - 81100 Caserta, Italia'
- |
    emma.daniello@unicampania.it\
    Department of Mathematics\
    Center for Mathematical Analysis, Geometry and Dynamical Systems,\
    Instituto Superior Técnico, University of Lisbon,\
    Av. Rovisco Pais, 1049-001, Lisboa, Portugal.\
    holiv@math.tecnico.ulisboa.pt
title: Huygens synchronisation of three clocks equidistant from each other
---

Introduction
============

Synchronisation among oscillators with some form of coupling has been called universal [@strogatz2004] and is prevalent concept in Nature [@Pit].

In 1665 Huygens, the inventor of the pendulum clock, observed synchronisation between two pendulum clocks [@Huy] hanging from the same support. The first observation of Huygens was about the two clocks hanging in the same wall beam of his house when he was lying in bed with some indisposition. He observed both phase and phase opposition as a final state of the coupled system. The second observation was made later, when Huygens hung the two clocks on a board sitting on two chairs.

The two systems observed by Huygens are quite different, therefore originating two lines of research completely separated. The later case has been studied in many papers [@Benn; @Col; @Frad; @Jova; @Col2; @Martens; @Oud; @Sen] by considering momentum conservation in the clocks-beam system, since the plank supporting the clocks is able to move. The system in that case is non-perturbative, there are three bodies that can move, i.e., three degrees of freedom: the pendula and the wood beam. The model is a classical mechanics paradigm, and in most of the works cited above the friction is considered viscous and not dry.

The first case, when the pendulum clocks are suspended at a very rigid house beam, therefore not able to move, has been approached in the works [@Abr2; @Abr; @Vass]. In this model the interaction is considered perturbative. When one of the pendulums suffers the internal impact from the clock escape mechanism, a small travelling wave perturbs the second one and vice-versa.

Recently, in [@OlMe] a theoretical model for this interaction was developed. In the same work, simulation and some experimental studies were carried on, that confirmed the validity of the proposed model when the support wall has infinite mass and, therefore, not moving. Naturally, the center of mass of the wall is not, in this case, a degree of freedom of the system and the coupling is, via very weak travelling waves, propagated in the rigid structure of the wall.

In this article, we use as a conceptual starting point [@OlMe] presenting a mathematical model where the coupling is assumed to be attained through the exchange of impacts between three identical oscillators, where each one of the clocks interacts with the two other clocks. The model presents the advantage of being independent of the physical nature of the oscillators, and thus it can be used in other oscillator systems where synchronisation and phase locking has been observed [@Pit].

The ideas for the model presented here originated from the Andronov [@And; @OlMe] model of the phase-space limit cycle of isolated clocks, and assumes the exchange of single impacts (travelling solitons, for this system) between the oscillators at a specific point of the limit cycle.

The fundamental hypotheses in this article are the existence of an asymptotically stable limit cycle for each oscillator and one very small interaction between each pair of clocks per cycle.

We point out that in [@OlMe] the authors obtained phase opposition, which is in line with the original Huygens observations [@Huy]. In this paper we obtain a particularly symmetric asymptotic state at which all the clocks remain at a phase difference of $\frac{2\pi}{3}$ between each other.

A natural step forward is to generalise the results obtained here to a larger set of oscillators and apply these results to bidimensional and tridimensional swarms of oscillators. The results of the present work can be used, namely, to study interacting insects or neuronal networks, that have been studied using the over-simplified integrate and fire models of Kuramoto [@Campbell1997; @Campbell1999; @Kuramoto1975; @Mirollo1990; @Strogatz2000].

The paper is organised in $5$ sections. In section $2$, we discuss the original model of the pendulum clock and we briefly recall the model for two identical clocks. In section $3$, we deduce the model for three identical clocks hanging at the same wall with mutual interactions. In section$\ 4$, we analyse the model, computing its symmetries and stabilities. In section $5 $, we draw conclusions and point out directions for future work.

Model for the synchronisation of two oscillators
================================================

Some background
---------------

For the sake of completeness we present here a very short theory of synchronisation for two oscillators exchanging small perturbations at each cycle. We consider identical oscillators. This theory can be applied to networks of identical oscillators, electronic oscillators and many other real world systems. We intend to consider, in future work, the case of slighly different oscillators, which give rise to regions of stability versus instability in the parameter space of these systems, i.e., Arnold Tongues [@boyland1986; @gilmore2011].

For basic, classical definitions and notions related to synchronisation, like phase and frequence, we follow and refer to [@Pit], and for concepts concerning general theory of dynamical systems like, for instance, limit cycle, we refer to [@arrowsmith1990]. In this paper, we always assume that an oscillator is a dynamical system having a limit cycle. We use the word clock when referring to a special type of oscillator described by the Andronov model [@OlMe].

Given a point $p_{0}$ in the limit cycle $\gamma$, the necessary time to return to $p_{0}$, after one round on the limit cycle, is the period $T_{0}$. A phase $\varphi$ is a real coordinate that describes the position of the representative point of the system on the limit cycle [@Nakao; @Pit].

Let $B_{\gamma}$ be the basin of attraction of the limit cycle. Consider the points outside the limit-cycle $\gamma$ but in $B_{\gamma}$. We extend the definition of phase to $B_{\gamma}$ as follows. We assign the same phase $\varphi$ to all points $p$ in $B_{\gamma}$ that converge to the same $p_{0}$ on the limit cycle $\gamma$ as $t\rightarrow\infty$, being the phase of $p_{0}$ precisely $\varphi$ [@gu1975]. The set of points $p$ that share the same phase is an *isochron curve*. If the oscillator’s states are on the same isochron at a given point in time, they continue to be on the same isochron in time [@gu1975; @Nakao]. When each clock suffers a perturbation, its state can go slightly off the limit cycle and generically jump to another isochron. Moreover, we assume that the limit cycles are structurally stable under small perturbations.

When we consider two oscillators, $1$ and $2$, with orbits on the limit-cycle or sufficiently near the limit cycle, each one has a particular phase, respectively $\varphi$ and $\psi$.

The study of the synchronisation of these two oscillators consists in establishing a dynamical system for the phase difference of the two oscillators.

We have two possible lines of research [@Pit]. The first is to consider the phase difference along continuous time, i.e., to look at the function $\phi\left(  t\right)  =\psi\left(  t\right)  -\varphi\left(  t\right)  $ for $t\in\left[  0,+\infty\right[  $. The second line of research, that we adopt in this paper, is to consider the phase difference $\phi_{n}=\psi_{n}%
-\varphi_{n}$ taken at discrete instants $n=0,1,2,\ldots$. In this paper, we consider exclusively this last approach.

There is phase synchronisation when the phase differences between the oscillators tend to a specific attractor. When this attractor is an isolated point, then there is phase locking. Naturally, richer coupled states can occur [@Martens]. The main goal for any theory of synchronisation is to obtain this phase difference dynamics and to establish the existence and nature of the attractor. In the case of Huygens observations, the attractor was the point $0$ or the point $\pi$ and the phase dynamics was unidimensional.

The Andronov model for an isolated clock
----------------------------------------

We recall here the model for the sake of completeness of this article. Assuming that dry friction predominates in the internal metal pieces of the clock and the viscous damping is not predominant, using the angular coordinate $q$, the differential equation governing the isolated pendulum clock is$$\ddot{q}+\mu\text{ }\operatorname*{sign}\dot{q}+q=0,$$ where $\mu>0$ is the dry friction coefficient and $\operatorname*{sign}\left(
x\right)  $ is the classical function taking the value $-1$ at $x<0$ and $1$ at $x>0$. In [@And] it was considered that, in each cycle, the escape mechanism gives to the pendulum a fixed amount of normalized kinetic energy $\frac{h^{2}}{2}$ so to compensate the loss of kinetic energy occurred because of the dry friction in each complete cycle. This transfer of kinetic energy is called a *kick*. The origin is fixed so that the kick is given precisely when $q=-\mu$. The phase portrait is shown in Fig. $1$.

\[h\]

[limitcycle.eps]{}

As in [@OlMe], with initial conditions $q\left(  t=0\right)  =-\mu$ and $\dot{q}\left(  t=0\right)  =v_{0}$, a Poincaré section [@Bir] is represented vol. II, page 268) as the half line $q=-\mu^{+}$ and $\dot{q}>0$ [@And]. The symbol $+$ means that we are considering that the section is taken immediately after the kick. Due to friction during a complete cycle, a loss of velocity $-4\mu$ occurs. By considering the velocity, $v_{n}=\dot
{q}\left(  2n\pi^{+}\right)  $, at the Poincaré section in each cycle, the non-linear discrete dynamical system [@And] is obtained $$v_{n+1}=\sqrt{\left(  v_{n}-4\mu\right)  ^{2}+h^{2}}\text{.} \label{recu}%$$ This equation has the asymptotically stable fixed point $$v_{f}=\frac{h^{2}}{8\mu}+2\mu\text{.}%$$ Any initial condition $v_{0}\in\left(  4\mu,+\infty\right)  $ is attracted to $v_{f}$. Each cycle corresponds to a phase increment of $2\pi$ and the phase $\varphi$ is linear with respect to $t$, precisely $$\varphi=2\pi t.$$ As already mentioned, the nature of limit cycle is not of fundamental importance when we consider the interaction of three identical clocks, as we shall see in the sequel. We have presented here the basis of our reasonings in the non-usual case when the computations of the limit cycle are explicit and the usual angular phase is a linear function of $t$.

Two interacting oscillators
---------------------------

We present briefly the conclusions of considering two pendulum clocks suspended at the same wall, in a simplified version of [@OlMe], where the clocks are assumed to have natural angular frequencies near each other but different. Here, we assume that the two clocks have the same angular frequency. When one clock receives the kick, the impact propagates in the wall slightly perturbing the second clock. The perturbation is assumed to be instantaneous since the time of travel of sound in the wall between the clocks is assumed very small compared to the period.

Consider two oscillators and index them by $i=1,2$. Each oscillator satisfies the differential equation $$\ddot{q}_{i}+\mu_{i}\text{ sign\ }\dot{q}_{i}+q_{i}=-\alpha_{i}\digamma\left(
q_{j}\right)  ,\text{ for }i,j=1,2\text{, }i\not =j\text{.}
\label{coupledandro}%$$ As in the *Andronov model*, the kinetic energy of each oscillator increases of the fixed amount $h_{i}$ when $q_{i}=-\mu_{i}$. The coupling term is the normalised force $-\alpha_{i}\digamma\left(  q_{j}\right)  $, where $\digamma$ is the interaction function and $\alpha_{i}$ a constant with acceleration dimension. Following [@OlMe], the effect of the interaction function $\digamma$ is considered to produce an increment $-\alpha$ in the velocity of each clock, leaving the position invariant when the other is struck by the energy kick. The reader finds the detailed treatment in [@OlMe]. Here we only recall some ideas from that article, for the sake of completeness and to make our three clocks model more simple and natural to deal with.

To describe and investigate the effect of the kicks, we construct a discrete dynamical system for the phase difference between the two clocks. We compute each cycle using as reference one of the clocks (the choice is irrelevant, since the model is symmetric). We choose, to fix ideas, clock 1 as the reference: whenever its phase reaches $0$ $\left(  \operatorname{mod}%
2\pi\right)  $, the number of cycles increases one unit from $n$ to $n+1$.

If there exists an attracting fixed point for that dynamical system, the phase locking occurs. As in [@OlMe], the assumptions are the following.

1.  Dry friction.

2.  The pendulums have the same natural angular frequency $\omega=1$.

3.  The perturbation in the momentum is always in the same vertical direction in the phase space [@Abr2; @Abr].

4.  Since the clocks have the same construction, the energy dissipated at each cycle of the two clocks is the same, $h_{1}=h_{2}=h$. The friction coefficient is the same for both clocks, $\mu_{1}=\mu_{2}=\mu$.

5.  The perturbative interaction is instantaneous. This is a reasonable assumption, since in general the perturbation propagation time between the two clocks is several orders of magnitude lower than the periods.

6.  The interaction is symmetric, the coupling has the same, very small, constant $\alpha$ when the clock $1$ acts on clock $2$, and conversely.

We compute at this point the phase difference when clock $1$ returns to the initial position. The secular repetition of perturbations leads the system with the two clocks in phase opposition as Huygens observed in 1665 [@Huy]. The discrete dynamical model that we deduce from [@OlMe] for the phase difference between the two clocks $\phi_{n}=\psi_{n}-\varphi_{n}$ is the Adler equation [@adler1946study; @Pit]$$\phi_{n+1}=\phi_{n}+\varepsilon\sin\phi_{n}, \label{Perturb}%$$ with a very small constant $\varepsilon=\frac{16\mu\alpha}{h^{2}}$. In the interval $\left[  0,2\pi\right[  $, there are two fixed points which are $\pi$ and $0$ respectively attracting and reppeling.

Equation (\[Perturb\]) is the starting point from where we begin, in the present paper, the study the three symmetric clocks in mutual interaction.

In any model with a perturbation of phase given by equation (\[Perturb\]) per cycle, i.e., Adler’s perturbation [@adler1946study; @Pit], despite being a physical clock (with Andronov model or any different model) or other type of oscillator, electric, quantic, electronic or biological, the theory presented here for three oscillators interacting by small impacts will be exactly the same, with the same conclusions.

Model for three pendulum clocks placed in the three vertices of an equilateral triangle
=======================================================================================

Hypotheses
----------

We consider three pendulum clocks suspended at the same wall, placed in the three vertices of an equilateral triangle, say the vertices are $A$, $B$, and $C$ and $B$ are the extreme points of the basis of the triangles.

\[ptb\]

[triangle.eps]{}

This geometric setting is purely conceptual. Any set of three dynamical systems receiving symmetric impacts from the other two will have the same type of response of the clocks depicted in the three vertices of an equilateral triangle.

Call the clocks placed in the three vertices $A$, $B$ and $C$, respectively, $O_{1}$, $O_{2}$ and $O_{3}$. When the clock $A$ receives the kick from the escape mechanism, the impact propagates in the wall slightly perturbing the other two clocks. As in [@OlMe], the perturbation is assumed to be instantaneous, since the time of travel of sound in the wall between the clocks is assumed very small compared to the period. As for the two clocks model discussed in [@OlMe], we make the following assumptions, now formulated for three clocks.

1.  The system has [dry friction ]{}[@And].

2.  \[work3\][The pendulums of clocks ]{}$O_{1},$ $O_{2}$ and $O_{3}$ [have respectively natural angular frequencies ${\omega}_{1}={\omega}%
    _{2}={\omega}_{3}=1$.]{}

3.  [The perturbation in the momentum is always in the same vertical direction in the phase space ]{} [@Abr2; @Abr].

4.  [The friction coefficient is the same for all the three clocks, ${\mu
    }_{1}={\mu}_{2}={\mu}_{3}={\mu}$. The energy dissipated at each cycle of the three clocks is the same, and the energy furnished by the escape mechanism to compensate the loss of energy to friction in each cycle is $h_{1}=h_{2}%
    =h_{3}=h$. ]{}

5.  [The perturbative interaction is instantaneous. This is a reasonable assumption, since in general the perturbation propagation time between two clocks is several orders of magnitude lower than the periods [@OlMe].]{}

6.  \[hyp6\][The interaction is symmetric. The couplings have the same constant $\alpha$ when one clock acts on another and conversely. In this model $\alpha$ is assumed to be very small.]{}

7.  \[hyp67\][Each perturbation from clock $i$ to clock $j$ (where $i,j\in\{1,2,3\}$ with $i\not =j$), when clock $i$ suffers its internal impact of kinetic energy $h^{2}$, gives rise to a small perturbative change of phase which is in first order a $2\pi$-periodic differentiable odd function $P$ of the real variable $\phi$$$P\left(  \phi\right)  =\varepsilon\sin\phi\text{,} \label{Perturb1}%$$ where $\phi={\phi}_{ij}$ is the phase difference between clock $i$ and clock $j$.]{}

The value of $\varepsilon$ is the above mentioned $\varepsilon=\frac
{8\mu\alpha}{h^{2}}$ from [[@OlMe] where $\mu$ is the dry friction coefficient, $\frac{h^{2}}{2}$ is the kinetic energy furnished by the internal escape mechanism of each clock once per cycle and $\alpha$ the interaction coefficient between the clocks. The greater the $\alpha$ is, the greater the mutual influence among the clocks. In this paper, we do not need to particularize ]{}$\varepsilon$, since we are not interested in doing experimental computations. Therefore, we are interested in the fundamental result of symmetry between three oscillators subject to very weak mutual symmetric interaction.

Most of the reasonings are independent on the form of the function $P\left(
\phi\right)  $, therefore we consider a general differentiable odd function of the real variable $\phi$, $P(\phi)$, for the development of the model, and consider it of the form (\[Perturb1\]) when we analyze the model in section 4.

Observe that $|\sin(x+\varepsilon\sin y)-\sin x|<{\varepsilon}$ when $\varepsilon$[ is assumed to be sufficiently small. Therefore, we restrict our model to first order. ]{}We consider all the values of variables and constants in IS units.

Construction of the model
-------------------------

We now construct a dynamical system using as reference the phase of the clock in the vertex A (= clock $O_{1}$). This reference is arbitrary: any of the clocks can be used as the reference clock with the same results at the end, since the system is symmetric. We compute the effects of all phase differences and perturbations when the clock at A makes a complete cycle returning to the initial position. Without loss of generality, we consider the next working hypotheses.

1.  \[work1\]The initial phase of clock at A at $t=0^{-}$ is zero, i.e., $\psi_{1}(0^{-})=0^{-}$, the minus ($-$) superscript means that at the instant $0^{-}$ clock $1$ is just about to receive the internal energy kick from its escape mechanism.

2.  \[work2\]We consider that the initial phases of the three clocks are: $\psi_{3}(0^{-})=\psi_{3}^{0}>\psi_{2}(0^{-})=\psi_{2}^{0}>0^{-}=\psi
    _{1}(0^{-})=\psi_{1}^{0}$.

3.  \[work4\]The perturbation satisfies the relation $P\left(
    x+Px\right)  \simeq Px$ in first order.

To obtain the desired model, we need to proceed through 6 steps, starting from the following initial conditions, that is the phase differences of all pairs of clocks.

In the sequel ${\psi}_{i}^{j}$ denotes the phase of clock $O_{i}$ at the $j-th$ step.

**INITIAL CONDITIONS**

The phase difference between $O_{3}$ and $O_{1}$ is $$(CA)_{0}={\psi}_{3}^{0}-{\psi}_{1}^{0}={\psi}_{3}^{0},$$ and the phase difference between $O_{1}$ and $O_{3}$ is symmetric, in the sense that $$(AC)_{0}={\psi}_{1}^{0}-{\psi}_{3}^{0}=-{\psi}_{3}^{0}=-(CA)_{0}.$$ The phase difference between $O_{2}$ and $O_{1}$ is $$(BA)_{0}={\psi}_{2}^{0}-{\psi}_{1}^{0}={\psi}_{2}^{0}$$ and the phase difference between $O_{1}$ and $O_{2}$ is $$(AB)_{0}={\psi}_{1}^{0}-{\psi}_{2}^{0}=-{\psi}_{2}^{0}=-(BA)_{0}.$$ The phase difference between $O_{3}$ and $O_{2}$ is $$(CB)_{0}={\psi}_{3}^{0}-{\psi}_{2}^{0}$$ and the phase difference between $O_{2}$ and $O_{3}$ is $$(BC)_{0}={\psi}_{2}^{0}-{\psi}_{3}^{0}=-(CB)_{0}.$$

**STEPS LEADING TO THE CONSTRUCTION OF THE MODEL**

**STEP 1:** first impact. Interactions of $O_{1}$ on $O_{2}$ and of $O_{1}$ on $O_{3}$, at $t=0$.

When the system in position A attains phase $0$ $(\operatorname{mod}2\pi)$ it receives a sudden supply of energy, for short a kick, from its escape mechanism, this kick propagates in the common support of the three clocks and reaches the other two clocks.

Now, the phase difference between $O_{3}$ and $O_{1}$ is corrected by the perturbative value $P$: $$\left(  CA\right)  _{I}=\left(  CA\right)  _{0}+P\left(  \left(  CA\right)
_{0}\right)  ={\psi}_{3}^{0}+P\left(  {\psi}_{3}^{0}\right)  .$$ The phase difference between $O_{1}$ and $O_{3}$ is $$\left(  AC\right)  _{I}=\left(  AC\right)  _{0}+P\left(  \left(  AC\right)
_{0}\right)  =-{\psi}_{3}^{0}+P\left(  -{\psi}_{3}^{0}\right)  =-(CA)_{I},$$ since $P$ must be an odd function of the mutual phase difference.

The phase difference between $O_{2}$ and $O_{1}$ is $$\left(  BA\right)  _{I}=\left(  BA\right)  _{0}+P\left(  \left(  BA\right)
_{0}\right)  ={\psi}_{2}^{0}+P\left(  {\psi}_{2}^{0}\right)  ,$$ and the symmetric phase difference between $O_{1}$ and $O_{2}$ is $$\left(  AB\right)  _{I}=\left(  AB\right)  _{0}+P\left(  \left(  AB\right)
_{0}\right)  =-{\psi}_{2}^{0}+P\left(  -{\psi}_{2}^{0}\right)  )=-\left(
BA\right)  _{I}.$$

The phase difference between $O_{3}$ and $O_{2}$ depends on $\left(
CA\right)  _{I}$ and $\left(  BA\right)  _{I}$ and it is $$\left(  CB\right)  _{I}=\left(  CA\right)  _{I}-\left(  BA\right)  _{I}={\psi
}_{3}^{0}-{\psi}_{2}^{0}+P({\psi}_{3}^{0})-P({\psi}_{2}^{0})=-(CA)_{I},$$

**STEP 2**: first natural time shift. The next clock to arrive at $2\pi^{-}$, from working hypothesis 3.2 (\[work2\]), is the clock ${O}_{3}$ at vertex $C$. The situation right before $O_{3}$ receives its kick of energy is when the phase of this clock is $2\pi^{-}$.

At this point we have$$%
\begin{cases}
\psi_{3}^{2} & ={2\pi}^{-}\\
\psi_{1}^{2} & =2\pi-(CA)_{I}=2\pi+(AC)_{I}=2\pi-\left(  {\psi}_{3}%
^{0}+P({\psi}_{3}^{0})\right) \\
\psi_{2}^{2} & =2\pi-\left(  CB\right)  _{I}=2\pi+\left(  BC\right)  _{I}%
=2\pi+{\psi}_{2}^{0}-{\psi}_{3}^{0}+P({\psi}_{2}^{0})-P({\psi}_{3}^{0}).
\end{cases}$$

**STEP 3**: second impact. Clock $O_{3}$ receives its internal kick, at the position $2\pi$.

Now, we have$$%
\begin{cases}
\psi_{3}^{3} & ={2\pi}\\
\psi_{1}^{3} & =\psi_{1}^{2}+P(\psi_{1}^{2})\\
& =2\pi-\left(  {\psi}_{3}^{0}+P({\psi}_{3}^{0})\right)  +P\left(
2\pi-\left(  {\psi}_{3}^{0}+P({\psi}_{3}^{0})\right)  \right) \\
& =2\pi-\left(  {\psi}_{3}^{0}+P({\psi}_{3}^{0})\right)  -P\left(  {\psi}%
_{3}^{0}+P({\psi}_{3}^{0})\right) \\
& \simeq2\pi-{\psi}_{3}^{0}-2P\left(  {\psi}_{3}^{0}\right) \\
\psi_{2}^{3} & =\psi_{2}^{2}+P(\psi_{2}^{2})\\
& =2\pi+{\psi}_{2}^{0}-{\psi}_{3}^{0}+P({\psi}_{2}^{0})-P({\psi}_{3}^{0})\\
& +P(2\pi+{\psi}_{2}^{0}-{\psi}_{3}^{0}+P({\psi}_{2}^{0})-P({\psi}_{3}^{0}))\\
& =2\pi+{\psi}_{2}^{0}-{\psi}_{3}^{0}+P({\psi}_{2}^{0})-P({\psi}_{3}^{0})\\
& +P({\psi}_{2}^{0}-{\psi}_{3}^{0}+P({\psi}_{2}^{0})-P({\psi}_{3}^{0}))\\
& \simeq2\pi+{\psi}_{2}^{0}-{\psi}_{3}^{0}+P\left(  {\psi}_{2}^{0}\right)
-P\left(  {\psi}_{3}^{0}\right)  +P\left(  {\psi}_{2}^{0}-{\psi}_{3}%
^{0}\right) \\
&
\end{cases}$$

**STEP 4**: second natural time shift. The next clock to arrive at $2\pi^{-}$, from working hypothesis3.2 (\[work2\]), is the clock $O_{2}$ at vertex $B$. The situation right before $O_{2}$ receives its kick of energy is when the phase of this clock is $2\pi^{-}$.

Then we have $$%
\begin{cases}
\psi_{2}^{4} & =2\pi^{-}\\
\psi_{1}^{4} & =\psi_{1}^{3}+2\pi-\psi_{2}^{3}\\
& \simeq2\pi-{\psi}_{3}^{0}-2P\left(  {\psi}_{3}^{0}\right)  +2\pi\\
& -\left(  2\pi+{\psi}_{2}^{0}-{\psi}_{3}^{0}+P\left(  {\psi}_{2}^{0}\right)
-P\left(  {\psi}_{3}^{0}\right)  +P\left(  {\psi}_{2}^{0}-{\psi}_{3}%
^{0}\right)  \right) \\
& =2\pi-{\psi}_{2}^{0}-P\left(  {\psi}_{2}^{0}\right)  -P\left(  {\psi}%
_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}^{0}\right) \\
\psi_{3}^{4} & =\psi_{3}^{3}+2\pi-\psi_{2}^{3}\\
& \simeq2\pi+{2\pi}-\left(  2\pi+{\psi}_{2}^{0}-{\psi}_{3}^{0}+P\left(  {\psi
}_{2}^{0}\right)  -P\left(  {\psi}_{3}^{0}\right)  +P\left(  {\psi}_{2}%
^{0}-{\psi}_{3}^{0}\right)  \right) \\
& \simeq2\pi-{\psi}_{2}^{0}+{\psi}_{3}^{0}-P\left(  {\psi}_{2}^{0}\right)
+P\left(  {\psi}_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}%
^{0}\right)  .
\end{cases}$$

**STEP 5**: third impact. Clock $O_{2}$ receives its internal energy kick. It reaches the position $2\pi$.

Then we have $$%
\begin{cases}
\psi_{2}^{5} & ={2\pi}\\
\psi_{3}^{5} & =\psi_{3}^{4} +P(\psi_{3}^{4})\\
& \simeq2\pi-{\psi}_{2}^{0}+{\psi}_{3}^{0}-P\left(  {\psi}_{2}^{0}\right)
+P\left(  {\psi}_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}%
^{0}\right) \\
& +P(2\pi-{\psi}_{2}^{0}+{\psi}_{3}^{0}-P\left(  {\psi}_{2}^{0}\right)
+P\left(  {\psi}_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}%
^{0}\right)  )\\
& \simeq2\pi-{\psi}_{2}^{0}+{\psi}_{3}^{0}-P\left(  {\psi}_{2}^{0}\right)
+P\left(  {\psi}_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}%
^{0}\right)  -P({\psi}_{2}^{0}-{\psi}_{3}^{0})\\
& =2\pi-{\psi}_{2}^{0}+{\psi}_{3}^{0}-P\left(  {\psi}_{2}^{0}\right)
+P\left(  {\psi}_{3}^{0}\right)  -2P\left(  {\psi}_{2}^{0}-{\psi}_{3}%
^{0}\right) \\
\psi_{1}^{5} & =\psi_{1}^{4}+P\left(  \psi_{1}^{4} \right) \\
& \simeq2\pi-{\psi}_{2}^{0}-P\left(  {\psi}_{2}^{0}\right)  -P\left(  {\psi
}_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}^{0}\right)  +\\
& P\left(  2\pi-{\psi}_{2}^{0}-P\left(  {\psi}_{2}^{0}\right)  -P\left(
{\psi}_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}^{0}\right)  \right)
\\
& \simeq2\pi-{\psi}_{2}^{0}-P\left(  {\psi}_{2}^{0}\right)  -P\left(  {\psi
}_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}^{0}\right)  -P({\psi
}_{2}^{0})\\
& =2\pi-{\psi}_{2}^{0}-2P\left(  {\psi}_{2}^{0}\right)  -P\left(  {\psi}%
_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}^{0}\right)  .
\end{cases}$$

**STEP 6 (the final)**: third natural time shift. The next clock to arrive at $2\pi^{-}$, from working hypothesis 3.2 (\[work2\]), is the clock $O_{1}$ at vertex $A$. The situation before $O_{1}$ receives its kick of energy is when the phase of this clock is $2\pi^{-}$, i.e., the cycles is complete.

At this point we are able to describe what happens to the phases after a complete cycle of the reference clock.

We have $$%
\begin{cases}
\psi_{1}^{6} & ={2\pi}^{-}\\
\psi_{2}^{6} & =\psi_{2}^{5}+2\pi-\psi_{1}^{5}\\
& \simeq2\pi+2\pi-\left(  2\pi-{\psi}_{2}^{0}-2P\left(  {\psi}_{2}^{0}\right)
-P\left(  {\psi}_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}%
^{0}\right)  \right) \\
& =2\pi+{\psi}_{2}^{0}+2P\left(  {\psi}_{2}^{0}\right)  +P\left(  {\psi}%
_{3}^{0}\right)  +P\left(  {\psi}_{2}^{0}-{\psi}_{3}^{0}\right)  ;\\
\psi_{3}^{6} & =\psi_{3}^{5}+2\pi-\psi_{1}^{5}\\
& \simeq2\pi-{\psi}_{2}^{0}+{\psi}_{3}^{0}-P\left(  {\psi}_{2}^{0}\right)
+P\left(  {\psi}_{3}^{0}\right)  -2P\left(  {\psi}_{2}^{0}-{\psi}_{3}%
^{0}\right)  +2\pi\\
& -\left(  2\pi-{\psi}_{2}^{0}-2P\left(  {\psi}_{2}^{0}\right)  -P\left(
{\psi}_{3}^{0}\right)  -P\left(  {\psi}_{2}^{0}-{\psi}_{3}^{0}\right)  \right)
\\
& =2\pi+{\psi}_{3}^{0}+P({\psi}_{2}^{0})+2P({\psi}_{3}^{0})-P({\psi}_{2}%
^{0}-{\psi}_{3}^{0}).
\end{cases}$$

Now, we compute the phase differences after the first cycle of $O_{1}$.

We have $$\begin{aligned}
(BA)_{I}  &  =-(AB)_{I}=\psi_{2}^{6}-\psi_{1}^{6}\\
&  \simeq2\pi+{\psi}_{2}^{0}+2P\left(  {\psi}_{2}^{0}\right)  +P\left(  {\psi
}_{3}^{0}\right)  +P\left(  {\psi}_{2}^{0}-{\psi}_{3}^{0}\right)  -2\pi\\
&  ={\psi}_{2}^{0}+2P\left(  {\psi}_{2}^{0}\right)  +P\left(  {\psi}_{3}%
^{0}\right)  +P\left(  {\psi}_{2}^{0}-{\psi}_{3}^{0}\right) \\
&  =(BA)_{0}+2P((BA)_{0})+P((CA)_{0})+P((BA)_{0}-(CA)_{0})\\
&\end{aligned}$$ and$$\begin{aligned}
&  (CA)_{I}\\
&  =-(AC)_{I}=\psi_{3}^{6}-\psi_{1}^{6}\\
&  =2\pi+{\psi}_{3}^{0}+P({\psi}_{2}^{0})+2P({\psi}_{3}^{0})-P({\psi}_{2}%
^{0}-{\psi}_{3}^{0})-2\pi\\
&  ={\psi}_{3}^{0}+P({\psi}_{2}^{0})+2P({\psi}_{3}^{0})-P({\psi}_{2}^{0}%
-{\psi}_{3}^{0})\\
&  =({(CA)}_{0})+P({(BA)}_{0})+2P({(CA)}_{0})-P((BA)_{0}-(CA)_{0})\\
&\end{aligned}$$ Hence, if we set $x=BA$ and $y=CA$, we obtain the system $$\left\{
\begin{array}
[c]{c}%
x_{1}=x_{0}+2P(x_{0})+P(y_{0})+P(x_{0}-y_{0})\\
y_{1}=x_{0}+P(x_{0})+2P({y}_{0})-P(x_{0}-y_{0}).
\end{array}
\right.$$

**THE MODEL**

By iterating the argument above, we get, for $n$ equal to the number of cycles described by $O_{1}$, the discrete dynamical system:$$\left\{
\begin{array}
[c]{c}%
x_{n+1}=x_{n}+2P(x_{n})+P(y_{n})+P(x_{n}-y_{n})\\
y_{n+1}=y_{n}+P(x_{n})+2P({y}_{n})-P(x_{n}-y_{n}).
\end{array}
\right.$$ If we write $$\left\{
\begin{array}
[c]{c}%
\varepsilon\varphi\left(  x,y\right)  =2P(x)+P(y)+P(x-y)\\
\varepsilon\gamma\left(  x,y\right)  =P(x)+2P({y})+P(y-x),
\end{array}
\right.$$ then we have $$\varphi\left(  x,y\right)  =\gamma\left(  y,x\right)  \text{,}%$$ and the iteration is a perturbation of the identity as $$\left[
\begin{array}
[c]{c}%
x_{n+1}\\
y_{n+1}%
\end{array}
\right]  =\left[
\begin{array}
[c]{cc}%
1 & 0\\
0 & 1
\end{array}
\right]  \left[
\begin{array}
[c]{c}%
x_{n}\\
y_{n}%
\end{array}
\right]  +\varepsilon\left[
\begin{array}
[c]{c}%
\varphi(x_{n},y_{n})\\
\varphi(y_{n},x_{n})
\end{array}
\right]  ,$$ that we can also write as $$X_{n+1}=F(X_{n})=X_{n}+\varepsilon\Omega(X_{n}), \label{Model}%$$ where $$X_{n+1}=\left[
\begin{array}
[c]{c}%
x_{n+1}\\
y_{n+1}%
\end{array}
\right]  ,$$ $$%
\begin{array}
[c]{c}%
F(X_{n})
\end{array}
=\left[
\begin{array}
[c]{cc}%
1 & 0\\
0 & 1
\end{array}
\right]  \left[
\begin{array}
[c]{c}%
x_{n}\\
y_{n},
\end{array}
\right]$$ and $$%
\begin{array}
[c]{c}%
\Omega(X_{n})
\end{array}
=\left[
\begin{array}
[c]{c}%
\varphi(x_{n},y_{n})\\
\varphi(y_{n},x_{n})
\end{array}
\right]  .$$

We now consider $P\left(  x\right)  =\varepsilon\sin x$, where $\varepsilon
=\frac{\alpha\mu}{8h^{2}}$ from hypothesis \[Perturb1\], explicitly,$$\begin{aligned}
\varphi\left(  x,y\right)   &  =2\sin x+\sin y+\sin\left(  x-y\right) \\
\gamma\left(  x,y\right)   &  =\sin x+2\sin y+\sin\left(  y-x\right)  .\\
&\end{aligned}$$

Analysis of the model
=====================

Fixed points and local stability
--------------------------------

In this section, we analyze the model (\[Model\]) obtained in the previous section. In a nutshell, in this section, we see that the system is differentiable and invertible in $S=\left[  0,2\pi\right]  \times\left[
0,2\pi\right]  $ when $\varepsilon>0$ is small. The perturbation map $\Omega\left(  x,y\right)  $ is periodic in $%
\mathbb{R}
^{2}$. This implies that the solution of the problem in the set $S$ is a dynamical system and not the usual semi-dynamical system associated with discrete time. That will provide a reasonable simple structure to the problem of the stability of fixed points and will enable to derive global properties. Moreover, we prove that for small $\varepsilon$ the set $S$ is invariant for the dynamics of $F$, meaning that the two phase differences of oscillators $O_{2}$ and $O_{3}$ relative to oscillator $O_{1}$ stay in the interval $\left[  0,2\pi\right[  $.

In particular, the map $\Omega$ has the zeros $\left(  \pi,\pi\right)  $, $\left(  \frac{2}{3} \pi,\frac{4}{3} \pi\right)  $ and $\left(  \frac{4}{3}
\pi,\frac{2}{3} \pi\right)  $ in the interior of the set $S=\left[
0,2\pi\right]  \times\left[  0,2\pi\right]  $, which are fixed points of the model $F$. There are also four trivial fixed points, $\left(  0,0\right)  $, $\left(  0,2\pi\right)  $, $\left(  2\pi,0\right)  $ and $\left(  2\pi
,2\pi\right)  $ at the corners of $S$, and the four fixed points $\left(
0,\pi\right)  $, $\left(  \pi,2\pi\right)  $, $\left(  2\pi,\pi\right)  $ and $\left(  \pi,0\right)  $ on the edges of $S$.

We now compute the Jacobian matrix $J\left(  x,y\right)  $ to establish the dynamical nature of the fixed points in the usual way.

We have $$J\left(  x,y\right)  =\left[
\begin{array}
[c]{cc}%
1 & 0\\
0 & 1
\end{array}
\right]  + \varepsilon\left[
\begin{array}
[c]{cc}%
2\cos x+\cos\left(  x-y\right)  & -\cos\left(  x-y\right)  +\ \cos y\\
\cos x-\cos\left(  x-y\right)  & \cos\left(  x-y\right)  +2\cos y
\end{array}
\right]  . \label{Jacobmatrix}%$$

We first consider the fixed points of $F$ in the interior of $S$. We start with $\left(  \frac{2}{3} \pi,\frac{4}{3} \pi\right)  $ and $\left(  \frac
{4}{3} \pi,\frac{2}{3} \pi\right)  $. The Jacobian is exactly the same $$\left[
\begin{array}
[c]{cc}%
1-3\frac{\sqrt{3}}{2}\varepsilon & 0\\
0 & 1-3\frac{\sqrt{3}}{2}\varepsilon
\end{array}
\right]  ,$$ meaning that those two points are locally asymptotically stable for $\varepsilon$ sufficiently small.

The Jacobian matrix of $F$ at $\left(  \pi\text{,}\pi\right)  $ is $$\left[
\begin{array}
[c]{cc}%
1-\varepsilon & -2\varepsilon\\
-2\varepsilon & 1-\varepsilon
\end{array}
\right]  ,$$ with eigenvalues $1-3\varepsilon$ and $1+\varepsilon$, which qualifies $\left(  \pi\text{,}\pi\right)  $ as a saddle point. The stable manifold has direction $\left(  1,1\right)  $, and the unstable manifold is tangent at $\left(  \pi\text{,}\pi\right)  $ to the vector $\left(  -1,1\right)  $.

We now consider now the points placed at the vertexes of $S$. The Jacobian matrix of $F$ at $\left(  0\text{,}0\right)  $, $\left(  0,2\pi\right)  $, $\left(  2\pi,0\right)  $ and $\left(  2\pi,2\pi\right)  $ is, for all of them, the following $$\left[
\begin{array}
[c]{cc}%
1+3\varepsilon & 0\\
0 & 1+3\varepsilon
\end{array}
\right]  ,$$ which qualifies all the vertexes of $S$ as a repellers.

On the vertical edges of $S$ we have the fixed points $\left(  0,\pi\right)
$, and $\left(  2\pi,\pi\right)  $, at which the Jacobian matrix of $F$ is $$\left[
\begin{array}
[c]{cc}%
1+\varepsilon & 0\\
2 \varepsilon & 1-3 \varepsilon
\end{array}
\right]  ,$$ which qualifies $\left(  0,\pi\right)  $, and $\left(  2\pi,\pi\right)  $ as saddle points. The stable manifold has the direction of the $y$ axis and the unstable manifold is tangent at $\left(  0\text{,}\pi\right)  $ and $\left(
2\pi,\pi\right)  $ to the vector $\left(  2,1\right)  $.

Finally, at the horizontal edges of $S$ we have the Jacobian matrix of $F$ at $\left(  \pi,0\right)  $, and $\left(  \pi,2\pi\right)  $ $$\left[
\begin{array}
[c]{cc}%
1-3 \varepsilon & 2 \varepsilon\\
0 & 1+ \varepsilon
\end{array}
\right]  ,$$ which qualifies $\left(  \pi,0\right)  $ and $\left(  \pi,2\pi\right)  $ again as saddle points. The stable manifold is the direction of the $x$ axis and the unstable manifold is tangent at $\left(  \pi,0\right)  $ and $\left(  \pi
,2\pi\right)  $ to the vector $\left(  1,2\right)  $.

The local analysis of the fixed points of $F$ reveals a very symmetric picture. When $\varepsilon>0$ is small ($0<\varepsilon<{\varepsilon}_{0}%
=\frac{1}{9}$ is good enough), $F$ is a small perturbation of the identity, $F(\partial{S})=\partial{S}$, the restriction of $F$ to the boundary of $S$, $\partial{S}$, is a bijection (see section 4 for more details), and the Jacobian determinant of $F$ is never null in the interior of $S$. Therefore, $F$ is invertible on S.

Heteroclinic connections and invariant sets
-------------------------------------------

We focus our attention on the existence of invariant subsets of $S$ for the dynamics of $F$. Additionally, we below prove that $S$ is itself an invariant set for the dynamics of $F$.

Recall that an *heteroclinic* (sometimes called a heteroclinic connection, or heteroclinic orbit) is a path in phase space which joins two different equilibrium points. In the sequel, by *sa-heteroclinic*, *rs-heteroclinic*, and *ra-heteroclinic*, we mean an heteroclinic orbit connecting a saddle point to an attractor, an heteroclinic orbit connecting a repeller to a saddle point, and an heteroclinic orbit connecting a repeller to an attractor, respectively.

Let $F$ be our model map in some set $T$ with two fixed points $p$ and $q$. Let $M_{u}\left(  F,p\right)  $ and $M_{s}\left(  F,q\right)  $ be the stable manifold and the unstable manifold ([@AlSaYo]: pages 78, 403) of the fixed points $p$ and $q$, respectively. Then, if by $M$ we denote the heteroclinic connecting $p$ and $q$, we have $$M\subseteq M_{s}\left(  F,p\right)  \cap M_{u}\left(  F,q\right)  .$$ In particular, $M$ is invariant, the $\alpha$-limit and $\omega$-limit sets of the points of $M$ is respectively $p$ and $q$ ([@AlSaYo]: page 331).

The other orbits, i.e., with initial conditions not in $M$, cannot cross the heteroclinic connections when the map $F$ is invertible. In that case, it would be violated the injectivity of the map. In the sequel, we study the heteroclinics that connect saddle points to the attractors. Those heteroclinics determine the nature of all the flow of the dynamical system in the plane, due to the invertible nature of $F$.

### Vertical heteroclinics

Consider the two vertical lateral edges of $S$, $s_{0}$ and $s_{1}$ that are the sets $s_{k}=\left\{  \left(  x,y\right)  \in S:\left(  x=2k\pi\right)
\wedge0\leq y\leq2\pi\right\}  $, $k=0,1$. Consider the image of these segments under $F$. If we write $F=(F_{1},F_{2})$, then $$%
\begin{cases}
F_{1}\left(  2k\pi,y\right)  & =2k\pi+\varepsilon\sin y+\varepsilon\sin\left(
-y\right)  =2k\pi\\
F_{2}\left(  2k\pi,y\right)  & =y+2\varepsilon\sin y+\varepsilon\sin
y=y+3\varepsilon\sin y,
\end{cases}$$ meaning that for $\varepsilon$ small enough the edges $s_{k}$, $k=0,1$, are invariant, as already mentioned in section 2. Because of the initial conditions, on each of the edges $s_{k}$, $k=0,1$, the dynamics is given by$$%
\begin{cases}
x_{n+1} & =2k\pi,\\
y_{n+1} & =y_{n}+3\varepsilon\sin y_{n}\text{,}%
\end{cases}
.$$

For $\varepsilon<\frac{1}{9}$, the map $g:\left[  0,2\pi\right]
\rightarrow\left[  0,2\pi\right]  $ defined by $g(t)=t+3\varepsilon\sin t$ is a homeomorphism from the interval $\left[  0,2\pi\right]  $ into itself, as we can see in figure \[homeo\]. Moreover, since there is an attracting fixed point of this map at $\pi$, the dynamics in the sets $s_{0}$ and $s_{1}$ can be splitted in two subsets where the dynamics is again invariant, which is not very important for our global discussion but establishes that the stable manifolds of the saddle points $\left(  0,\pi\right)  $ and $\left(  2\pi
,\pi\right)  $ are, exactly and respectively, the sets $s_{0}$ and $s_{1}$

\[ptb\]

[Fig3.eps]{}

We have just shown that both $s_{0}$ and $s_{1}$ contain two heteroclinic connections: in $s_{0}$, the line segment $s_{0}^{-}$ from $\left(
0,0\right)  $ to $\left(  0,\pi\right)  $ and $s_{0}^{+}$ from $\left(
0,2\pi\right)  $ to $\left(  0,\pi\right)  $; and in $s_{1}$, the line segment $s_{0}^{-}$ from $\left(  2\pi,0\right)  $ to $\left(  2\pi,\pi\right)  $ and $s_{0}^{+}$ from $\left(  2\pi,2\pi\right)  $ to $\left(  2\pi,\pi\right)  $. The total number of vertical $rs$-heteroclines is $4.$

### Horizontal heteroclinics

Consider the two horizontal top and bottom edges of $S$, $r_{0}$ and $r_{1}$, that are the sets $r_{k}=\left\{  \left(  x,y\right)  \in S:0\leq x\leq
2\pi\wedge\left(  y=2k\pi\right)  \right\}  $, $k=0,1$. Consider the image of these segments under $F$. As before, if we write $F=(F_{1},F_{2})$, then $$%
\begin{cases}
F_{1}\left(  x,2k\pi\right)  & =x+3\varepsilon\sin x\\
F_{2}\left(  x,2k\pi\right)  & =2k\pi,
\end{cases}$$ meaning that, for $\varepsilon$ small enough, the edges $r_{k}$, $k=0,1$, are invariant. Because of the initial conditions, on each of the edges $s_{k}$, $k=0,1$, the dynamics is given by$$%
\begin{cases}
x_{n+1} & =x_{n}+3\varepsilon\sin x_{n}\\
y_{n+1} & =2k\pi\text{.}%
\end{cases}$$ For $\varepsilon<\frac{1}{9}$, the map $g:\left[  0,2\pi\right]
\rightarrow\left[  0,2\pi\right]  $, defined by $g(t)=t+3\varepsilon\sin t$, is the same occurred before, now involved in the dynamics in the invariant edges $r_{0}$ and $r_{1}$. The stable manifolds of $\left(  \pi,0\right)  $ and $\left(  \pi,2\pi\right)  $ are again, respectively, the edges $r_{0}$ and $r_{1}$.

Arguing as for $s_{0}$ and $s_{1}$, we have that both, $r_{0}$ and $r_{1}$, contain two analogous heteroclinic connections.

We have just proved, in detail, that the boundary of $S$ is an invariant set. More is true: each edge of $S$ is an invariant set.

Since the map $F$ is invertible, the initial conditions in the interior of $S
$, $S^{0}$, cannot cross the invariant boundary $\partial S=s_{0}\cup
s_{1}\cup r_{0}\cup r_{1}$, meaning that ${S}^{0}$ is an invariant set. This means, in particular, that for equal clocks there will be no secular drift of phase differences of the three clocks, the delays and advances are contained in the set $S=\left[  0,2\pi\right]  \times\left[  0,2\pi\right]  $.

The total number of horizontal $rs$-heteroclines is $4$. The total number of $rs$-heteroclinics in the boundary of $S$ is $8$.

### Diagonal heteroclinics

Finally, we now show that, $S^{o}$, the interior set of $S$, can be splitted in two subsets, $S_{U}$ and $S_{D}$, $U$ for up and $D$ for down, where the dynamics is again invariant. Consider now the set $$\Delta=\left\{  \left(  x,y\right)  \in S:y=x\text{, }x\in\left[
0,2\pi\right]  \right\}  ,$$ the diagonal of $S$ connecting $\left(  0,0\right)  $ to $\left(  2\pi
,2\pi\right)  $. The image of a point of $\Delta$ by $F$ is now $$%
\begin{cases}
F_{1}\left(  x,x\right)  & =x+3\varepsilon\sin x,\\
F_{2}\left(  x,x\right)  & =x+3\varepsilon\sin x.
\end{cases}$$ Hence, the same homeomorphism $g$ as before appears again. We repeat the same reasonings as before and deduce that $\Delta$ is invariant under $F$, and it splits $S^{o}$ in two open sets: the triangle above it and the triangle below it. Moreover, the stable manifold of the saddle point $\left(  \pi,\pi\right)
$ is the set $\Delta$.

This also proves the existence of two heteroclinics in $\Delta$, connecting $\left(  0,0\right)  $ to $\left(  \pi,\pi\right)  $ and $\left(  2\pi
,2\pi\right)  $ to $\left(  \pi,\pi\right)  $, respectively. The total number of $rs$-heteroclines is now $10$, respectively $8$ on the edges and $2$ on the main diagonal $\Delta$, all of them connecting repellers to saddles.

Consider now the other diagonal of $S$, i.e., the set $$\tilde{\Delta}=\left\{  \left(  x,y\right)  \in S:y=2\pi-x\text{, }x\in\left[
0,2\pi\right]  \right\}  .$$ The image of a point of $\tilde{\Delta}$ under $F$ now is $$%
\begin{cases}
F_{1}\left(  x,y\left(  x\right)  \right)  & =x+\varepsilon\sin x+\varepsilon
\sin2x,\\
F_{2}\left(  x,y\left(  x\right)  \right)  & =2\pi-\left(  x+\varepsilon\sin
x+\varepsilon\sin2x\right)  .
\end{cases}$$ Hence, $\tilde{\Delta}$ is invariant.

The map $h_{1}:\left[  0,2\pi\right]  \rightarrow\left[  0,2\pi\right]  $, defined as $h_{1}(t)=t+\varepsilon\sin t+\varepsilon\sin2t$, is a homeomorphism with $5$ fixed points from $\left[  0,2\pi\right]  $ to itself (see figure \[homeo2\]).

We repeat the same reasonings as before and deduce that the set $\tilde
{\Delta}$ splits the interior set $S^{o}$ again in two open sets: the triangle above and the triangle below. So, now we have splitted $S^{0}$ in four small triangles.

There are four heteroclinic connections in $\tilde{\Delta}$, one connecting the repeller $\left(  0,2\pi\right)  $ to the attractor $\left(  \frac{2\pi
}{3},\frac{4\pi}{3}\right)  $ (ra-heteroclinic), two $sa$-heteroclinics connecting the saddle point $\left(  \pi,\pi\right)  $ to the attractors $\left(  \frac{2\pi}{3},\frac{4\pi}{3}\right)  $ and $\left(  \frac{4\pi}%
{3},\frac{2\pi}{3}\right)  $, and, finally, the last heteroclinic on this diagonal set is the one that connects the repeller $\left(  2\pi,0\right)  $ to the attractor $\left(  \frac{4\pi}{3},\frac{2\pi}{3}\right)  $ (ra-heteroclinic). The total number of sa-heteroclinics is now $2$.\[ptb\]

[Fig4.eps]{}

We proceed with the same line of reasoning for the other $sa$-heteroclinics.

Consider now the set $$d_{1}=\left\{  \left(  x,y\right)  \in S:y=\pi+\frac{x}{2}\text{, }x\in\left[
0,\frac{2\pi}{3}\right]  \right\}$$ and the map $F$ applied to the points of $d_{1}$:$$%
\begin{cases}
F_{1}\left(  x,y\left(  x\right)  \right)  & =x+2\varepsilon\sin
x-2\varepsilon\sin\left(  \frac{x}{2}\right)  ,\\
F_{2}\left(  x,y\left(  x\right)  \right)  & =\pi+\frac{1}{2}\left(
x+2\varepsilon\sin x-2\varepsilon\sin\left(  \frac{x}{2}\right)  \right)  .
\end{cases}$$ The points of $d_{1}$ stay in $d_{1}$ under the action of $F$, proving that this set also is invariant. The function $h_{2}:[0,2\frac{\pi}{3}%
]\rightarrow\lbrack0,2\frac{\pi}{3}]$, defined as $h_{2}(t)=t+2\varepsilon\sin
t-2\varepsilon\sin\left(  \frac{t}{2}\right)  $, is a homeomorphism, from which we can readily see that the dynamics in $d_{1}$ is quite simple. The graph of this homeomorphism can be seen in figure \[homeo3\]. There is one $sa$-heteroclinic from the saddle at $\left(  0,\pi\right)  $ to the attractor $\left(  \frac{2\pi}{3},\frac{4\pi}{3}\right)  $. Actually, there is another heteroclinic in the segment connecting the repeller $\left(  2\pi,2\pi\right)
$ to the attractor $\left(  \frac{2\pi}{3},\frac{4\pi}{3}\right)  $, but this is not an sa-heteroclinic. Up to now, we have $3$ $sa$-heteroclinic connections.

\[ptb\]

[Fig5.eps]{}

Consider now the set $c_{1}=\left\{  \left(  x,y\right)  \in S:y=2x\text{,
}x\in\left[  \frac{2\pi}{3},\pi\right]  \right\}  $ and the map $F$ applied to the points of $c_{1}$$$\begin{aligned}
F_{1}\left(  x,y\right)   &  =x+\varepsilon\sin x+\varepsilon\sin2x,\\
F_{2}\left(  x,y\right)   &  =y+2\varepsilon\sin x+2\varepsilon\sin2x,\end{aligned}$$ the points of $c_{1}$ stay in $c_{1}$ under $F$, proving that this set is invariant. Actually, the segment would be invariant if we extended $x$ to the interval $\left[  0,\pi\right]  $, but we are not interested in heteroclinics from repellers to attractors. Moreover, the dynamics is given by a restriction of $h_{1}$ to the interval $\left[  \frac{2\pi}{3},\pi\right]  $. In this interval there are only two fixed points, the attractor $\frac{2\pi}{3}$ and the repeller $\pi$. This procedure adds one more $sa$-heteroclinic to the global picture. So, we have found, up to now, $4$ $sa$-heteroclinics.

In $S_{D}$, we consider $$c_{2}=\left\{  \left(  x,y\right)  \in S:y=2\left(  x-\pi\right)  \text{,
}x\in\left[  \pi,\frac{4\pi}{3}\right]  \right\}$$ and $$d_{2}=\left\{  \left(  x,y\right)  \in S:y=\frac{x}{2}\text{, }x\in\left[
\frac{4\pi}{3},2\pi\right]  \right\}  .$$ Following exactly the same reasonings as before, we obtain two more $sa$-heteroclinics, one connecting $\left(  \pi,0\right)  $ to the attractor $\left(  \frac{4\pi}{3},\frac{2\pi}{3}\right)  $ and the other connecting $\left(  2\pi,\pi\right)  $ to the same attractor.

Phase portrait
--------------

The total number of $sa$-heteroclinics is $6$. All of them are straight segments. The other $8$ $sa$-heteroclinics split the set $S$ in six invariant sets as can be seen in Figure \[Portrait\], where the red curves represent saddle-node heteroclines. The flow curves represented in the phase portrait. Since the map $F$ is invertible, no orbit can cross either the red curves, blue curves or black flow curves. There are only two attractors and the dynamics, due to the invertible nature of the map $F$ and its large symmetry, is relatively simple: in every invariant set in the plane, the restriction maps are again homeomorphisms and the flow curves must follow, by continuity, the heteroclinic connections on the outer boundaries of each invariant set.

Consequently, only the orbits on the outer edges and main diagonal, i.e., in the set $s_{0}\cup s_{1}\cup r_{0}\cup r_{1}\cup d$ are not attracted to the two attractors $\left(  \frac{2\pi}{3},\frac{4\pi}{3}\right)  $ and $\left(
\frac{4\pi}{3},\frac{2\pi}{3}\right)  $. The upper attractor $\left(
\frac{2\pi}{3},\frac{4\pi}{3}\right)  $ attracts the points in the open upper triangle $S_{U}$ with converse results for the lower attractor $\left(
\frac{4\pi}{3},\frac{2\pi}{3}\right)  $ in $S_{D}$. The full picture can be seen in figure \[Portrait\].

Conclusions and future work
===========================

In this paper we have proved that three oscillators, mutually interacting with symmetric coupling, converge to a final symmetric locked state with mutual phase differences of $\frac{2\pi}{3}$, this can happen in two different settings, clockwise or counterclockwise, depending on the initial conditions.

This very symmetrical final locked state induces us to consider the conjecture that $n$ oscillators weakly interacting with all the others $n-1$ oscillators will reach a final state with mutual phase differences of $\frac{2\pi}{n}$ clockwise or counterclockwise distributed.

In future work, already in preparation, we shall discuss the same phenomenon with slightly different natural angular frequencies $\omega_{1}$, $\ \omega_{2}$ and $\omega_{3}$ and, in particular, the existence and form of Arnold Tongues [@boyland1986; @gilmore2011].

As done for [@OlMe], it would be interesting to check experimentally our model, to see if the real world matches the theoretical predictions.

#### Aknowledgements {#aknowledgements .unnumbered}

The author ED was partially supported by the program Erasmus+. The author HMO was partially supported by FCT/Portugal through the project UID/MAT/04459/2013.

Appendix {#appendix .unnumbered}
========

In this Appendix, as another way to determine the nature of the two attractors, that is their global asymptotical stability, we point out the existence, in the sets $S_{U}$ and $S_{D}$, of two Liapounov functions, $V_{U}$ and $V_{D}$, respectively [@lasalle1976stability].

Consider, first, the invariant set $S_{U}$. Define $V_{U}: S_{U}
\rightarrow\mathbb{R}$ as follows: $$V_{U}\left(  x,y\right)  =\left(  x-\frac{2\pi}{3}\right)  ^{2}+\left(
y-\frac{4\pi}{3}\right)  ^{2}-\left(  x-\frac{2\pi}{3}\right)  \left(
y-\frac{4\pi}{3}\right)  \text{.}$$ The *discrete orbital derivative* inside the invariant set $S_{U}$ is $$DF\left(  x,y\right)  =V_{U}\left(  F\left(  x,y\right)  \right)  -V_{U}
\left(  x,y\right)  .$$

We have $$\begin{aligned}
{DF\left(  x,y\right)  }  &  =V_{U}\left(  F\left(  x,y\right)  \right)
-V_{U} \left(  x,y\right) \\
&  =V_{U}(x+\varepsilon\phi(x,y),y+\varepsilon\gamma(x,y))-V_{U}(x,y)\\
&  ={\varepsilon}^{2}[{\phi}^{2}(x,y)+{\gamma}^{2}(x,y)-{\phi}(x,y)\cdot
{\gamma}(x,y)]\\
&  +\epsilon\lbrack(x-\frac{2}{3}\pi)(2\phi(x,y)-\gamma(x,y))+(y-\frac{4}%
{3}\pi)(2\gamma(x,y)-\phi(x,y)],\\
&\end{aligned}$$ where, recall that $$%
\begin{cases}
\varphi\left(  x,y\right)  & =2\sin x+\sin y+\sin\left(  x-y\right) \\
\gamma\left(  x,y\right)  & =\sin x+2\sin y+\sin\left(  y-x\right)
=\phi(y,x).
\end{cases}$$ An easy computation shows that $$%
\begin{cases}
2\varphi\left(  x,y\right)  -\gamma(x,y) & =3(\sin x+\sin\left(  x-y\right)
)\\
2\gamma\left(  x,y\right)  -\phi(x,y) & =3(\sin y-\sin\left(  x-y\right) \\
\phi(x,y)+\gamma(x,y) & =3(\sin x+\sin y).
\end{cases}$$ Hence, $$\begin{aligned}
\frac{DF(x,y)}{\varepsilon}  &  =\varepsilon\lbrack{(\phi(x,y)+\gamma
(x,y))}^{2}-3\phi(x,y)\cdot\gamma(x,y)]\\
&  =3\varepsilon\lbrack3{(\sin x+\sin y)}^{2}-\phi(x,y)\cdot\gamma(x,y)]\\
&  +(x-\frac{2}{3}\pi)\cdot(\sin x+\sin(x-y))+(y-\frac{4}{3}\pi)\cdot(\sin
y-\sin(x-y))\\
&\end{aligned}$$

By using numerical analysis, we conclude that this discrete orbital derivative, $DF$, is non-positive in $S_{U}$ for small $\varepsilon$, more precisely, zero for the fixed points and negative elsewhere.

Analogously, by the same method, the function $V_{D}:S_{D}\rightarrow R$ defined as

**$$V_{D}\left(  x,y\right)  =\left(  x-\frac{4\pi}{3}\right)  ^{2}+\left(
y-\frac{2\pi}{3}\right)  ^{2}-\left(  x-\frac{4\pi}{3}\right)  \left(
y-\frac{2\pi}{3}\right)  \text{,}%$$**

is a Liapounov function on $S_{D}$.

[99]{}

R. Abraham. , pages 49–78. World Scientific, Singapore, 1991.

R. Abraham and A. Garfinkel. The dynamics of synchronization and phase regulation. 2003. http://www.ralph-abraham.org/articles/Blurbs/blurb111.shtml.

R. Adler. A study of locking phenomena in oscillators. , 34(6):351–357, 1946.

K. T. Alligood, T. D. Sauer, and J. A. Yorke. . Springer, 1997.

A. A. Andronov, A. A. Vitt, and S. E. Khaikin. . Pergammon Press, Oxford, New York, 1959/1963/1966.

D. K. Arrowsmith, C. M. Place, C. Place, et al. . Cambridge university press, 1990.

M. Bennett, M. Schatz, H. Rockwood, and K. Wiesenfeld. Huygen’s clocks. , 458(2019):563–579, 2002.

G. D. Birkhoff. . American Mathematical Society, Providence, Rhode Island, 1950.

P. L. Boyland. Bifurcations of circle maps: Arnol’d tongues, bistability and rotation intervals. , 106(3):353–381, 1986.

S. R. Campbell. . PhD thesis, The Ohio State University, 1997.

S. R. Campbell, D. L. Wang, and C. Jayaprakash. Synchrony and desynchrony in integrate-and-fire oscillators. , 11(7):1595–1619, 1999.

K. Czolczynski, P. Perlikowski, A. Stefanski, and T. Kapitaniak. Huygen’s odd sympathy experiment revisited. , 07(21):2047–2056, 2011.

A. Fradkov and B. Andrievsky. Synchronization and phase relations in the motion of two-pendulum system. , 6(42):895–901, 2007.

R. Gilmore and M. Lefranc. . WILEY-VCH Verlag GmbH [&]{} Co. KGaA, Weinheim, Germany, 2 edition, 2011.

J. Guckenheimer. Isochrons and phaseless sets. , 1(3):259–273, 1975.

C. Huygens. . Societe Hollandaise Des Sciences, Martinus Nijho, La Haye, 1895.

V. Jovanovic and S. Koshkin. Synchronization of huygens’ clocks and the poincare method. , 12(331):2887–2900, 2012.

M. Kapitaniak, K. Czolczynski, P. Perlikowski, A. Stefanski, and T. Kapitaniak. Synchronization of clocks. , 1(517):1–69, 2012.

Y. Kuramoto. Self-entrainment of a population of coupled non-linear oscillators. In *International symposium on mathematical problems in theoretical physics*, pages 420–422. Springer, 1975.

J. P. LaSalle. , volume 25. Siam, 1976.

E. A. Martens, S. Thutupalli, A. Fourrièrec, and O. Hallatschek. Chimera states in mechanical oscillator networks. , 26(110):10563–10567, 2013.

R. E. Mirollo and S. H. Strogatz. Synchronization of pulse-coupled biological oscillators. , 50(6):1645–1662, 1990.

H. Nakao. Phase reduction approach to synchronisation of nonlinear oscillators. , 57(2):188–214, 2016.

H. M. Oliveira and L. V. Melo. Huygens synchronization of two clocks. , 5(11548):1–12, 2015. doi: 10.1038/srep11548.

W. T. Oud, H. Nijmeijer, and A. Y. Pogromsky. , volume 336 of *Lecture Notes in Control and Information Science*, pages 191–203. Springer, Berlin Heidelberg, 2006.

A. Pikovsky, M. Rosenblum, and J. Kurths. , volume 12 of *Cambridge Nonlinear Science Series*. Cambridge University Press, Cambridge, 1 edition, 5 2003.

M. Senator. Synchronization of two coupled escapement-driven pendulum clocks. , 3–5(291):566–603, 2006.

S. Strogatz. . Penguin UK, 2004.

S. H. Strogatz. From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators. , 143(1-4):1–20, 2000.

J. Vassalo-Pereira. A theorem on phase-locking in two interacting clocks (the [H]{}uygens effect). In A. Avez, A. Blaquiere, and A. Marzollo, editors, *Dynamical Systems and Microphysics: Geometry and Mechanics*, pages 343–352, New York, London, 1982. Academic Press.
