# Bifurcations of heteroclinic contours in two-parameter planar systems:
Overview and explicit examples

Yuri A. Kuznetsov and Joost Hooyman  Mathematical Institute, Utrecht
University,
P.O. Box 80010, 3508 TA Utrecht,
The Netherlands

((to be inserted by publisher))

###### Abstract

Smooth planar vector fields containing two hyperbolic saddles may possess
contours formed by heteroclinic connections between these saddles. We present
an overview of the bifurcations of these contours based on papers by J.W. Reyn
and A. Dukov. Additionally, two new explicit polynomial systems containing
such contours are derived, which are studied using the bifurcation software
matcont and are shown to exhibit the theoretically predicted phenomena,
including series of heteroclinic connections.

###### keywords:

Heteroclinic contours, global bifurcations, planar systems

##  1 Introduction

Consider a vector field in the plane containing two hyperbolic saddles and
suppose an orbit leaving one saddle connects to the other saddle and vice
versa. This configuration of saddles and orbits is called a heteroclinic
contour (other names used: polygone, polycycle). This paper deals with
bifurcations occurring in generic two-parameter families of vector fields
containing such contours at some critical parameter values.

Heteroclinic contours can be subdivided into two categories, namely monodromic
and non-monodromic. For both cases there already exist separate theoretical
studies. Results on the monodromic case have been published by Reyn [1980] and
are quite well known. For example, they are included in the exposition Bazykin
_et al._ [1989] and reproduced in the encyclopedic work Arnold _et al._
[1994]. A much more recent publication by Dukov [2018] is devoted to the non-
monodromic case.

In both cases, depending on the eigenvalues of the saddles, either one limit
cycle or two limit cycles of opposite stability are generated as the result of
homoclinic bifurcations. When two limit cycles exist these may collide to form
a semi-stable limit cycle at a fold (or saddle-node) of cycles bifurcation.

The non-monodromic case exhibits an additional phenomenon of flashing (or
sparkling) heteroclinic connections, which may be described as an infinite
series of bifurcation curves corresponding to heteroclinic connections from
one saddle to the other making arbitrarily many windings around one of the
saddles. This phenomenon was already discovered in 1981 by Malta & Palis
[1981] but in a less general (one-parameter) context. Two-parameter
perturbation of the non-monodromic case provides an example of a situation
which is only two-dimensional but nevertheless possesses a bifurcation diagram
of considerable complexity, where an infinite number of bifurcation curves
emanate from a codim 2 point. Notice, however, that a possibility of such
parametric complexity near a non-monodromic contour was already briefly
mentioned in Bazykin _et al._ [1989] (and probably known to experts earlier).

The first aim of this paper is to collect the results of both papers and
present them in a unified manner by describing generic bifurcation diagrams in
both the monodromic and the non-monodromic cases.

The second aim of this paper is to derive two polynomial systems which contain
heteroclinic connections of the monodromic and non-monodromic kind
respectively. These systems may be perturbed to bifurcate as predicted in the
theoretical part. The method to obtain the monodromic example is adapted from
a method to derive a polynomial system containing a homoclinic connection by
Sandstede [1997]. The non-monodromic example is derived by modifying a
reversible system containing a homoclinic connection. After these explicit
polynomial systems are obtained we use the standard numerical software package
matcont Dhooge _et al._ [2003, 2008]; De Witte _et al._ [2012] to continue
relevant phase objects in two parameters, thus numerically producing
bifurcation diagrams of these systems.

In Appendix A, we summarize for reader’s convenience the major facts on
flashing heteroclinic bifurcations, assuming that basic facts on codim 1
bifurcations of cycles and connecting orbits of planar vector fields are known
(see, e.g. Kuznetsov [2004]). In Appendix 25 we show that the splitting of the
heteroclinic connections in the monodromic example is regular.

##  2 Heteroclinic contours of codim 2

Consider a planar vector field

$\left\\{\begin{array}[]{rcl}\dot{x}&=&f(x,y),\\\
\dot{y}&=&g(x,y),\end{array}\right.$ (1)

where $f$ and $g$ are smooth functions from $\mathbb{R}\times\mathbb{R}$ to
$\mathbb{R}$. Assume the system (1) has two hyperbolic equilibria of saddle
type which we label $L$ and $M$. Let $\lambda_{s}<0<\lambda_{u}$ and
$\mu_{s}<0<\mu_{u}$ be the eigenvalues of $L$ and $M$, respectively. Introduce
the corresponding _saddle indices_ by

$\lambda:=-\frac{\lambda_{s}}{\lambda_{u}},~{}~{}\mu:=-\frac{\mu_{s}}{\mu_{u}}.$

Also suppose that this system has solutions $\gamma_{1}(t)$ satisfying

$\lim_{t\to-\infty}\gamma_{1}(t)=L\quad\text{and}\quad\lim_{t\to\infty}\gamma_{1}(t)=M,$

and $\gamma_{2}(t)$ satisfying

$\lim_{t\to-\infty}\gamma_{1}(t)=M\quad\text{and}\quad\lim_{t\to\infty}\gamma_{1}(t)=L.$

We see that the orbits of $\gamma_{1}(t)$ and $\gamma_{2}(t)$ are heteroclinic
connections between $L$ and $M$ in opposite directions. Together with the
saddles they form a heteroclinic contour.

(a) Monodromic (b) Non-monodromic

Figure 1: The two types of heteroclinic contours

Depending on the position of the stable and unstable manifolds of $L$ and $M$
there are two distinct types of heteroclinic contours called _monodromic_ and
_non-monodromic_ , see Fig. 1.

To understand this choice of terminology, note that the monodromic contour may
be approached as a limit set by solutions inside the contour. It is also
possible to contain all stable and unstable manifolds inside the contour, in
which case the contour may be approached as a limit set for solutions outside
the contour. This case is also monodromic but is topologically equivalent to
the depicted one and, thus, will not be considered. In the non-monodromic
case, the contour can not be approached as a limit set, neither from the
inside nor the outside.

We shall see that the bifurcation diagrams in both cases depend strongly on
the values of $\lambda$, $\mu$ and also $\lambda\mu$. Especially, it is
important whether these quantities are smaller or larger than 1. Recall that
saddles with index greater than one are called _dissipative_. In principle,
the following generic subcases are possible:

1 | $\lambda<1$ | $\mu>1$ | $\lambda\mu<1$
---|---|---|---
2 | $\lambda<1$ | $\mu>1$ | $\lambda\mu>1$
3 | $\lambda>1$ | $\mu<1$ | $\lambda\mu>1$
4 | $\lambda>1$ | $\mu<1$ | $\lambda\mu<1$
5 | $\lambda>1$ | $\mu>1$ | $\lambda\mu>1$
6 | $\lambda<1$ | $\mu<1$ | $\lambda\mu<1$

However, reversing time and interchanging the roles of $L$ and $M$, we see
that there are only two essentially different subcases to consider. More
specifically, reversing time, we see that 1 goes to 3 and 5 to 6. When
changing the roles of $L$ and $M$, 1 goes to 4 and 2 to 3. Thus if we study
only the subcases 1 and 6, the results for all others follow immediately. To
conclude, we consider only subcases $\lambda<1$, $\mu<1$ and $\lambda<1$,
$\mu>1$ with $\lambda\mu<1$ in both of them.

##  3 Bifurcations of heteroclinic contours

We are now ready to describe bifurcations of planar monodromic and non-
monodromic heteroclinic contours in generic two-parameter systems. Consider a
smooth vector field

$\left\\{\begin{array}[]{rcl}\dot{x}&=&f(x,y,\alpha),\\\
\dot{y}&=&g(x,y,\alpha),\end{array}\right.$ (2)

where $f$ and $g$ are smooth functions from
$\mathbb{R}\times\mathbb{R}\times\mathbb{R}^{2}$ to $\mathbb{R}$. We assume
that at $\alpha=(0,0)$ the system (2) has a monodromic or non-monodromic
heteroclinic contour described in Sec. 2. Each case is treated separately.

###  3.1 Bifurcation diagrams in the monodromic cases

In a neighborhood of $L$, we introduce two cross-sections, $\Sigma_{L}$ and
$\Pi_{L}$, on the incoming and outgoing heteroclinic connection respectively,
see Fig. 2. On these cross sections, we introduce coordinates $\xi_{L}$ and
$\eta_{L}$, again respectively. We define these coordinates in such a way that
the points $\xi_{L}=0$ and $\eta_{L}=0$ coincide with the points on the
intersection of the stable and unstable manifolds and the corresponding cross-
section, and such that positive values of the coordinates correspond with
points inside the heteroclinic contour.

Figure 2: The monodromic contour with all cross-sections and coordinates

Analogously, in a neighborhood $U_{M}$ of $M$, we introduce cross-sections
$\Sigma_{M}$ and $\Pi_{M}$ furnished with coordinates $\xi_{M}$ and $\eta_{M}$
with analogous properties.

Using the cross-sections, we define the splitting parameters
$\beta=(\beta_{1},\beta_{2})$ which measure the splitting of the heteroclinic
connection along $\Sigma_{L}$ and $\Sigma_{M}$ respectively.

Summarising the results for both monodromic subcases, we may now sketch their
bifurcation diagrams in Figs. 3 and 4. We denote the curves corresponding to
homoclinic connections $P_{L}$ and $P_{M}$. In the second subcase, the curve
in the parameter plane corresponding to the semi-stable (double) limit cycle
will be denoted by $F$. Curves corresponding to heteroclinic connections
between the saddles are denoted with the letter $H$ and are labeled by the
saddle from which they leave. Thus the curve corresponding to a heteroclinic
connection from $L$ to $M$ is labeled with $H_{L}$.

Choosing initial parameter values on the curve $P_{L}$ and varying the
parameters through the bifurcation diagram in Fig. 3 counterclockwise, we
observe the following changes of the phase portrait in the first subcase. We
begin with a small homoclinic connection of the saddle $L$. This homoclinic
connection is broken outwards and an unstable limit cycle is generated in
accordance with the fact that $\lambda<1$. Then, a heteroclinic connection
from $L$ to $M$ is formed and destroyed, and thereafter the same happens to a
heteroclinic connection from $M$ to $L$. Eventually, the unstable limit cycle
disappears via a small homoclinic connection of $M$, which is then broken
inwards. This happens in accordance with $\lambda<1$. Under further variation
of the parameters, heteroclinic connections from $L$ to $M$ and back are
formed and destroyed before the original homoclinic connection of $L$ is
regained. In short, one might describe these changes as an interplay between
two homoclinic bifurcations.

In the second subcase, we perform the same analysis. Starting with parameter
values on $P_{L}$ and varying them counterclockwise through the bifurcation
diagram in Fig. 4, we begin with a small homoclinic connection of $L$ which is
destroyed forming an unstable limit cycle. This behaviour is not different
from the first subcase, which was to be expected as it still holds that
$\lambda<1$. The evolution does not differ from the first subcase until a
small homoclinic connection of $M$ is formed. In this process, the unstable
limit cycle does not disappear but is contained within the homoclinic
connection. This happens in accordance with the fact that $\mu>1$.
Furthermore, when this connection is broken inwards another limit cycle, now
stable, is generated enclosing the original unstable one.

Figure 3: $\lambda<1$ and $\mu<1$ Figure 4: $\lambda<1$ and $\mu>1$

Eventually, these limit cycles collide, forming a single double limit cycle
which is stable from the outside but unstable from the inside. This semi-
stable and nonhyperbolic limit cycle is immediately destroyed and further
variation of the parameters only creates heteroclinic connections from $L$ to
$M$ and back, before the original homoclinic connection of $L$ is again
formed. In short, one might describe these changes as the interplay between
two homoclinic bifurcations resulting in a fold of cycles bifurcation.

The above results can be obtained by considering a model Poincaré map that is
constructed as follows. We may define approximate singular maps
$\Delta_{L}:\Sigma_{L}\to\Pi_{L}$ and $\Delta_{M}:\Sigma_{M}\to\Pi_{M}$ by

$\Delta_{L}(\xi_{L})=\xi_{L}^{\lambda}\quad\text{and}\quad\Delta_{M}(\xi_{M})=\xi_{M}^{\mu}$

near the saddles and truncated regular maps $Q:\Pi_{L}\to\Sigma_{M}$ and
$R:\Pi_{M}\to\Sigma_{L}$ by

$Q(\eta_{L})=\beta_{2}+\theta_{1}\eta_{L}\quad\text{and}\quad
R(\eta_{M})=\beta_{1}+\theta_{2}\eta_{M},$

where $\theta_{1,2}$ are positive constants, near the heteroclinic
connections. By composing these maps we arrive at a model Poincaré map
$P:\Sigma_{L}\to\Sigma_{L}$, which becomes

$\begin{split}P(\xi_{L})&=(R\circ\Delta_{M}\circ Q\circ\Delta_{L})(\xi_{L})\\\
&=\beta_{1}+\theta_{2}(\beta_{2}+\theta_{1}\xi_{L}^{\lambda})^{\mu}.\end{split}$
(3)

It has been shown by Reyn [1980] that taking into account higher-order terms
in the singular and regular maps in (3) does not change topology of the
bifurcation diagrams in both subcases.

###  3.2 Bifurcation diagrams in the non-monodromic cases

We now proceed to the non-monodromic case. The treatment will be analogous to
that of the monodromic case. Thus, we begin by describing homoclinic
connections, then move to limit cycles and their bifurcations, and finally to
heteroclinic connections.

The construction of the model maps in the non-monodromic case is virtually the
same as in the monodromic case. However, to obtain a singular map describing
the behaviour of a solution near the saddle, the cross-sections at $L$ must be
oriented in the opposite direction to allow solutions to pass near the saddles
on the positive side of the cross-sections.

The new configuration of the cross sections is pictured in Fig. 5.
Consequently, the splitting parameters at $L$ are now also measured
differently.

Figure 5: The non-monodromic contour with all cross-sections and coordinates

The bifurcation diagrams in the two non-monodromic cases are sketched in Figs.
6 and 7. Once more, we will be able to find two curves, $P_{L}$ and $P_{M}$,
corresponding to homoclinic connections at $L$ and $M$ which border a wedge in
the parameter-plane.

In the first subcase ($\lambda<1$ and $\mu<1$) we find that outside the wedge
no limit cycles are present. If we continuously vary our parameters in this
region, we can enter the wedge by crossing its border either through $P_{L}$
or $P_{M}$. Breaking either one of the associated homoclinic connections
generates an unstable limit cycle via a homoclinic bifurcation. This limit
cycle exists in the wedge between $P_{L}$ and $P_{M}$. If we vary the
parameters further, a homoclinic connection of the opposing saddle is created
and then destroyed, also destroying the unstable limit cycle.

In the second subcase ($\lambda<1$ and $\mu>1$) the wedge between $P_{L}$ and
$P_{M}$ still contains an unstable limit cycle that is generated by breaking
the homoclinic connection of $L$. Varying the parameters further, a homoclinic
connection appears for $M$, now with $\mu>1$. As a result of this difference,
breaking of the this connection does not destroy the unstable limit cycle, but
generates an additional stable limit cycle enclosing it. Eventually, these
limit cycles collide, momentarily forming a nonhyperbolic semi-stable double
limit cycle at a cyclic fold bifurcation at the curve $F$.

We recognize that all these results are similar to the monodromic case.
However, the non-monodromic case also presents a new phenomenon: non-trivial
heteroclinic connections, i.e. heteroclinic connections which appear when both
splitting parameters are non-zero. These connections must wind around $L$ at
least once before reaching $L$ or $M$. Actually, in both subcases, there exist
two infinite series of curves corresponding to heteroclinic connections
between the saddles, with each connection making an increasing number of
turns. Thus, we encounter the phenomenon of _flashing heteroclinic
connections_ , both at the bifurcation of a semi-stable limit cycle as well as
at a homoclinic connection (see Appendix A). In the parameter plane, the
infinite series of curves corresponding to these non-trivial heteroclinic
connections accumulate, in the first subcase at $P_{L}$ and $P_{M}$,
corresponding to the homoclinic connections, and in the second subcase at
$P_{L}$ and $F$.

Summarising the results by Dukov [2018] on both non-monodromic subcases, we
present the bifurcation diagrams in Figs. 6 and 7. The curves corresponding to
heteroclinic connections are indicated along with the number of turns around
the saddle $L$, shown as a bracketed superscript. For practical reasons, we
only supply phase portraits for heteroclinic connections making at most one
turn. The accumulating curves corresponding to the connections making more
turns are drawn dashed.

Choosing parameter values on the curve $H_{L}^{(0)}$ and varying the
parameters through the bifurcation diagram in Fig. 6 counterclockwise, we
observe the following changes of the phase portrait. Initially, there is a
heteroclinic connection from $L$ to $M$ which is destroyed. After variation of
the parameters, a second heteroclinic connection is formed which winds around
the saddle $L$ once before reaching $M$. Varying the parameters further
results in an infinite series of flashing heteroclinic connections from $L$ to
$M$ making an increasing number of turns around $L$ before eventually a big
homoclinic connection of $L$ is formed. This big homoclinic connection is
broken inwards, generating an unstable limit cycle in accordance with the fact
that $\lambda<1$. This limit cycle disappears when a small homoclinic
connection of $M$ is formed, in accordance with $\mu<1$. Breaking this
homoclinic connection results again in an infinite series of flashing
heteroclinic connections, now from $M$ to $L$ and making a decreasing number
of turns. After further variation of the parameters the original heteroclinic
connection from $L$ to $M$ is regained.

Figure 6: $\lambda<1$ and $\mu<1$ Figure 7: $\lambda<1$ and $\mu>1$

In short, the changes of this situation may be described as an interplay
between two homoclinic bifurcations in the presence of two series of flashing
heteroclinic connections.

The analysis of the second subcase (see Fig. 7) yields the same result up to
the appearance of a small homoclinic connection of $M$. In accordance with the
fact that $\mu>1$, the existing unstable limit cycle is not destroyed in this
process. As in the monodromic case, the destruction of this small homoclinic
orbit generates a stable limit cycle enclosing the stable one. Eventually
these two cycles collide, momentarily forming a double limit cycle which is
stable from the outside but unstable from the inside. In this case, an
infinite series of heteroclinic connections from $M$ to $L$ accumulates at
this semi-stable nonhyperbolic limit cycle. Thus, when this limit cycle is
destroyed we again encounter an infinite series of flashing heteroclinic
connections from $M$ to $L$ making a decreasing number of turns around $L$.
The further behaviour is the same as in the first subcase. These changes may
be described as the interplay between two homoclinic bifurcations yielding a
fold of cycles bifurcation in the presence of two series of flashing
heteroclinic connections.

All above results can be obtained by considering a model Poincaré map that is
constructed similar to the monodromic case. The approximate singular maps
$\Delta_{L}:\Sigma_{L}\to\Pi_{L}$ and $\Delta_{M}:\Sigma_{M}\to\Pi_{M}$ can be
expressed exactly in the same way as before, namely

$\Delta_{L}(\xi_{L})=\xi_{L}^{\lambda}\quad\text{and}\quad\Delta_{M}(\xi_{M})=\xi_{M}^{\mu}.$

On the contrary, the truncated regular maps $Q:\Pi_{L}\to\Sigma_{M}$ and
$R:\Pi_{M}\to\Sigma_{L}$ become

$Q(\eta_{L})=\beta_{2}-\theta_{1}\eta_{\lambda}\quad\text{and}\quad
R(\eta_{M})=\beta_{1}-\theta_{2}\eta_{M}.$

Thus, in the non-monodromic case we obtain the model Poincaré map
$P:\Sigma_{L}\to\Sigma_{L}$ as

$P(\xi_{L})=\beta_{1}-\theta_{2}(\beta_{2}-\theta_{1}\xi_{L}^{\lambda})^{\mu}.$
(4)

Note that in terms of the model maps we can change between the monodromic and
non-monodromic cases through the transformation
$(\theta_{1},\theta_{2})\mapsto(-\theta_{1},-\theta_{2})$.

It has been shown by Dukov [2018] that taking into account higher-order terms
in the singular and regular maps in (4) does not change topology of the
bifurcation diagrams in both subcases.

##  4 Two polynomial examples

Here we present two systems of polynomial ordinary differential equations
which illustrate some of the theoretical results described in the previous
section. These explicit examples allow for numerical computation of the
various phase objects and bifurcation diagrams using matcont Dhooge _et al._
[2003, 2008]; De Witte _et al._ [2012], which can be compared to their
theoretically predicted counterparts.

Our strategy to find explicit monodromic examples – originally applied by
Sandstede [1997] to construct homoclinic bifurcation examples – will be to
look for an algebraic variety which has the shape of a monodromic contour and
then to look for vector fields which are tangent to this variety. For a non-
monodromic example a different construction is applied. It is based on
perturbing a reversible system containing a big homoclinic connection.

###  4.1 A monodromic example

####  4.1.1 Derivation of the example

To construct a vector field containing a monodromic heteroclinic contour, we
will look for an algebraic variety which has a shape similar to such a
contour. For this example we choose the union of a parabola with the $x$-axis.
Specifically, we choose the parabola described by $y=x(1-x)$, see Fig. 8.

01$x$$y$ Figure 8: The union of the parabola $y=x(1-x)$ with the $x$-axis.

This algebraic variety can be realised as the zero-set of the function

$G(x,y)=y(y-x(1-x))$ (5)

and we write $\Gamma=G^{-1}(0)$.

Our next step is to look for two polynomial functions $f(x,y)$ and $g(x,y)$
such that the vector field

$\begin{cases}\dot{x}=f(x,y)\\\ \dot{y}=g(x,y)\end{cases}$ (6)

is tangent to $\Gamma$. For this to happen it must hold that

$\langle\nabla G(x,y),\begin{pmatrix}f(x,y)\\\ g(x,y)\end{pmatrix}\rangle=0$
(7)

for all $(x,y)\in\Gamma$. The equation (7) will yield a condition from which
we can derive preliminary expressions for $f(x,y)$ and $g(x,y)$. As we will
soon see, it will suffice to use polynomials of second degree. Thus, we set

$\begin{array}[]{rcl}f(x,y)&=&a_{10}x+a_{01}y+a_{20}x^{2}+a_{11}xy+a_{02}y^{2},\\\
g(x,y)&=&b_{10}x+b_{01}y+b_{20}x^{2}+b_{11}xy+b_{02}y^{2}.\end{array}$ (8)

Then, (7) together with $G(x,y)=0$ yields the system of two polynomial
equations

$\begin{array}[]{rcl}b_{10}x^{2}+(a_{10}-2b_{10}+b_{01})xy+(a_{01}-2b_{01})y^{2}&-&\\\
(b_{10}-b_{20})x^{3}&-&\\\ (2a_{10}+2b_{20}-b_{11}+b_{01}-a_{20})x^{2}y&+&\\\
(a_{11}-2b_{11}-2a_{01}+b_{02})xy^{2}+(a_{02}-2b_{02})y^{3}&-&\\\
b_{20}x^{4}-(2a_{20}+b_{11})x^{3}y-(2a_{11}+b_{02})x^{2}y^{2}&-&\\\
2a_{02}xy^{3}&=&0,\\\ y^{2}-xy+x^{2}y&=&0,\end{array}$ (9)

where the coefficients $a_{jk}$ and $b_{jk}$ are yet unknown.

To study this system we use the software package Maple to calculate its
Gröbner basis using the graded reverse lexicographic order for the
$x^{j}y^{k}$-monomials. This yields a total of four basis polynomials
including $G(x,y)$ itself. To look for solutions of (9), we choose the basis
polynomial of the minimal degree and require all its coefficients to be
identically zero. This yields the system of linear equations

$\left\\{\begin{array}[]{rcl}b_{20}&=&0,\\\ b_{10}-b_{20}&=&0,\\\
b_{10}&=&0,\\\ a_{02}&=&0,\\\ 2a_{01}+a_{11}-2a_{20}+b_{11}&=&0,\\\
a_{10}+2b_{10}+a_{20}+2b_{11}&=&0,\\\ a_{02}+2a_{11}-b_{02}&=&0,\\\
2a_{10}+a_{01}-b_{01}+a_{20}+2b_{20}&=&0,\end{array}\right.$ (10)

which is solved whenever we set

$\displaystyle a_{10}$ $\displaystyle=$ $\displaystyle-a_{20},$ $\displaystyle
b_{01}$ $\displaystyle=$ $\displaystyle a_{01}+a_{10},$ $\displaystyle b_{11}$
$\displaystyle=$ $\displaystyle-2a_{01}-2a_{10}-a_{11},$ $\displaystyle
b_{02}$ $\displaystyle=$ $\displaystyle 2a_{11},$

and all remaining coefficients to zero. With this choice, all basis
polynomials except $G$ become identically zero. We see that this leaves us
with three free parameters remaining. Introducing

$\displaystyle a$ $\displaystyle:=$ $\displaystyle a_{10},$ $\displaystyle b$
$\displaystyle:=$ $\displaystyle a_{01},$ $\displaystyle c$ $\displaystyle:=$
$\displaystyle a_{11},$

and substituting the result into $f(x,y)$ and $g(x,y)$ in (8), we now obtain
the vector field

$\begin{cases}\dot{x}=ax+by+cxy-ax^{2},\\\
\dot{y}=(a+b)y-(2a+2b+c)xy+2cy^{2},\end{cases}$ (11)

for which $\Gamma$ is invariant for any choice of $a$, $b$ and $c$. This can
now be verified directly, since the left-hand side of (7) contains $G(x,y)$ as
factor.

As expected, this vector field always has equilibria $L=(0,0)$ and $M=(1,0)$.
However, the types of these equilibria depend on the values of $a$, $b$ and
$c$. The eigenvalues at $(0,0)$ are $a$ and $a+b$ and the eigenvalues at
$(1,0)$ are $-a$ and $-a-b-c$.

As a first example, we set $a=1$, $b=-3$ and $c=3/2$. We find that
$\lambda_{u}=1$ and $\lambda_{s}=-2$ so
$\lambda=-\lambda_{s}/\lambda_{u}=2>1$, while $\mu_{u}=1/2$ and $\mu_{s}=-1$
so $\mu=-\mu_{s}/\mu_{u}=2>1$. Obviously, $\lambda\mu>1$ implying that the
heteroclinic contour is stable from the inside. We recognize that these saddle
indices correspond to the first subcase treated in Section 3.1, after the time
reversal. A phase portrait of this system is given in Fig. 9.

Figure 9: Phase portrait of (11) for $a=1$, $b=-3$ and $c=3/2$ Figure 10:
Phase portrait of (11) for $a=1$, $b=-3$ and $c=1/2$

As a second example, we set $a=1$, $b=-3$ and $c=1/2$. Then it remains that
$\lambda=2>1$ but now $\mu_{u}=3/2$ and $\mu_{s}=-1$ so
$\mu=-\mu_{s}/\mu_{u}=2/3<1$ and $\lambda\mu=4/3>1$. We recognize that we are
in the situation of the second subcase treated in Section 3.1. A phase
portrait of this system is shown in Fig. 10.

In order to perform bifurcation analysis we still need to add two parameter-
dependent perturbing terms to (11) which break the heteroclinic connections
separately. As these connections can be described by the equations $y=x(1-x)$
and $y=0$ this can be achieved by choosing parameters $\alpha$ and
$\varepsilon$ and adding the perturbing terms $\alpha(y-x(1-x))$ and
$\varepsilon y$ to (11). It is clear that the first perturbation does not
affect the parabola while the second perturbation does not affect the straight
line. We add one perturbing term to $\dot{x}$ and one to $\dot{y}$. This
yields

$\left\\{\begin{array}[]{rcl}\dot{x}&=&ax+by+cxy-ax^{2}+\varepsilon y,\\\
\dot{y}&=&(a+b)y-(2a+2b+c)xy+2cy^{2}\\\
&&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+~{}\alpha(y-x(1-x)).\end{array}\right.$
(12)

####  4.1.2 Numerical analysis of the example

Using the Matlab software package matcont Dhooge _et al._ [2003, 2008]; De
Witte _et al._ [2012], we perform numerical bifurcation analysis of the system
(12) in both subcases. More specifically, we will be able to locate and
continue homoclinic connections and cycles which will allow us to compute the
bifurcation diagrams of both subcases.

In Appendix 25 we show that for the chosen values of the parameters $a$, $b$
and $c$, the splitting of the heteroclinic connections under the given
perturbations is regular. This is done by showing that the relevant Melnikov
integrals are nonvanishing.

In Section 3.1, a general sketch of these diagrams is presented (see Fig. 3).
The numerical study of these polynomial systems will allow us to compare this
general diagram based on the behaviour of model maps to an actual example of
such a diagram.

This analysis in the first case yields the bifurcation diagram displayed in
Fig. 11. This analysis in the second case yields the bifurcation diagram
displayed in Fig. 12.

To complete our analysis of these systems we compute representative phase
portraits for parameter values in all components of the bifurcation diagrams.
The contents of these phase portraits were already predicted in the previous
part. They are displayed in Figs. 13 and 14.

Figure 11: The parameter plane diagram in the first subcase ($a=1$, $b=-3$ and
$c=3/2$) Figure 12: The parameter plane diagram in the second subcase ($a=1$,
$b=-3$ and $c=1/2$)

Figure 13: The full bifurcation diagram for the first subcase Figure 14: The
full bifurcation diagram for the second subcase

###  4.2 A non-monodromic example

####  4.2.1 Finding an example

The starting point is the following reversible system111The system is
invariant under the involution $(t,x,y)\mapsto(-t,-x,y)$.

$\left\\{\begin{array}[]{rcl}\dot{x}&=&y,\\\
\dot{y}&=&x+xy-x^{3},\end{array}\right.$ (13)

that has a saddle $(0,0)$, a stable focus $(-1,0)$, and an unstable focus
$(1,0)$. The saddle at the origin has a ‘big’ homoclinic orbit, while both
foci are located inside the domain bounded by the homoclinic loop as seen in
Fig. 15.

Figure 15: System (13) containing a big homoclinic loop Figure 16: The
symmetric non-monodromic heteroclinic contour in the reversible system (14)
with $\gamma=\gamma_{0}$. Figure 17: Partial bifurcation diagram of (15) with
$\gamma=2.7$: heteroclinic bifurcation curves (green), primary homoclinic
curves (blue).

The one-parameter modification of (13) that preserves the reversibility,

$\left\\{\begin{array}[]{rcl}\dot{x}&=&y(y+\gamma),\\\
\dot{y}&=&x+xy-x^{3},\end{array}\right.$ (14)

has an additional saddle at $(0,-\gamma)$. This saddle is connected to the
saddle at $(0,0)$ by two heteroclinic orbits when $\gamma_{0}=2.5315\ldots$
(see Fig. 16). Of course, both saddles of (14) are neutral.

To break reversibility, introduce the following two-parameter deformation of
(14):

$\left\\{\begin{array}[]{rcl}\dot{x}&=&\varepsilon x+y(y+\gamma),\\\
\dot{y}&=&\alpha y+x+xy-x^{3}.\end{array}\right.$ (15)

This system still has a saddle at $(0,0)$, while the second saddle moves off
the $y$-axis. In order to obtain a generic non-monodromic heteroclinic
contour, all three parameters of (15) have to vary. Since the existence of
such heteroclinic contour is a codim 2 phenomenon, we expect a curve in the
$(\alpha,\varepsilon,\gamma)$-space passing through $(0,0,\gamma_{0})$, along
which the contour persists.

####  4.2.2 Numerical analysis of the example

The main curves of the bifurcation diagram of (15) in the
$(\alpha,\varepsilon)$-plane at $\gamma=2.7$ is shown in Fig. 17. Since the
transformation $(t,x,\alpha,\varepsilon)\mapsto(-t,-x,-\alpha,-\varepsilon)$
preserves the system (15), its bifurcation diagram is invariant under the
inversion $(\alpha,\varepsilon)\mapsto(-\alpha,-\varepsilon)$.

Figure 18: Partial bifurcation diagram of (15) near the lower codim 2 point
$C_{1}$: primary and secondary heteroclinic bifurcation curves $H_{L,M}^{(k)}$
(green), primary homoclinic curves $P_{L,M}$(blue).

There are two non-monodromic heteroclinic contours at the intersection points
$C_{1,2}$ of the heteroclinic curves $H_{L,M}^{(0)}$. These points are the
end-points of two primary homoclinic curves $P_{L,M}$, which also intersect at
$(\alpha,\varepsilon)=(0,0)$, where two independent homoclinic orbits to both
saddles exist. Since the system is reversible at $\alpha=\varepsilon=0$, both
homoclinic orbits are symmetric. All orbits in the annulus bounded by these
two homoclinic orbits are periodic.

The lower heteroclinic contour point is found to be

$C_{1}=(\alpha,\varepsilon)=(0.422432\ldots,-0.452007\ldots).$

The corresponding phase portrait is shown in Fig. 19. The saddles at these
parameter values are

$L=(0,0)~{}~{}~{}{\rm and}~{}~{}~{}M=(-0.57039\ldots,-2.6009\ldots),$

with eigenvalues

$\lambda_{s}=-1.7151\ldots,\ \ \lambda_{u}=1.6856\ldots$

and

$\mu_{s}=-2.8436\ldots,\ \ \mu_{u}=2.2436\ldots,$

respectively. The corresponding indices are

$\lambda=-\frac{\lambda_{s}}{\lambda_{u}}=1.0175\ldots>1$

and

$\mu=-\frac{\mu_{s}}{\mu_{u}}=1.2674\ldots>1.$

Notice that both saddles in this heteroclinic contour are dissipative.

Figure 19: Non-monodromic heteroclinic contour in system (14) with
$\gamma=2.7$ at codim 2 point $C_{1}$.

The upper heteroclinic contour point

$C_{2}=(\alpha,\varepsilon)=(-0.422432\ldots,0.452007\ldots).$

The corresponding phase portrait is shown in Fig. 20. The saddles at these
parameter values are

$L=(0,0)~{}~{}~{}{\rm and}~{}~{}~{}M=(0.570399\ldots,-2.6009\ldots),$

with the eigenvalues

$\lambda_{s}=-1.6856\ldots,\ \ \lambda_{u}=1.7151\ldots$

and

$\mu_{s}=-2.2436\ldots,\ \ \mu_{u}=2.8436\ldots,$

respectively. The corresponding indices are

$\lambda=-\frac{\lambda_{s}}{\lambda_{u}}=0.982\ldots<1$

and

$\mu=-\frac{\mu_{s}}{\mu_{u}}=0.789\ldots<1.$

Figure 20: Non-monodromic heteroclinic contour in system (15) with
$\gamma=2.7$ at codim 2 point $C_{2}$.

Figure 21: Phase portraits for critical values of the parameters

##  5 Discussion

The theoretical results collected in this paper can be derived from the study
of the truncated model maps. However, it must be shown that adding higher
order terms to the truncated model maps does not change the results in generic
cases, so that the leading terms are indeed sufficient to obtain all relevant
theoretical results. Both papers Reyn [1980]; Dukov [2018], on which our
presentation is based, address this issue, however a unified treatment is
still absent.

There are two approaches leading to the same truncated Poincaré maps that
differ by construction of singular correspondence maps near saddles. The first
approach is based on fact that there is a local invertible $C^{1}$ coordinate
transformation that maps solutions of a nonlinear planar system near a
hyperbolic saddle onto the solutions of its linearization, see Deng [1989] or
Kuznetsov [2004] for a more geometric treatment. This allows to compute the
singular map exactly as $\xi\mapsto\xi^{\lambda}$, where $\xi$ is a coordinate
in the unstable direction and $\lambda$ is the saddle index. Since for the
analysis of cycle bifurcations further smoothness of the Poincaré map is
required, one needs to impose more genericity conditions. Actually, there is a
$C^{k}$ linearizing transformation with any $k>1$, if a finite number $N(k)$
of lower-order resonance conditions on the eigenvalues are not satisfied
Sternberg [1958]. Thus the singular map is generically as smooth as required
for $\xi>0$, provided the system (1) is sufficiently smooth.

There is an alternative approach Shilnikov _et al._ [2001] that leads to the
same truncated Poincaré maps but avoids a reduction to a linear system near
the saddles and, thus, does not impose extra conditions. In this approach, one
first introduces new coordinates near the saddle, such that the system remains
nonlinear but the stable and unstable manifolds locally coincide with the axes
(and some other simplifications are done). Then by considering a boundary-
value problem, one proves that the singular map near the saddle has the form
$\xi\mapsto\xi^{\lambda}+o(\xi^{\lambda})$. Composing such maps with the
regular global maps, gives the same truncated Poincaré maps as in the first
approach. This construction leads to a small finite loss of smoothness
relative to that of the original planar system (1). Both papers Reyn [1980];
Dukov [2018] use the first approach (linearization), while the second one
seems to be less restrictive.

In order to complete our collection of examples demonstrating the theoretical
results of this thesis, one more example needs to be added. Of the non-
monodromic contours, we were only able to produce an example of the simpler
subcase (with $\lambda<1$ and $\mu<1$). Addition of a non-monodromic example
of the other subcase (with $\lambda<1$ but $\mu>1$), displaying a fold of
cycles bifurcation will complete the stock of examples.

The bifurcations discussed in this paper all have a semi-local nature. That
is, they appear in a neighborhood of a heteroclinic contour. It is often the
case that bifurcations which have a global or semi-local nature in systems
with few parameters appear locally in systems with more parameters. For
example, in Bazykin _et al._ [1989]; Dumortier _et al._ [1991] the monodromic
case is shown to appear close to one case of a degenerate Bogdanov-Takens
bifurcation of codim 3. For the non-monodromic case, an analogous result
remains to be found. We expect that it must have codimension bigger than 3 and
that the linear part of the associated critical normal form is identically
zero.

Flashing heteroclinic connections The most exotic bifurcation phenomenon
encountered in the paper is the appearance of flashing (or sparkling)
heteroclinic connections. In fact, its appearance is one of the most important
reasons to study bifurcations of non-monodromic heteroclinic contours, as it
results in a bifurcation diagram which is surprisingly complex for a system
which is only planar. Here we summarize for the reader’s convenience the major
facts on these bifurcations, closely following Dukov [2018].

As noted in the introduction, this bifurcation was discovered in 1981 by I.P.
Malta and J. Palis in Malta & Palis [1981], see also Shilnikov _et al._
[2001]. The flashing heteroclinic connections refer to the generation of an
infinite series of heteroclinic connections between two saddles which make an
increasing number of turns around one of the saddles before reaching the
other.

The situation in which this bifurcation was first discovered was associated
with a non-hyperbolic limit cycle. Of the two saddles involved in this
bifurcation, one is located inside the limit cycle while the other is located
on the outside222Of course, more equilibria must be encircled by the cycle..
The non-hyperbolic limit cycle serves as the limit set of a solution leaving
the outside saddle for increasing $t$ while simultaneously it is approached
for decreasing $t$ by a solution approaching the other saddle for increasing
$t$. Thus, there are solutions winding around the cycle from both the outside
and inside.

The following theorem shows that when this limit cycle is destroyed, an
infinite series of heteroclinic connections is generated which wind around the
inside saddle an increasing number of times (see Fig. 23).

In the theorem we consider the parameter $\alpha\in\mathbb{R}$ involved in the
fold of cycles bifurcation. Thus, for $\alpha=0$ there exists a semi-stable
limit cycle as described above and for $\alpha>0$ this cycle is destroyed.

######  Theorem .1.

Consider a vector field containing two saddles and a semi-stable limit cycle
as described above. Then there exists a sequence
$(\alpha_{n})_{n\in\mathbb{N}}$ with $\alpha_{n}>0$ for all $n\in\mathbb{N}$
and $\alpha_{n}\to 0$ such that for every $n\in\mathbb{N}$, the vector field
corresponding to the parameter value $\alpha_{n}$ contains a heteroclinic
connection making more turns around the inner saddle then the previous one.

######  Proof .2.

For $\alpha=0$, the vector field contains a semi-stable limit cycle which
separates the saddles $L$ and $M$. On this cycle we select a cross-section
$\Sigma$ on which we define a coordinate $x$ directed outwards (see Fig. 22).
Suppose that the limit cycle and its cross-section intersect at $x=0$. The
orbit moving towards the saddle $L$ reaches the cross-section at $x_{1}<0$ and
the orbit leaving from $M$ at $x_{2}>0$.

Figure 22: Placing a cross-section on the cycle

When we pick $\alpha>0$ and study the vector field still furnished with the
cross-section $\Sigma$, we see that we get parameter dependent values of
$x_{1}(\alpha)<0<x_{2}(\alpha)$ and that the semi-stable limit cycle is
destroyed. Thus, the parameter-dependent Poincaré map $P_{\alpha}$ has no
fixed points. Therefore, we have that $P_{\alpha}(x)<x$ so after enough
iterations of the Poincaré map we will find $n\in\mathbb{N}$ and
$m\in\mathbb{N}$ such that $P_{\alpha}^{(-n)}(x_{1}(\alpha))>0$ and
$P_{\alpha}^{(m)}(x_{2}(\alpha))<0$.

In contrast to this, we always have that $P_{\alpha}(x_{1}(0))<0$ and
$P_{\alpha}(x_{2}(0))>0$. Thus, assuming that $n\geq m$, and letting
$\alpha\to 0$ we will find $\alpha_{n}$ such that
$P_{\alpha_{n}}^{(-n)}(x_{1}(\alpha_{n}))=P_{\alpha_{n}}^{(n)}(x_{2}(\alpha_{n}))=0$
which implies a heteroclinic connection.

(a) No turns (b) One turn

(c) Two turns

Figure 23: Heteroclinic connections making an increasing number of turns

The occurrence of flashing heteroclinic connections is not limited to
situations involving a cyclic fold bifurcation. They may also appear in a
similar situation involving the bifurcation of a homoclinic connection. More
specifically, we consider the situation in which there exist two saddles and a
homoclinic connection at one of the saddles. The other saddle is located
inside this saddle. We consider a system depending on the splitting parameter
$\beta$ of the homoclinic connection. Suppose the saddle at which the
homoclinic connection exists has saddle index smaller than one, so for
$\beta>0$ the homoclinic connection is broken and no limit cycle is generated.

######  Theorem .3.

Consider a vector field containing two saddles and a homoclinic connection as
described above. Then there exists a sequence $(\beta_{n})_{n\in\mathbb{N}}$
with $\beta_{n}>0$ for all $n\in\mathbb{N}$ and $\beta_{n}\to 0$ such that for
every $n\in\mathbb{N}$, the vector field corresponding to the parameter values
$\beta_{n}$ contains a heteroclinic connection making more turns around the
inner saddle then the previous one.

######  Proof .4.

The vector field for $\beta=0$ contains a homoclinic connection at $M$ and $L$
lies on the inside of this connection. The solutions leaving $L$ have the
homoclinic connection as its limit set. On this connection we define a cross-
section $\Sigma$ on which we define a coordinate $x$ directed inwards (see
Fig. 24). Suppose that the homoclinic connection (that is, the solution
leaving and approaching $M$) and its cross-section intersect at $x=0$ and that
the solution leaving the saddle $L$ reaches the cross-section at $x_{1}>0$.

Figure 24: Placing a cross-section on the homoclinic connection

When we pick $\beta>0$ and study the vector field, still furnished with the
cross-section $\Sigma$, we see that we get parameter dependent values of
$x_{1}(\beta)>0$ and that the homoclinic connection is destroyed. Moreover,
the parameter dependent Poincaré map $P_{\beta}$ has no fixed points.
Therefore, we have that $P_{\beta}(x)<x$ so after enough iterations of the
Poincaré map we will find $n\in\mathbb{N}$ such that
$P_{\beta}^{(n)}(x_{1}(\beta))<0$.

In contrast to this, we always have that $P_{\beta}(x_{1}(0))>0$. Thus,
letting $\beta\to 0$ we will find $\beta_{n}$ such that
$P_{\beta_{n}}^{(n)}(x_{1}(\beta_{n}))=0$ which corresponds to a heteroclinic
connection.

(a) No turns (b) One turn

(c) Two turns

Figure 25: Heteroclinic connections making an increasing number of turns

Regularity of splitting In this appendix we show that the splitting of the
heteroclinic connections in the planar system given by (12) is regular by
verifying that the relevant Melnikov integrals are nonzero when evaluated
along the heteroclinic connections. An introduction to this technique can be
found, e.g., in Kuznetsov [2004].

Writing a smooth planar system dependent on a parameter $p\in\mathbb{R}$ as

$\left\\{\begin{array}[]{rcl}\dot{x}&=&f(x,y,p),\\\
\dot{y}&=&g(x,y,p),\end{array}\right.$ (16)

the Melnikov integral becomes

$M_{p}(0)=\int_{-\infty}^{\infty}\varphi(t)\psi(t)dt\,,$

where

$\varphi(t)=\exp\left[-\int_{0}^{t}\left(\frac{\partial f}{\partial
x}+\frac{\partial g}{\partial y}\right)d\tau\right]$

and

$\psi(t)=f\frac{\partial f}{\partial p}-g\frac{\partial f}{\partial p}$

have to be evaluated along the heteroclinic solutions of the system at the
critical parameter value $p=0$.

The heteroclinic connection in (12) along $y=0$ splits under variation of the
parameter $\alpha$. Setting $p=\alpha$ and $y=0$ we find that

$\psi(t)=-x(t)^{2}(1-x(t))^{2}<0$

for all $t\in\mathbb{R}$. As $\varphi(t)>0$ for all $t$ as well, we conclude
that $M_{\alpha}(0)\neq 0$ independent on the parameters $a$, $b$ and $c$ of
(12).

The connection along $y=x(1-x)$ splits under variation of the parameter
$\varepsilon$ so we set $p=\varepsilon$ and $y=x(1-x)$. In this case, we
transform the integrals with respect to $t$ to integrals in $x$ over the
domain $[0,1]$, and proceed numerically. We explicitly work out the case
$c=1/2$. The remaining case is similar. Using the first equation in (12), we
obtain that

$\psi(t(x))\,dx=x(-2x^{2}+3x-1)\,dx.$

Similarly, the integrand involved in $\varphi(t(x))$ becomes

$\frac{5x^{2}-8x+2}{x^{3}-5x^{2}+4x}\,dx.$

Choosing the particular solution satisfying $x(0)=1/2$ we obtain

$\varphi(t(x))=\exp\left[-\int_{1/2}^{x}\frac{5\xi^{2}-8\xi+2}{\xi^{3}-5\xi^{2}+4\xi}\,d\xi\right].$

As $\lim_{t\to\infty}x(t)=0$, the total integral becomes

$M_{\varepsilon}(0)=-\int_{0}^{1}\varphi(t(x))\psi(t(x))dx$

Numerical evaluation yields $M_{\varepsilon}(0)=-0.0177888$. Analogously, we
find that $M_{\varepsilon}(0)=-0.0866532$ when $c=3/2$, showing that the
splitting happens regularly in both cases.

## References

  * Arnold _et al._ [1994] Arnold, V. I., Afrajmovich, V. S., Il’yashenko, Y. S. & Shil’nikov, L. P. [1994] _Dynamical Systems V. Bifurcation theory and catastrophe theory_ , Encyclopaedia of Mathematical Sciences, Vol. 5 (Springer-Verlag, Berlin), [A translation of Current Problems in Mathematics. Fundamental Directions. Vol. 5, USSR Academy of Sciences, VINITI, Moscow, 1986, in Russian].
  * Bazykin _et al._ [1989] Bazykin, A. D., Kuznetsov, Yu. A. & Khibnik, A. I. [1989] _Portraits of Bifurcations: Bifurcation diagrams of dynamical systems on the plane_ , Current Life, Science and Technology: Series “Mathematics and Cybernetics”, Vol. 89 (“Znanie”, Moscow), [in Russian].
  * De Witte _et al._ [2012] De Witte, V., Govaerts, W., Kuznetsov, Yu. A. & Friedman, M. [2012] “Interactive initialization and continuation of homoclinic and heteroclinic orbits in MATLAB.” _ACM Trans. Math. Softw._ 38, 34.
  * Deng [1989] Deng, B. [1989] “The Šil’nikov problem, exponential expansion, strong $\lambda$-lemma, $C^{1}$-linearization, and homoclinic bifurcation,” _J. Differential Equations_ 79, 189–231.
  * Dhooge _et al._ [2003] Dhooge, A., Govaerts, W. & Kuznetsov, Yu. A. [2003] “MATCONT: A MATLAB package for numerical bifurcation analysis of ODEs.” _ACM Trans. Math. Softw._ 29, 141–164.
  * Dhooge _et al._ [2008] Dhooge, A., Govaerts, W., Kuznetsov, Yu. A., Meijer, H. G. & Sautois, B. [2008] “New features of the software matcont for bifurcation analysis of dynamical systems,” _Mathematical and Computer Modelling of Dynamical Systems_ 14, 147–175.
  * Dukov [2018] Dukov, A. [2018] “Bifurcations of the ‘heart’ polycycle in generic 2-parameter families,” _Transactions of the Moscow Mathematical Society_ 79, 1.
  * Dumortier _et al._ [1991] Dumortier, F., Roussarie, R., Sotomayor, J. & Żoładek, H. [1991] “Bifurcations of planar vector fields. Nilpotent singularities and Abelian integrals.” _Lect. Notes Math._ 1480, viii + 226.
  * Kuznetsov [2004] Kuznetsov, Yu. A. [2004] _Elements of Applied Bifurcation Theory_ , 3rd ed. (Springer-Verlag, New York).
  * Malta & Palis [1981] Malta, I. P. & Palis, J. [1981] “Families of vector fields with finite modulus of stability,” _Dynamical Systems and Turbulence, Warwick 1980_ , eds. Rand, D. & Young, L.-S. (Springer Berlin Heidelberg, Berlin, Heidelberg), ISBN 978-3-540-38945-3, pp. 212–229.
  * Reyn [1980] Reyn, J. W. [1980] “Generation of limit cycles from separatrix polygons in the phase plane,” _Geometrical Approaches to Differential Equations_ , ed. Martini, R. (Springer Berlin Heidelberg, Berlin, Heidelberg), pp. 264–289.
  * Sandstede [1997] Sandstede, B. [1997] “Constructing dynamical systems having homoclinic bifurcation points of codimension two,” _Journal of Dynamics and Differential Equations_ 9, 269–288.
  * Shilnikov _et al._ [2001] Shilnikov, L. P., Shilnikov, A. L., Turaev, D. & Chua, L. O. [2001] _Methods of Qualitative Theory in Nonlinear Dynamics. Part II_ , World Scientific Series on Nonlinear Science. Series A: Monographs and Treatises, Vol. 5 (World Scientific Publishing Co., Inc., River Edge, NJ).
  * Sternberg [1958] Sternberg, S. [1958] “On the structure of local homeomorphisms of euclidean $n$-space. II,” _Amer. J. Math._ 80, 623–631.


# Tunable magnetic alignment between trapped exciton-polariton condensates

H. Ohadi<EMAIL_ADDRESS>Y. del Valle-Inclan Redondo  A. Dreismann
Department of Physics, Cavendish Laboratory, University of Cambridge,
Cambridge CB3 0HE, United Kingdom  Y. G. Rubo  Instituto de Energías
Renovables, Universidad Nacional Autónoma de México, Temixco, Morelos, 62580,
Mexico  F. Pinsker  Clarendon Laboratory, University of Oxford, Parks Road,
Oxford OX1 3PU, United Kingdom  S. I. Tsintzos  Foundation for Research and
Technology–Hellas, Institute of Electronic Structure and Laser, 71110
Heraklion, Crete, Greece  Z. Hatzopoulos  Foundation for Research and
Technology–Hellas, Institute of Electronic Structure and Laser, 71110
Heraklion, Crete, Greece  Department of Physics, University of Crete, 71003
Heraklion, Crete, Greece  P. G. Savvidis  Foundation for Research and
Technology–Hellas, Institute of Electronic Structure and Laser, 71110
Heraklion, Crete, Greece  Department of Materials Science and Technology,
University of Crete, 71003 Heraklion, Crete, Greece  J. J. Baumberg
<EMAIL_ADDRESS>Department of Physics, Cavendish Laboratory, University of
Cambridge, Cambridge CB3 0HE, United Kingdom

###### Abstract

Tunable spin correlations are found to arise between two neighboring trapped
exciton-polariton condensates which spin-polarize spontaneously. We observe a
crossover from an antiferromagnetic- to a ferromagnetic pair state by reducing
the coupling barrier in real-time using control of the imprinted pattern of
pump light. Fast optical switching of both condensates is then achieved by
resonantly but weakly triggering only a single condensate. These effects can
be explained as the competition between spin bifurcations and spin-preserving
Josephson coupling between the two condensates, and open the way to polariton
Bose-Hubbard ladders.

The development of spin-charge lattice models for understanding strongly-
correlated states of matter is a successful theme of modern quantum physics.
This has driven the desire to model and probe complex condensed matter
phenomena using highly-controlled systems, such as ultracold atoms Bloch  _et
al._ (2008), photons Szameit and Nolte (2010); Hafezi  _et al._ (2013), or
superconducting junctions You and Nori (2011). Exciton-polariton (polariton)
lattices have emerged as an alternative system Carusotto and Ciuti (2013);
Jacqmin  _et al._ (2014) with unique properties. Due to their strongly
dissipative and nonlinear nature, many-body polariton gases can reach steady
states which are remarkably different from their equilibrium case Le Boité
_et al._ (2013). Moreover, they have peculiar spin properties Leyder  _et al._
(2007); Lagoudakis  _et al._ (2009); Sala  _et al._ (2015) and exhibit
spontaneous magnetization (emitting circularly polarized light) above a
critical bifurcation threshold Ohadi  _et al._ (2015), analogous to the weak
lasing regime Aleiner  _et al._ (2012). In this Letter we study the basic
building block of a polariton spin lattice: two optically trapped spin-
polarized condensates which are tunably coupled. We demonstrate that trapped
out-of-equilibrium polariton condensates can exhibit Ising-like behavior
related to spin bifurcations. The two condensate system investigated here is
shown to correspond to one plaquette of a bosonic ladder Atala  _et al._
(2014), and allows demonstration of a crossover in the competition between
Josephson coupling and spin bifurcation. These features have not been seen in
any other system to date.

Polariton condensates are coherent many-body states Deng  _et al._ (2002);
Kasprzak  _et al._ (2006); Balili  _et al._ (2007); Baumberg  _et al._ (2008),
which can be confined in potentials Wertz  _et al._ (2010); Tosi  _et al._
(2012a); Cristofolini  _et al._ (2013) and interact with each other via
Josephson junctions Lai  _et al._ (2007); Lagoudakis  _et al._ (2010); Tosi
_et al._ (2012b); Abbarchi  _et al._ (2013). For a pair of interacting trapped
spin-polarized condensates, their polarization states are expected to couple.
However, the driven-dissipative and nonlinear nature of polariton condensates
makes the underlying coupling mechanism considerably richer than that in the
conventional Ising case, leading to exotic forms of magnetism where the
orientation and strength of coupling is not determined by the sign of the
interaction.

Figure 1: (a) Coupling between condensates 1,2 controlled by Josephson
tunneling $J$, spin-coupling ($+\leftrightarrow-$) within condensates
controlled by energy splitting $\varepsilon$. (b) Schematic double condensate
trap from 10 pump beams. (c,d) Experimental spin states seen for (c) AF- and
(d) F-coupled condensates. In each case two possible states exist, with actual
state chosen randomly upon each realization.

We achieve tuning between ferromagnetic (F) or antiferromagnetic (AF) coupling
by directly modulating the tunneling barrier, adjusting either the height of
the barrier or the separation between the condensates. We show how optical
switching of the spin state of one condensate results in fast switching of the
state of the neighboring condensate. Our result is a key step towards using
trapped polariton condensates for the realization of interacting bosons in a
driven-dissipative spinor Bose-Hubbard model Le Boité  _et al._ (2013);
Carusotto and Ciuti (2013).

The Bose-Hubbard model with polaritons can be studied in other systems where,
for example, the sample is etched to form micro-pillar arrays Abbarchi  _et
al._ (2013); Jacqmin  _et al._ (2014); Sala  _et al._ (2015), metal films are
deposited on the surface of the cavity Lai  _et al._ (2007); Kim  _et al._
(2008, 2011), or surface acoustic waves are applied Cerda-Méndez  _et al._
(2010). While the confinement potential is then separate from the condensate
gain, the trapped exciton-polaritons studied in the present letter have the
advantage that the confinement, even at a single site level, is versatile and
can be adjusted on the fly. Particularly for larger arrays of condensates,
this capability to tune the different barriers is vital. Moreover, the
interaction with the reservoir particles is reduced, crucial for spin
stability.

Exciton-polaritons (polaritons) are quasiparticles formed by the strong
coupling of excitons in semiconductor quantum wells with photons in the
microcavity in which they are embedded Kavokin  _et al._ (2007). We create
optically trapped polariton condensates Cristofolini  _et al._ (2013);
Askitopoulos  _et al._ (2013); Dreismann  _et al._ (2014) by nonresonant
linearly-polarized continuous wave (CW) excitation of a membrane microcavity
(for details see SI. 1). Driven by their repulsive excitonic interactions,
polaritons travel away from the pump region and feed the zero-momentum ground
state at the center of the optical trap. Once the density exceeds the
condensation threshold, a macroscopically coherent condensate forms in the
trap center. Polaritons in quantum-well microcavities have two $\pm 1$ (spin
up or down) projections of their total angular momentum along the structure
growth axis, which correspond to right- and left-circularly polarized photons
emitted by the cavity (of intensities $i_{\pm}$).

Trapped polariton condensates spontaneously exhibit a high degree of circular
polarization, or magnetization $M=s_{z}=(i_{+}-i_{-})/(i_{+}+i_{-})$, above a
critical spin-bifurcation threshold as a result of energy and dissipation
splitting of their linear polarizations Ohadi  _et al._ (2015). We operate
above this spin-bifurcation threshold, which means that the trapped condensate
is spin-polarized in either spin-up $\ket{\uparrow}$ or spin-down
$\ket{\downarrow}$ states. These spin-polarized condensates emit nearly
circularly polarized ($|M|>$85%) light, which can be measured using
conventional polarimetry. In each realization we excite the sample for
$200\text{\,}\mathrm{\SIUnitSymbolMicro s}$ and measure the condensate
$M_{n}(t)$, where $n\in\\{1,2\\}$ denotes the left and right condensate (see
SI. 2). The optical excitation is patterned using a spatial light modulator
(SLM) into the shape of a double-hexagon as shown by dashed circles in Fig.
1(b-d) such that a spin-polarized condensate is formed at the center of each
hexagon. The middle ‘barrier’ pump spots [orange in Fig. 1(b)] between the two
traps are weaker than the outer spots (intensity ratio $\simeq 75\%$) to allow
inter-trap tunneling of polaritons.

We can spatially squeeze or stretch the traps and change the condensates
separation without changing the barrier pump intensity. To achieve this the
barrier pump spots of the double-hexagon trap are fixed while shifting the
location of the other spots. When the separation of the condensates maxima is
greater than $\ell_{c}=$13.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, the two
condensates independently pick a spin-up or spin-down state. However, when
this is decreased to $\ell=0.90\ell_{c}$, we observe AF coupling, where the
condensates spontaneously collapse into either $\ket{\uparrow\downarrow}$ or
$\ket{\downarrow\uparrow}$ states in each realization [Fig. 1(c), see also SI.
2]. Further decreasing the separation to $0.74\ell_{c}$, we observe F coupling
where the condensates pick either of $\ket{\uparrow\uparrow}$ or
$\ket{\downarrow\downarrow}$ states randomly in each realization [Fig. 1(d)].
As in the case of a single condensate Ohadi  _et al._ (2015), these states
remain stable for many seconds at $5\text{\,}\mathrm{K}$, and do not depend on
the position on the sample, the geometrical pattern of the pump spots or the
power above the spin-bifurcation threshold.

Figure 2: 2D histogram of correlations in measured double condensate spins
$(M_{1},M_{2})$ over $1000$ realizations in the (a) antiferromagnetic and (b)
ferromagnetic coupling regimes.

For each set of trapping conditions, $1000$ realizations are created and for
each we measure $M$ for the two condensates to perform a statistical analysis
of the condensate-pair spin correlation. Each polarization-resolved
realization is recorded for $200\text{\,}\mathrm{\SIUnitSymbolMicro s}$ by a
camera, allowing us to map the 2D histogram of $(M_{1},M_{2})$ and resolving F
and AF situations (Fig. 2). Absolute correlations of $|C|>0.99$ are found for
the condensate spins in both coupling regimes. Increasing the condensate
separation to $\ell_{c}$ reduces this correlation to 0.09, confirming that
condensate spins then become uncoupled. It is important to note that we
observe partial phase coherence between the condensates in both AF and F
regimes, and the condensates are at equal energies (see SI. 3).

Figure 3: (a,b) Snapshot in time of left and right condensate $\sigma_{+}$
intensities for (a) $u_{r}$=0.80 and (b) $u_{r}$=0.65, showing spontaneous
correlated flipping of both condensates spins in (a) AF and (b) F regimes. (c)
Magnetic phase diagram, with measured correlation $C$ between magnetization of
condensates 1,2 as a function of barrier $u_{r}$ (circles, dotted line is
guide for the eye). Simulations give the correlation $C$ as a function of
$J/\varepsilon$ (red line).

To accurately map the magnetic phase diagram of the system the influence of
barrier on condensate spin correlation is investigated. Instead of changing
the separation of the condensates by changing the trap geometry, we vary their
barrier potential. For this the intensity of the barrier pump spots [orange,
Fig. 1(a)] is changed. This allows finer control over the coupling interaction
than changing the separation of the condensates, which is discretized due to
single-pixel shifting of the SLM. To better observe the correlations, we
induce spin flips by increasing the spin noise via spatially broadening the
pump spots, which increases the overlap of the condensates and the pump. This
increases the spontaneous spin-flip rate of the condensates (here set to
$\sim$10 flips/$\mu\mathrm{s}$, see also SI. 4), allowing faster and more
reliable correlation measurements within the stability time of our setup.
Since the energy blueshift from pump-injected excitons is proportional to the
pump intensity Askitopoulos  _et al._ (2013), the ratio of the intensity of
barrier pump spots to the other spots, $u_{r}$, is a reliable measure of the
relative barrier height. The blueshift above the condensate energy at the
saddle point of the barrier is $u_{r}U_{0}$, where
$U_{0}\simeq$200\text{\,}\mathrm{\SIUnitSymbolMicro eV}$$ (see also SI. 3).
Selecting the right-circularly polarized ($\sigma_{+}$) emission, we spatially
resolve the left and right condensates recording their intensities at each
$u_{r}$ for $2\text{\,}\mathrm{ms}$ using photomultipliers. A typical trace
for $u_{r}=0.8$ [Fig. 3(a)] shows that the double condensates flip randomly
between the two AF states with a switching time limited by our measurement
resolution ($\sim 5\text{\,}\mathrm{ns}$). Reducing $u_{r}$ to 0.65 shows now
flipping between two F-states [Fig. 3(b)]. For each barrier potential we
record 20 traces (each lasting $75\text{\,}\mathrm{\SIUnitSymbolMicro s}$) and
calculate the average correlation of the condensates spins, plotting this as a
function of $u_{r}$ [Fig. 3(c)]. We observe a clear transition from the
ferromagnetic to antiferromagnetic state at $u_{r}\simeq 0.72$. Increasing
$u_{r}$ further results in zero coupling.

The uncoupling of the condensates when increasing their separation
$\ell/\ell_{c}=0.9\rightarrow 1.0$ (while the intensity of the shared pump
spots remains constant) implies that the shared reservoir between two
condensates does not play a significant role here. Our trapped condensates
form with $\bar{k}$=0 ($\Delta
k=$0.4\text{\,}\mathrm{rad}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$),
where the transverse-electric and transverse-magnetic splitting vanishes
Shelykh  _et al._ (2004). As a result the optical spin-Hall effect Kavokin
_et al._ (2005) is negligible in trapped condensates, with spin torque rates
much smaller than tunnelling rates thus preserving spin during the Josephson
process 111At
$k=$0.4\text{\,}\mathrm{rad}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$
the half-wavelength of the spin precession observed in Ref. 44 is $\sim
375\text{\,}\mathrm{\SIUnitSymbolMicro m}$, which is more than 2 orders of
magnitude larger than the condensate separations we have here.. On the other
hand, observation of phase coherence between the condensates signifies that
the spin coupling must be mediated by a coherent mechanism, as described by
our theory below.

Our description of above effects is based on the theory for a single trapped
condensate Ohadi  _et al._ (2015), which is extended to include the Josephson
coupling Wouters and Carusotto (2007); Shelykh  _et al._ (2008); Read  _et
al._ (2010) between the two condensates. The order parameter for each exciton-
polariton condensate is a two-component complex vector
$\Psi_{n}=[\psi_{n+},\psi_{n-}]^{\mathrm{T}}$ and $\psi_{n+}$ and $\psi_{n-}$
are the spin-up and spin-down wave functions. The components of the order
parameter define the measurable condensate pseudospin
$\mathbf{S}_{n}=(1/2)(\Psi_{n}^{\dagger}\cdot\bm{\sigma}\cdot\Psi_{n})$, and
the normalized spin vector $\mathbf{\hat{s}}_{n}=\mathbf{S}_{n}/S_{n}$, where
$\sigma_{x,y,z}$ are the Pauli matrices. The order parameters evolve according
to the driven dissipative equation

$\displaystyle\begin{split}\mathrm{i}\frac{d\Psi_{n}}{dt}=&-\textstyle{\frac{\mathrm{i}}{2}}g(S_{n})\Psi_{n}-\frac{\mathrm{i}}{2}(\gamma-\mathrm{i}\varepsilon)\sigma_{x}\Psi_{n}\\\
&+\textstyle{\frac{1}{2}}\left[(\alpha_{1}+\alpha_{2})S_{n}+(\alpha_{1}-\alpha_{2})S_{nz}\sigma_{z}\right]\Psi_{n}\\\
&-\textstyle{\frac{1}{2}}J\Psi_{3-n}.\end{split}$ (1)

Here $g(S_{n})=\Gamma-W+\eta S_{n}$ is the pumping-dissipation balance,
$\Gamma$ is the (average) dissipation rate, $W$ is the incoherent in-
scattering (or ‘harvest’ rate), and $\eta$ captures the gain-saturation term
Keeling and Berloff (2008). This gain saturation depends on the total
occupation of the condensate (treated more generally in Ref. Ohadi  _et al._
(2015)). $X$ (horizontal) and $Y$ (vertical) linearly-polarized single-
polariton states are split in energy by $\varepsilon$ and dissipation rate by
$\gamma$. The repulsive interaction constant for polaritons with the same spin
is $\alpha_{1}$, and the interaction constant for polaritons with opposite
spins is $\alpha_{2}$. Finally, $J>0$ is the extrinsic spin-preserving
Josephson coupling between the left and the right condensate, which we have
introduced here to account for coherent coupling.

It is important to note the fundamental differences between the spin coupling
mechanism demonstrated here and that seen in closed systems such as atoms
trapped in optical lattices Struck  _et al._ (2011); Greif  _et al._ (2013).
In equilibrium spin systems such as those described by the Ising model, the
coupling is achieved via the minimization of the total energy. In the driven-
dissipative system described here, the minimization of energy does not play a
direct role since the system is out of equilibrium. Here the spin alignment is
a direct result of spin bifurcation correlation, which itself is a product of
pumping, dissipation and nonlinearity, and coherent exchange of particles
between the two condensates.

We perform dynamical simulations to calculate spin correlations. We calculate
the steady state of the coupled Eq. 1 for $10000$ realizations with random
initial conditions 222Similar to Ref. 11, the parameters used for 0D
simulations are: $W=$0.18\text{\,}\mathrm{p}\mathrm{s}^{-1}$$;
$\Gamma=$0.1\text{\,}\mathrm{p}\mathrm{s}^{-1}$$;
$\varepsilon=$0.04\text{\,}\mathrm{p}\mathrm{s}^{-1}$$;
$\gamma=0.2\varepsilon$;
$\alpha_{1}=$0.01\text{\,}\mathrm{p}\mathrm{s}^{-1}$$;
$\alpha_{2}=-0.5\alpha_{1}$; $\eta=$0.02\text{\,}\mathrm{p}\mathrm{s}^{-1}$$
and plot the final circular polarization correlation of the two condensates at
different Josephson coupling rates [Fig. 3(c), red line]. For small Josephson
amplitudes there is AF coupling of the condensate states. The stable AF
configuration is characterized by $\psi_{2+}=-\psi_{1-}$ and
$\psi_{2-}=-\psi_{1+}$, so that Eq. (1) is reduced to the two-component
problem with renormalized splitting between $X$\- and $Y$-polarized states
$\varepsilon^{\prime}=\varepsilon-J$. Correspondingly, the AF state loses
stability at $\varepsilon=J$ when $\varepsilon^{\prime}$ changes sign. The AF
state is first converted to the limit cycle motion by the Hopf bifurcation and
then to pseudo-chaotic behavior with further increase of $J$. Simultaneously,
the stable F configuration with $\psi_{2+}=\psi_{1+}$ and
$\psi_{2-}=\psi_{1-}$ emerges for $J>\varepsilon$. This behavior is in good
agreement with the experiment: we observe the AF coupling at high barrier
heights where the tunneling rate between two condensates (Josephson coupling)
is small and same-site coupling between spins dominates [Fig. 1(a), white
arrows]. On the other hand when the barrier height is low, F-coupling is seen
since the Josephson coupling dominates [Fig. 1(a), black arrows] 333We note
that apart from most probable pure F and AF configuration, other attractors
exist for the system (1) even away from the crossover region
$J\approx\varepsilon$. In particular, the F configuration for $J>\varepsilon$
can be also achieved with the broken $1\leftrightarrow 2$ symmetry.. Each
interacting pair of condensates then forms a plaquette, and longer chains of
condensates will form a bosonic ladder Atala  _et al._ (2014), which will in
future be useful to probe with magnetic fields van der Zant  _et al._ (1992).
We note that paired condensates emit light at the same frequency, as it is
also found experimentally.

The chaotic dynamics in the AF-F crossover region of the present system is
different from the chaos in resonantly excited pair of exciton-polariton
condensates studied before Solnyshkov  _et al._ (2009). The trapped exciton-
polariton condensates here are excited _incoherently_ , yielding chaos in an
autonomous dynamical system. The details of this chaos and its effects on the
emission from the trapped condensates will be studied elsewhere.

Figure 4: (a) Spin switching for AF-coupled condensates. The
$\sigma_{+}$-polarized gate is applied on the right condensate (shown above).
The condensate-pair switches from $\ket{\uparrow\downarrow}$ to
$\ket{\downarrow\uparrow}$. (b) Spin switching for F-coupled condensates. The
condensate-pair switches from $\ket{\downarrow\downarrow}$ to
$\ket{\uparrow\uparrow}$. The top panels show the gate laser intensity profile
G . The gate laser is turned on at $t=0$.

Finally, we demonstrate the resonant switching of the coupled spin states
experimentally. We resonantly excite one of the condensates with a narrow
linewidth CW diode laser, which we refer to as the gate (G). The membrane
microcavity allows resonant excitation from the back side of the cavity
without requiring filtering of the laser backscatter. The gate laser as well
as the pump can be switched on or off by acousto-optic modulators (AOM) with a
rise time of $\sim 40\text{\,}\mathrm{ns}$. We resonantly excite the right
condensate $20\text{\,}\mathrm{\SIUnitSymbolMicro s}$ after we turn on the
pump laser. The gate is applied on the right-hand condensate and is right-
circularly polarized ($\sigma_{+}$). Applying a cross-polarized gate switches
the polarization state of the condensate, as previously shown Ohadi  _et al._
(2015). Here, we observe that resonant switching of the spin of one condensate
also switches the spin of the other coupled condensate, both in the
antiferromagnetic and ferromagnetic regimes (Fig. 4). The condensates remain
in the switched state after the gate laser is switched off, due to the
bistable nature of the spin states. We note that the resolution-limited
condensate switching time reported here is an order of magnitude shorter than
that of the gate, clearly showing that (as for single condensates Ohadi  _et
al._ (2015)) the switching process is nonadiabatic. Here we are able to switch
the coupled spin states by transiently injecting minority spins which are only
1% of the condensates majority spin. Our theoretical description reproduces
the fast spin flips observed in the experiment when noise is injected into Eq.
1. Our theory also shows that in both AF- and F-regimes, switching the spin
state of condensate 1 also switches the spin of condensate 2, both remaining
switched after the pulse is turned off (see SI. 5).

In conclusion, we demonstrate a tunable spin coupling mechanism for trapped
polariton condensates. A transition from antiferromagnetic to ferromagnetic
coupling is seen as the potential barrier between the condensates decreases.
Our results correspond well to the interplay between spin-bifurcation and
Josephson coupling in theory. Finally resonant switching of the spin states in
both ferromagnetic and antiferromagnetic regimes is shown.

_Acknowledgments_ —Authors acknowledge Andrew Ramsay, Dieter Jaksch and Tomi
Johnson for fruitful discussions. This work was supported by EPSRC Grant No.
EP/G060649/1, EP/L027151/1, EU Grant No. INDEX 289968, ERC Grant No. LINASS
320503, the Leverhulme Trust Grant No. VP1-2013-011, Spanish MEC (Grant No.
MAT2008-01555), the Greek GSRT ARISTEIA Apollo program and Fundación La Caixa,
and Mexican CONACYT Grant No. 251808. FP acknowledges financial support
through a Schrödinger Fellowship at the University of Oxford and NQIT project
(EP/M013243/1). The data corresponding to the figures in this paper can be
found at https://www.repository.cam.ac.uk/handle/1810/253679.

## References

  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Szameit and Nolte (2010) A. Szameit and S. Nolte, J. Phys. B: At. Mol. Opt. Phys. 43, 163001 (2010).
  * Hafezi  _et al._ (2013) M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor, Nat. Photon. 7, 1001 (2013).
  * You and Nori (2011) J. Q. You and F. Nori, Nature 474, 589 (2011).
  * Carusotto and Ciuti (2013) I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).
  * Jacqmin  _et al._ (2014) T. Jacqmin, I. Carusotto, I. Sagnes, M. Abbarchi, D. D. Solnyshkov, G. Malpuech, E. Galopin, A. Lemaître, J. Bloch, and A. Amo, Phys. Rev. Lett. 112, 116402 (2014).
  * Le Boité  _et al._ (2013) A. Le Boité, G. Orso, and C. Ciuti, Phys. Rev. Lett. 110, 233601 (2013).
  * Leyder  _et al._ (2007) C. Leyder, M. Romanelli, J. P. Karr, E. Giacobino, T. C. H. Liew, M. M. Glazov, A. V. Kavokin, G. Malpuech, and A. Bramati, Nat. Phys. 3, 628 (2007).
  * Lagoudakis  _et al._ (2009) K. G. Lagoudakis, T. Ostatnický, A. V. Kavokin, Y. G. Rubo, R. André, and B. Deveaud-Plédran, Science 326, 974 (2009).
  * Sala  _et al._ (2015) V. G. Sala, D. D. Solnyshkov, I. Carusotto, T. Jacqmin, A. Lemaître, H. Terças, A. Nalitov, M. Abbarchi, E. Galopin, I. Sagnes, _et al._ , Phys. Rev. X 5, 011034 (2015).
  * Ohadi  _et al._ (2015) H. Ohadi, A. Dreismann, Y. G. Rubo, F. Pinsker, Y. del Valle-Inclan Redondo, S. I. Tsintzos, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Phys. Rev. X 5, 031002 (2015).
  * Aleiner  _et al._ (2012) I. L. Aleiner, B. L. Altshuler, and Y. G. Rubo, Phys. Rev. B 85, 121301 (2012).
  * Atala  _et al._ (2014) M. Atala, M. Aidelsburger, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Nat Phys 10, 588 (2014).
  * Deng  _et al._ (2002) H. Deng, G. Weihs, C. Santori, J. Bloch, and Y. Yamamoto, Science 298, 199 (2002).
  * Kasprzak  _et al._ (2006) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, _et al._ , Nature 443, 409 (2006).
  * Balili  _et al._ (2007) R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West, Science 316, 1007 (2007).
  * Baumberg  _et al._ (2008) J. Baumberg, A. Kavokin, S. Christopoulos, A. Grundy, R. Butté, G. Christmann, D. Solnyshkov, G. Malpuech, G. Baldassarri Höger von Högersthal, E. Feltin, _et al._ , Phys. Rev. Lett. 101, 136409 (2008).
  * Wertz  _et al._ (2010) E. Wertz, L. Ferrier, D. D. Solnyshkov, R. Johne, D. Sanvitto, A. Lemaître, I. Sagnes, R. Grousson, A. V. Kavokin, P. Senellart, G. Malpuech, and J. Bloch, Nat. Phys. 6, 860 (2010).
  * Tosi  _et al._ (2012a) G. Tosi, G. Christmann, N. G. Berloff, P. Tsotsis, T. Gao, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Nat. Phys. 8, 190 (2012a).
  * Cristofolini  _et al._ (2013) P. Cristofolini, A. Dreismann, G. Christmann, G. Franchetti, N. G. Berloff, P. Tsotsis, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Phys. Rev. Lett. 110, 186403 (2013).
  * Lai  _et al._ (2007) C. W. Lai, N. Y. Kim, S. Utsunomiya, G. Roumpos, H. Deng, M. D. Fraser, T. Byrnes, P. Recher, N. Kumada, T. Fujisawa, and Y. Yamamoto, Nature 450, 529 (2007).
  * Lagoudakis  _et al._ (2010) K. G. Lagoudakis, B. Pietka, M. Wouters, R. André, and B. Deveaud-Plédran, Phys. Rev. Lett. 105, 120403 (2010).
  * Tosi  _et al._ (2012b) G. Tosi, G. Christmann, N. G. Berloff, P. Tsotsis, T. Gao, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Nat. Comm. 3, 1243 (2012b).
  * Abbarchi  _et al._ (2013) M. Abbarchi, A. Amo, V. G. Sala, D. D. Solnyshkov, H. Flayac, L. Ferrier, I. Sagnes, E. Galopin, A. Lemaître, G. Malpuech, and J. Bloch, Nat. Phys. 9, 275 (2013).
  * Kim  _et al._ (2008) N. Y. Kim, C.-W. Lai, S. Utsunomiya, G. Roumpos, M. Fraser, H. Deng, T. Byrnes, P. Recher, N. Kumada, T. Fujisawa, and Y. Yamamoto, phys. stat. sol. (b) 245, 1076 (2008).
  * Kim  _et al._ (2011) N. Y. Kim, K. Kusudo, C. Wu, N. Masumoto, A. Löffler, S. Höfling, N. Kumada, L. Worschech, A. Forchel, and Y. Yamamoto, Nat. Phys. 7, 681 (2011).
  * Cerda-Méndez  _et al._ (2010) E. A. Cerda-Méndez, D. N. Krizhanovskii, M. Wouters, R. Bradley, K. Biermann, K. Guda, R. Hey, P. V. Santos, D. Sarkar, and M. S. Skolnick, Phys. Rev. Lett. 105, 116402 (2010).
  * Kavokin  _et al._ (2007) A. V. Kavokin, J. Baumberg, G. Malpuech, and F. P. Laussy,  _Microcavities_ (Oxford Univ. Press, Oxford, 2007).
  * Askitopoulos  _et al._ (2013) A. Askitopoulos, H. Ohadi, A. V. Kavokin, Z. Hatzopoulos, P. G. Savvidis, and P. G. Lagoudakis, Phys. Rev. B 88, 041308 (2013).
  * Dreismann  _et al._ (2014) A. Dreismann, P. Cristofolini, R. Balili, G. Christmann, F. Pinsker, N. G. Berloff, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Proc. Natl. Acad. Sci. U.S.A. 111, 8770 (2014).
  * Shelykh  _et al._ (2004) I. Shelykh, K. V. Kavokin, A. V. Kavokin, G. Malpuech, P. Bigenwald, H. Deng, G. Weihs, and Y. Yamamoto, Phys. Rev. B 70, 035320 (2004).
  * Kavokin  _et al._ (2005) A. Kavokin, G. Malpuech, and M. Glazov, Phys. Rev. Lett. 95, 136601 (2005).
  * Note (1) At $k=$0.4\text{\,}\mathrm{rad}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ the half-wavelength of the spin precession observed in Ref. 44 is $\sim 375\text{\,}\mathrm{\SIUnitSymbolMicro m}$, which is more than 2 orders of magnitude larger than the condensate separations we have here.
  * Wouters and Carusotto (2007) M. Wouters and I. Carusotto, Phys. Rev. Lett. 99, 140402 (2007).
  * Shelykh  _et al._ (2008) I. A. Shelykh, D. D. Solnyshkov, G. Pavlovic, and G. Malpuech, Phys. Rev. B 78, 041302 (2008).
  * Read  _et al._ (2010) D. Read, Y. G. Rubo, and A. V. Kavokin, Phys. Rev. B 81, 235315 (2010).
  * Keeling and Berloff (2008) J. Keeling and N. G. Berloff, Phys. Rev. Lett. 100, 250401 (2008).
  * Struck  _et al._ (2011) J. Struck, C. Ölschläger, R. L. Targat, P. Soltan-Panahi, A. Eckardt, M. Lewenstein, P. Windpassinger, and K. Sengstock, Science 333, 996 (2011).
  * Greif  _et al._ (2013) D. Greif, T. Uehlinger, G. Jotzu, L. Tarruell, and T. Esslinger, Science 340, 1307 (2013).
  * Note (2) Similar to Ref. 11, the parameters used for 0D simulations are: $W=$0.18\text{\,}\mathrm{p}\mathrm{s}^{-1}$$; $\Gamma=$0.1\text{\,}\mathrm{p}\mathrm{s}^{-1}$$; $\varepsilon=$0.04\text{\,}\mathrm{p}\mathrm{s}^{-1}$$; $\gamma=0.2\varepsilon$; $\alpha_{1}=$0.01\text{\,}\mathrm{p}\mathrm{s}^{-1}$$; $\alpha_{2}=-0.5\alpha_{1}$; $\eta=$0.02\text{\,}\mathrm{p}\mathrm{s}^{-1}$$.
  * Note (3) We note that apart from most probable pure F and AF configuration, other attractors exist for the system (1) even away from the crossover region $J\approx\varepsilon$. In particular, the F configuration for $J>\varepsilon$ can be also achieved with the broken $1\leftrightarrow 2$ symmetry.
  * van der Zant  _et al._ (1992) H. S. J. van der Zant, F. C. Fritschy, W. J. Elion, L. J. Geerligs, and J. E. Mooij, Phys. Rev. Lett. 69, 2971 (1992).
  * Solnyshkov  _et al._ (2009) D. D. Solnyshkov, R. Johne, I. A. Shelykh, and G. Malpuech, Phys. Rev. B 80, 235303 (2009).
  * Kammann  _et al._ (2012) E. Kammann, T. C. H. Liew, H. Ohadi, P. Cilibrizzi, P. Tsotsis, Z. Hatzopoulos, P. G. Savvidis, A. V. Kavokin, and P. G. Lagoudakis, Phys. Rev. Lett. 109, 036404 (2012).

Supplemental Information

##  1 Experimental methods

The cavity’s top (bottom) distributed Bragg reflector (DBR) is made of 32 (35)
pairs of Al0.15Ga0.85As/AlAs layers of
$57.2\text{\,}\mathrm{nm}$/$65.4\text{\,}\mathrm{nm}$. Four sets of three
$10\text{\,}\mathrm{n}\mathrm{m}$ GaAs quantum wells (QW) separated by
$10\text{\,}\mathrm{n}\mathrm{m}$ thick layers of Al0.3Ga0.7As are placed at
the maxima of the cavity light field. The 5$\lambda$/2
($583\text{\,}\mathrm{nm}$) cavity is made of Al0.3Ga0.7As. The microcavity
sample is chemically etched from the substrate side to form
$300\text{\,}\mathrm{\SIUnitSymbolMicro m}$ diameter membranes allowing
optical access from the back of the sample for resonant excitation. Note the
phenomena reported in this paper do not depend on this etching, which is
merely to conveniently explore the resonant gating regime. The sample shows
condensation under nonresonant excitation tsotsis_lasing . The CW pump is a
single-mode Ti:Sapphire laser tuned to the first Bragg mode $\sim
100\text{\,}\mathrm{meV}$ above the condensate energy, and is linearly
polarized. In order to create short time-scale realizations, the pump is
amplitude-modulated by an AOM (rise time $\simeq$40\text{\,}\mathrm{ns}$$). A
SLM is used to spatially pattern the pump beam into a double hexagon using the
MRAF algorithm MRAF . The separation of the pump spots can only be changed
within the single-pixel shift limit of the SLM; however, the intensity of the
individual spots can be changed nearly arbitrarily. A 0.4 NA objective is used
for imaging the pattern onto the sample. The pulsed gate, which is applied on
the back side of the membrane, is a circularly-polarized CW single-mode
extended cavity diode laser, which is amplitude modulated using a second AOM.
A cooled CCD and a $0.55\text{\,}\mathrm{m}$ spectrometer is used for imaging
and energy resolving the photoluminescence (PL). PL is analyzed using a
quarter-waveplate and a Wollaston prism in front of the camera. Two
photomultiplier tubes (PMT) and a fast oscilloscope ($1\text{\,}\mathrm{GHz}$)
are used for time-resolving the polarization-resolved PL of the left and right
condensates. The CCD, the pump/gate laser and the oscilloscope are all
synchronized together using arbitrary function generators
($20\text{\,}\mathrm{MHz}$).

##  2 real-space

Figure S1: (a) Pump intensity profiles in AF(left) and F-regimes (right) (b)
Polarization-resolved real-space PL in AF (left) and F regimes (right). (c)
Intensity profiles of the total emission in F and AF regimes.

We can squeeze or stretch the trap without affecting the barrier pump spots
[Fig.S1(a)]. Using a Wollaston prism in front of the CCD, we measure both
circular polarization components simultaneously for many realizations, as
shown in Fig. S1. Here, each realization is a short-time exposure of the
sample ($t\simeq$200\text{\,}\mathrm{\SIUnitSymbolMicro s}$$). The degree of
circular polarization $M$ is given by $(i_{+}-i_{-})/(i_{+}+i_{-})$, where
$i_{{\pm}}$ is the intensity image of the right- and left-circularly polarized
emission. The emitted polarization is not exactly circular due to the spin
bifurcation mechanism which shows the stable states have a small fraction of
linear polarization ohadi_spontaneous_2015_SI . The spatial separation of the
two condensates is extracted from the line profiles in Fig. S1(c).

##  3 Phase and Energy

Figure S2: (a) Michelson interferometer. (b) Interference pattern and
extracted phase in AF (left) and F-regime (right).

To experimentally study the phase coherence between the two condensates, a
modified Michelson interferometer is used, as shown in Fig. S2. The circularly
polarized PL is 5 times expanded in one arm and is superimposed onto that from
the other arm. By doing this the PL of one condensate overlaps the PL from
both condensates. The vertical displacement of the beams from the two arms on
a final lens before the CCD results in fringes. Observation of fringes shows
that we have partial phase coherence between the two condensates ($8\%<$
1${}^{\text{st}}$ order deg of coherence $<25\%$). The phase can be extracted
by selecting the first-order diffraction of the Fourier transformed images and
transforming it back into real-space images. We observe near zero phase
difference in the circularly polarized basis of both AF and F regimes.

Phase fluctuations in exciton-polariton condensates are shown to be due to a
combination of condensate number fluctuations and the interaction of the
condensate with the excitonic reservoir love_intrinsic_2008 . Both effects
perturb the phases of the two condensates independently because the two
condensates have separate reservoirs. On the other hand, the coherent exchange
of particles via Josephson coupling between the two condensates tries to lock
their phases. The mutual coherence is then set by the interplay between these
two effects. The phase coherence is zero in the case where the condensates are
completely uncoupled and gradually increases as the Josephson interaction
increases.

Figure S3: Phase histogram for very small Josephson coupling (a), stable AF
coupling (b), during the transition (c), and stable F-coupling (d).

Our dynamical simulation also predicts a nearly zero phase difference between
the same spin components in both F and AF regimes. Fig. S3 shows a histogram
of the phase difference between same spin components of the left and right
condensates ($\phi_{1+}-\phi_{2+}$) as a function of the Josephson coupling
rate $J$. For small Josephson couplings ($J=0.025\varepsilon$), F-states with
zero phase difference appear, but their fraction rapidly decreases as $J$
increases until two stable AF-states with broken phase symmetry appear
($J=0.625\varepsilon$). At the transition point ($J=0.9\varepsilon$) the AF
states destabilize and as $J$ increases further one stable F-state appears
($J=1.375\varepsilon$).

Figure S4: (a) Energy profile of the condensates (b) Measured blueshift at a
cross section of the hexagonal trap (shown by dotted line in inset) for
$u_{r}=0.95$.

Observation of interference fringes suggest that the two condensates have
equal energies. This is confirmed by energy resolving the condensate emission.
We observe that the two condensates emit light at the same frequency. In the
experiment the energies of the two condensates are close
($<$$100\text{\,}\mathrm{\SIUnitSymbolMicro eV}$) even when they are uncoupled
because experimentally the two traps are symmetric, i.e., the pattern
intensity is the same for both traps. Then the coherent coupling between the
two traps pulls the emission frequencies together, in a similar fashion to the
injection locking of two laser oscillators stover_locking_1966 .

By measuring the blueshift of polaritons across the saddle points of the
double-hexagon trap, we can extract the potential profile, as shown in Fig.
S4. The potential barrier between the two traps is a linear function of the
middle pump spots intensity and is $\sim
200\text{\,}\mathrm{\SIUnitSymbolMicro eV}$ at $u_{r}=1$.

##  4 Pump-induced spin flips

In previous work ohadi_spontaneous_2015_SI , we demonstrated how the
interaction of the condensate with pump-induced reservoir results in the
reduction of the average condensate magnetization $M$. In that work we noted:

> “In a similar fashion to thermal noise, an overlapping reservoir can also
> induce spin noise in the condensate. The spin noise causes condensate spin
> flips, which result in the reduction of the time-averaged circular
> polarization. Theoretically, this can be studied by introducing a noise term
> similar to Eq. (9), but instead of depending on temperature, the noise
> intensity depends on the overlapping reservoir density.”

Here, we quantify the influence of the reservoir on the spin-flip rate, as
shown in Fig. S5. We measure the rate of condensate spin-flips over
$1\text{\,}\mathrm{ms}$ as a function of pump spot diameter. We observe that
the condensate spin-flip rate increases exponentially as the overlap between
the reservoir and condensate increases and plateaus, most probably due to
finite resolution of our detection system. To measure the spin correlations as
a function of trap barrier, we broaden the pump spots to
$\sim$10\text{\,}\mathrm{f}\mathrm{l}\mathrm{i}\mathrm{p}\mathrm{s}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$$,
marked by a dotted line in Fig. S5.

Figure S5: Spin-flip rate as a function of pump spot diameter. The dotted line
marks the rate where the correlated spin oscillations were measured. Black
line is a guide for the eye.

##  5 resonant excitation and spin noise

Figure S6: The time evolution of two condensates for the case when
$J=0.8\varepsilon$ (b) and $J=1.2\varepsilon$ (c) is shown. The condensate-
pairs switch their spin states when a gate pulse with an opposite spin is
applied on condensate 1. The gate pulse is a spin-up square pulse that is
turned on at $t=$4\text{\,}\mathrm{ns}$$ and lasts for
$1\text{\,}\mathrm{ns}$. Parameters as in text with $|F|/\sqrt{N}=7.5\times
10^{-3}$, where $N$ is the occupation of the condensates with $\omega_{g}$ at
the condensate energy. Figure S7: Simulated noise-induced spin flips for AF-
($J=0.8\varepsilon$) and F-regimes ($J=1.2\varepsilon$).

To investigate resonant excitation and spin noise we modify our driven
dissipative equation to

$\displaystyle\begin{split}\mathrm{i}\frac{d\Psi_{n}}{dt}=&-\textstyle{\frac{\mathrm{i}}{2}}g(S_{n})\Psi_{n}-\frac{\mathrm{i}}{2}(\gamma-\mathrm{i}\varepsilon)\sigma_{x}\Psi_{n}\\\
&+\textstyle{\frac{1}{2}}\left[(\alpha_{1}+\alpha_{2})S_{n}+(\alpha_{1}-\alpha_{2})S_{nz}\sigma_{z}\right]\Psi_{n}\\\
&-\textstyle{\frac{1}{2}}J\Psi_{3-n}\qquad\text{(Josephson coupling)}\\\
&+F_{n}(\omega_{g},t)\quad\;\;\;\text{(resonant excitation)}\\\
&+f_{n}(t).\qquad\;\quad\text{(spin noise)}\end{split}$ (S1)

Here, $F$ is a two-component function representing resonant excitation at
frequency $\omega_{g}$. To demonstrate spin-switching within our theoretical
framework, we set the resonant gate excitation $F$ to a right-circularly
polarized ($\sigma_{+}$) pulse which is applied to condensate 1 at
$t=$4\text{\,}\mathrm{ns}$$ for a duration of $1\text{\,}\mathrm{ns}$. In both
the AF-regime and the F-regime, when the spin state of condensate 1 is
switched, condensate 2 also switches spin, as shown in Fig. S6. Moreover, both
condensates stay in the switched state after the pulse is turned off.

To demonstrate noise-induced spin flips a spin noise function $f_{n}$ is added
to Eq. S1. Here, $f_{n}(t)=[f_{n+},f_{n-}]^{T}$ where $f_{\sigma}(t)$ with
$\sigma=\pm 1$ is a realization of Gaussian random processes with zero mean
$\langle f_{n\sigma}(t)\rangle=0$ and $\delta$-like two-point correlation
function:

$\displaystyle\langle f_{n\sigma}(t)f_{n\sigma^{\prime}}(t^{\prime})\rangle$
$\displaystyle=0,$ (S2) $\displaystyle\langle
f_{n\sigma}(t)f^{*}_{n\sigma^{\prime}}(t^{\prime})\rangle$
$\displaystyle=2D\delta_{n,n^{\prime}}\delta_{\sigma,\sigma^{\prime}}\delta(t-t^{\prime}),$
(S3)

where $D$ is the total noise intensity. Fig. S7 shows the influence of spin
noise on the coupled condensates. We observe that in both regimes, noise-
induced spin flips in one condensate switch the spin of the neighboring
condensate, in agreement with the experiment.

## References

  * (1) P. Tsotsis, P. S. Eldridge, T. Gao, S. I. Tsintzos, Z. Hatzopoulos, and P. G. Savvidis, N. J. Phys 14, 023060 (2012).
  * (2) M. Pasienski and B. DeMarco, Opt. Express 16, 2176 (2008).
  * (3) H. Ohadi, A. Dreismann, Y. G. Rubo, F. Pinsker, Y. del Valle-Inclan Redondo, S. I. Tsintzos, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Phys. Rev. X 5, 031002 (2015).
  * (4) A. P. D. Love, D. N. Krizhanovskii, D. M. Whittaker, R. Bouchekioua, D. Sanvitto, S. Al Rizeiqi, R. Bradley, M. S. Skolnick, P. R. Eastham, R. André, and Le Si Dang, Phys. Rev. Lett. 101, 067404 (2008)].
  * (5) H. L. Stover and W. H. Steier, Appl. Phys. Lett. 8, 91 (1966)