Tired of Misattribution, Modeling Player Fatigue in the NBA

Austin Stephen1, Matthew Yep2, Grace Fain3

Advised by Ron Yurko4, Maksim Horowitz5

University of Wyoming1, University of California Berkeley2,University of
Oklahoma3 Carnegie Mellon University4, Atlanta Hawks5

##  Abstract

The prevailing belief propagated by NBA league observers is that the workload
of the NBA season dramatically influences a player’s performance1,2,3,4,5. We
offer an analysis of cross game player fatigue that calls into question the
empirical validity of these claims. The analysis is split into three
observational studies on prior NBA seasons. First, to offer an analysis
generalized to the whole league, we conduct an examination of relative
workloads with in-game player tracking data as a proxy for exertion. Second,
to introduce a more granular perspective, we conduct a case study of the
effectiveness of load management for Kawhi Leonard. Third, to extend the
analysis to a broader set of fatigue sources, we examine the impact of
schedule features structurally imposed on teams. All three analyses indicate
the impact of cumulative player fatigue on game outcomes is minimal. As a
result, we assert in the context of the measures already taken by teams to
maintain the health of their players,6 the impact of fatigue in the NBA is
overstated.

##  Status Quo NBA

Throughout the past decade, addressing player fatigue at a schedule level has
been and continues to be a prevalent concern for both NBA teams and the league
office.

###  Density Reduction

One way the league office has reduced the burden placed on athletes is by
removing occurrences where teams play four games in five days. In the 2012-13
season, there were 76 instances of four in five, in 2016-17 such instances
were reduced to 23, and by the 2017-18 season, four games in five days became
extinct from the schedule. This was achieved by strategically starting the
season a week earlier, such that increasing the season length from 170 to 177
days allowed for a less condensed schedule7. Such efforts have also aimed to
reduce the number of three games in four days as well as back-to-back games,
however these events are still present throughout the schedule.

As of the 2020-21 NBA season, the league office also introduced a new feature
to the schedule: two-game series, where visiting teams play the home team two
times in a row at the same arena. These series are played over the course of
either two or three days. Not only did they reduce traveling and impose
restrictions to reduce COVID-19 exposure, they also decreased the amount of
one game road trips in the schedule.

However, despite the league’s efforts to design less densely packed schedules,
many coaches such as Greg Popovich and Nick Nurse have still adopted load
management strategies where they rest their key players during the regular
season. This is problematic to the NBA league office, especially when these
star-less games are nationally televised. Not only are these games less
entertaining for fans, they consequently result in substantially lower
television ratings8.

To address this, in recent years the league office has instituted a resting
policy that prohibits teams from resting healthy players in high profile,
nationally televised games. In addition, barring unusual circumstances, teams
should not rest multiple healthy players in the same game, nor should they
rest healthy players on the road. The fines for violating these resting
policies sit at a minimum of $100,0009.

###  Schedule Fairness

Many NBA teams located on the west coast and southern borders travel
substantially more miles throughout a given season compared to northeastern
teams. Consequently, these teams also cross many more time zones to get to and
from games.

Fig. 1

The map above illustrates the cumulative flight distance traveled for each NBA
team in the 2014-2015 season. The sizes of the logos are scaled using
percentiles, showing that the larger the logo, the more that team had to
travel. The Portland Trail Blazers had to fly over 62,000 miles that season,
while the Cleveland Cavaliers only traveled 35,000 miles.

This raises an important question, to both NBA teams and the league office, as
to whether or not the teams who travel more face a structural disadvantage due
to scheduling. Aggregated travel could inflict cumulative fatigue, and having
to cross time zones could also very well disrupt players’ circadian rhythms
and inhibit their performance and ability to recover.

##  Related Work

###  Defining fatigue in basketball

To maintain consistency with prior work examining fatigue in basketball10, we
define fatigue as “an acute impairment of performance that includes both an
increase in the perceived effort necessary to exert a desired force and the
eventual inability to produce this force”11. This definition is inclusive of
fatigue from any source not discriminating against muscular cognitive or any
other qualifier attached to the source of the fatigue. On average, during a
game of basketball, players travel 5–6 km at physiological intensities above
lactate threshold and 85% of maximal heart rate12. This leaves no question
basketball competitions levy a degree of strain on players that given
insufficient recovery could induce fatigue under this definition.

###  Fatigue’s and athlete performance

There is an extensive body of work examining the physiological demands placed
on athletes in basketball at all levels. A literature review of twenty-five
publications found evidence that players incur fatigue as the game progresses,
albeit to differing degrees depending on position12. “Higher-level players”,
the category NBA players would fall into, sustain greater workloads than
lower-level players at the same positions12. This suggests the NBA is a prime
candidate for examining fatigue accumulation over a season.

Nonetheless, the literature also creates a basis for questioning prevailing
beliefs on fatigue’s influence on game outcomes because several research
groups have failed to find a relationship between the muscular properties of
athletes and player performance in basketball. An analysis of eight semi-
professional basketball players found “trivial-small associations between all
external and internal workload variables and player efficiency”13.
Furthermore, an article from Barca Innovation Hub examining training workload
found no connection between a higher external load and better performance in-
game14.

As opposed to in-game fatigue, analyses on structural fatigue sources suggest
fatigue does influence performance. A research group using t-tests and
analysis of variance (ANOVA) found a significant difference of means in
traveling West to East asserting this influences a player’s circadian
rhythm15. It is well-established that air travel in general influences
athletes and one research group asserts these impacts should materialize in
the NBA16.

One of the biggest performance concerns with athlete fatigue is potential for
increased injuries with higher workloads. Work examining sports in general has
shown that combining an understanding of the training dose and external
influences helps players prevent injury and improve performance17. However,
there are inconsistent results across the literature about the impact of
fatigue on injury occurrence in basketball. A paper examining thirty-three NBA
players in an observational retrospective cohort study found an inverse
relationship between training workload and injuries18. Yet, a work using over
70,000 NBA games and 1,600 injury events showed the likelihood of injury
increased by 2.87% for every ninety-six minutes played and decreased by 15.96%
for every day of rest19.

###  Gap in the literature

There are three main gaps in the literature targeted by the three sections of
this paper. First, we are unaware of any research that measures relative
player workloads and performance at a league wide level. This allows
examination of exertion and cumulative fatigue related to in-game performance.
These ideas are addressed in the section of the paper titled In-Game Workload
Induced Fatigue. Second, we are unaware of any case studies on the
effectiveness of strategies that rest players to avoid fatigue accumulation,
referred to as load management colloquially. This is explored in the section
titled Case study of Kawhi Leonard. Lastly, while there has been work on
structural fatigue in basketball15,16, we propose a more robust procedure for
isolating the impact of structural fatigue that is less vulnerable to latent
variables and offers interpretations on the relevance of observed impacts to
game outcome. This work is detailed in the section titled Structural Fatigue.

##  Data

To access data about player and team performance, as well as potential sources
of fatigue, we used a set of existing R packages maintained by the basketball
analytics community. These packages include nbastatr20, airball21, and NBAr22.
They provided access to box score statistics about player performance and some
more advanced metrics such as game-net rating and Player Impact Estimate at a
game and season level. To measure in-game fatigue we used a player’s total in-
game distance traveled via computer vision generated by Second Spectrum and
courtesy of nbastatr20.

##  In-game workload induced fatigue

In-game workload induced fatigue examines player in-game workload and
performance. The analysis of in-game workload induced fatigue consists of
three additional conceptual divisions. First, an examination of how player
exertion correlates with performance. Second, relating an accumulation of
player workload over a window of time and player exertion in an observation
game. Third, examining how accumulation of player workload over a window of
time influences player performance. These are located in sections 5.3, 5.4 and
5.5 respectively.

###  Theoretical foundations

A player’s relative in-game distance is used as the proxy for player workload.
The formal definition of work from kinematics includes an object’s
displacement and the force acting on the object. Where W is work, F is force,
and s is displacement.

$W=Fs$

The force to move an object of constant mass is also constant. Assuming a
player’s mass is constant over a season, changes in distance for a game
captures the relative change in that player’s work for the game. In other
words, work is agnostic to the duration and cardinality of the player’s motion
and leveraging this allows measurement of player exertion with only tracking
data of player displacement. It is important to acknowledge the literature
establishes that in basketball different types of motions are associated with
slightly different degrees of exertion due to the biomechanics associated with
them23,24. Furthermore, players pushing against each other during games,
moving in the z axis via jumps (of particular relevance to basketball), and
the tracking data being generated by a computer vision model all introduce
additional noise. However, there is no publicly available tracking data with
the granularity to account for these sources of noise nor is there a guarantee
they even cause enough noise to justify the complexity of designing a system
of equivalence. Therefore, studying a player’s total displacement in the xy
plane relative to their own baseline offers the most effective picture of
relative changes in work available.

###  Exertion and performance methodology

An important assumption baked into the hypothesis that a fatigued player’s
performance suffers is that exerting more effort in a game positively affects
performance. This relationship is non-trivial because the scoring dynamics of
basketball are highly tied to athlete skill level offering some isolation from
exertion not present in sports like track.

The relationship between exertion and performance is modeled using simple
linear regression as there was found no indication a more flexible method was
necessary. In the models, a player’s deviation from their season average in a
set of performance measures is the response and the player’s deviation from
their season average in distance traveled divided by minutes played is the
predictor.

$PMD=\beta_{0}+\beta_{1}(DSD)$

Where PMD is performance measure deviation from season mean, DSD is distance
per second of playtime deviation from season mean. As discussed in the
Theoretical foundations, a player’s deviations from their season average
treats player mass as constant and correlates the relative change in work with
a relative change in performance. The distance measured is the total player
movement for the game (regardless if they have the ball). If the performance
statistic measures offensive performance, such as field goal percentage, then
only offensive distance divided by offensive minutes is used. This ideally
tailors the players exertion in that subarea of the game to their performance
in that subarea as well. It is important to note this procedure was not used
when measuring total workload for fatigue discussed in the next section. The
distance divided by playtime is necessary because many cumulative statistics
like rebounds and steals would naturally accumulate as a player is in the game
for more time which is not necessarily indicative of better performance.

Each observation in this analysis is a player in a game, and two subsets of
this data were used to improve the quality of the findings. The first subset
was composed of players who in the game four or more minutes, and the second
was twelve or more minutes (NBA games are 48 minutes long). Excluding players
in the game for fewer than four or twelve minutes reduces the variability in
player performance measures. For example, a player who is substituted into the
game for a couple minutes and misses their one shot but in the next game makes
their single shot would be a 100 percent change in shooting percentage. These
small time periods offers a very noisy measurement of that player’s
performance. Filtering by minutes, as opposed to a different arbitrary cutoff,
introduces less selection bias towards players with certain behavioral
patterns. For example, a cutoff on number of shot attempts would introduce a
bias towards players who take lots of shots. Two subsets are reported because
a more strict limit on the number of minutes in the game lowers variability in
performance measures at the tradeoff of only measuring players who are in the
game for a larger percentage of the time. The results section highlights the
importance of this procedure to the empirical findings.

###  Results exertion and performance

There are two major takeaways from the analysis offered in this section with
full results in Fig. 2 and Fig. 3. The first takeaway is deviations in player
exertion appear to be loosely positively correlated with their performance in
the game. Points scored is arguably the most important performance metric
because the points a player scores is most closely tied to the game’s outcome.
Points scored shows a strong correlation with player exertion for both
subsets, p-values of .01 and 1e-12 for the four and twelve minute subsets
respectively. This offers evidence of a positive relationship between exertion
and scoring. Furthermore, the correlation was positive and statistically
significant for field goal percentage, offensive rebounds, and steals across
both subsets.

Fig. 2

Fig. 3

Despite these relationships, the results present a reasonable degree of
ambiguity about how exertion relates to more granular measures of performance.
There was no correlation in at least one subset between turnovers, three point
shooting, offensive rebounds, and free throws. Furthermore, there was a
significant inverse relationship with the number of offensive rebounds (more
offensive rebounds is considered good in basketball) leaving a literal
interpretation that greater exertion hurts performance in this capacity.
Furthermore, the response is highly sensitive to the subset of the data used.
P-values varied by several orders of magnitude across the subsets. Despite
achieving significant p-values on some measures the dramatic change in results
across subsets has two potential causes. First, it could suggests their is a
different underlying relationship between the these different areas of the
data or second, the large number of observations is leading models to over
assert the prominence of weak relationships, an established problem with large
sample sizes and the p-value 25. Either interpretation raises questions about
how any conclusions regarding exertion and performance will generalise.

The second major takeaway is exertion accounts for a low percentage of
variation in performance across all measures and time constraints. Players in
the game for longer than 12 minutes and field goal percentage had the largest
adjusted-r squared value (.0013) of all measured factors. However, even this
accounted for less than 1 percent of the variation in the response. This has
two potential interpretations. First, player exertion accounts for a very
small, arguably negligible, amount of the variation in player performance.
Alternatively second, measuring player exertion in this fashion introduces too
much noise to capture the magnitude of its relationship with performance.
Without additional experimentation it is challenging to argue for either
hypothesis beyond speculation, however, there is work in the literature that
has shown there is some ambiguity attached to the notion simply exerting more
work means an athlete will perform better in basketball14,18.

###  Post distance overload regress

If players experience a performance detriment associated with fatigue, then
subsets of the season where a player’s workload increases should be followed
by a regression. Extracting this relationship required examination of windows
of the NBA schedule at the player level. An observation included a player in
the game for more than 15 minutes for the 2015 season (n= 19,751). The
playtime cutoff of 15 minutes intends to exclude players who are in the game
for a short period of time as it is unlikely they experience measurable cross-
game fatigue. To examine a potential association between prior workload and
observed game performance a simple linear regression was used. The previous x
days difference from the players season average in distance was treated as the
predictor. The distance in the game difference from season average immediately
following the window of time was treated as the response. Fig. 4 contains a
scatter plot of these results for window sizes of 3, 5, 7, and 10 days with
the regression line.

Fig. 4

This analysis offers evidence contrary to what the fatigued player hypothesis
would suggest. Players that covered more distance in prior games covered more
distance in the observation game as well across all windows. These
relationships had $\leq$2e-16 p-values and adjusted r-squared values showing
about 5 percent of the variance in the response is explained by the window.
One hypothesis for such a strong correlation is this measurement strategy ties
too closely to player’s usage rate which maintains local temporal variations
as their role on the team evolves over the season. This can be adjusted for by
using a playtime adjusted response which is agnostic to a players usage rate.
The trade-off of this analysis is measuring distance playtime adjusted is no
longer a proxy for player work. Scatter plots of the relationships for windows
sizes 5 and 7 are in Fig.6 , and show how the positive association disappears,
supporting this hypothesis.

Fig. 5

In fact, the $\beta_{1}$ coefficient for the linear regressions became
negative with p-values of .003 and 1.1 e-6 for windows of five and seven days
respectively. This indicates an albeit small significant negative association.
A concern with such a small the negative relationship is it could simply be
regression to the mean. Furthermore, making the playtime adjustment violates
the assumptions set in the empirical motivations as it no longer measures
cumulative fatigue. This makes the interpretation of this result in the
context of cumulative fatigue impossible.

###  Post distance overload performance

Lastly, searching for any highly compelling evidence of players experiencing a
workload regression, we ignored all intermediate measures and directly
compared prior workload to observations of game outcomes. Fig. 6 and 7 shows
the results of treating performance measures as the response and windows
remain the predictor.

In line with the prior results the relationship between these performance
measures and fatigue is ambiguous. Field goal percentage achieves a
significant p-value in the 3 day window but no others shown in Fig. 6.
Furthermore, with points there is a positive association between prior
distance, meaning a literal interpretation would suggest players who are more
fatigued score more points. Similar to the first windowed analyses that found
a positive association between prior workload and observation workload, the
most likely hypothesis for the relationship is the usage rate increases as
players are performing better. Therefore, if a player is scoring more points
in prior games they are given more playtime and in turn travel more distance.

Fig. 6

Fig. 7

###  In-game workload big picture

Across all of the approaches to modeling fatigue influences on performance
there is a consistent narrative. Regardless how the data is subset,
normalized, or what is treated as the response there appears to be an
inconsistent and largely ambiguous relationship between fatigue and
performance.

##  Case Study of Kawhi Leonard

Load management is a strategy employed by NBA coaches where they deliberately
rest players, often star players, from playing a game, with the goals of
reducing their risk for injury and preserving them for the playoffs. While
these are ultimately long run goals, the case study below explores how
strategies of load management can pay off in the short run. In particular, it
delves into whether or not Kawhi Leonard, the NBA player most notorious for
practicing load management, plays better in the games immediately following
periods of extra rest.

Fig. 8

Each point in Fig. 8 represents a game in Leonard’s career, with the blue
points showing instances where he was load managed before the game versus the
gray points showing no extra days of rest. The response variable is Player
Impact Estimate; PIE is the NBA’s metric for player performance, factoring in
both defense and offense. The dashed line across the middle represents an
exceptional PIE of 0.15, which was the 95th percentile in the most recent NBA
season.

Using a difference of means test, Kawhi was no more likely to perform better
after being load managed. This is visually apparent in Fig. 8. The majority of
all of the blue points are above the dashed line, suggesting that Leonard
plays exceptionally well after resting for an extra day. However, these blue
points are not convincingly above the gray points, as many of the gray points
are also located above this dashed line. Even in games following no extra
rest, Leonard performs at an elite level for his team. Thus there exists no
substantial evidence to support that Leonard necessarily performs better in
the games following load management versus the games where he did not rest
extra. This suggests that load management is not particularly effective in the
short run, at least for Kawhi Leonard.

##  Structural Fatigue

Structural fatigue encapsulates any factors outside of a team’s control, in
other words, induced by the structure of the schedule. For example, teams have
no influence over how far they have to travel to play a game or how many games
they play in a week. We aggregated travel data from the 2010 through 2019 NBA
seasons, courtesy of airball21. This data consisted of between-game potential
features of interest such as time zone changes, tipoff time, flight distance,
days of rest and schedule density. This section is an examination of how these
schedule-induced-factors influence performance at a team level.

###  Model Construction

The first step to developing the model was to find a strong measure for
capturing game outcome. Compared to alternatives such as score differential or
a binary measure like win loss, game net rating offered the most detailed and
accurate measurement of the result of a game. Game net rating is the
difference between a team’s offensive and defensive rating for a given game,
and both of these ratings are calculated per 100 possessions. Having the
response variable standardized to be per 100 possessions allowed for stronger
comparisons across NBA seasons as playstyles and points per game evolve over
time.

Measuring the effects that schedule induced factors have on game net ratings
required normalizing for the strength of the teams. There was significant
exploration for the strongest proxy of team strength. The simplest model uses
the difference of the teams’ end of season net rating. Attempting to improve
on this, models that temporally weight games based on recency, nonparametric
models including random forest, and dimension reduction techniques over a
large set of team information using elastic net were all tested. These models
all achieved improvements within 2 standard errors in 10 fold cross validation
of team net rating. As this paper was not about modeling team strength, an
active area of sports research, we selected end of season net rating in a
linear regression model as a proxy for team strength. We acknowledge improving
the proxy for team strength is an excellent direction to expand the results of
this work.

After determining a proxy to account for differences in team strength, the
next step was developing a model in two stages. The first stage regressed game
net rating on the net rating difference strength proxy. Then by calculating
the residuals of this first regression and regressing them on all of the
travel metrics, the second stage effectively teased out what portion of the
variance in game net rating the travel metrics could explain that the strength
proxy could not. The results can be seen in the following figures.

###  First model results

$game\>net\>rating=\beta_{0}+\beta_{1}(net\>rating\>di\\!f\\!ference)$

The proxy for team strength difference, net rating difference, yielded an
adjusted R-squared of 0.2262. The data was structured such that the response
was relative to the visiting team, meaning the $\beta_{0}$ term captures that
on average, visiting teams are disadvantaged by about -2.81 in game net
rating. The 1coefficient presents that for each point in net rating that a
given team had more than their opponent, their game net rating for that game
is expected to increase by 0.968.

###  Second model results

$residuals=\beta_{0}+\beta_{i}(schedulefactor_{i})$

Fig. 10

The regression results convey that out of all of the travel metrics, days of
rest difference, traveling three hours back, and the observed game being the
third game in four days were the only variables that had statistically
significant effects on game net rating. The coefficient on rest difference
presents that the visiting team’s game net rating is expected to increase by
0.35 for each day of rest they have more than their opponent. Traveling three
hours back also had a detrimental effect; if a team had to travel from the
east coast all the way west, their game net rating is expected to go down by
about -1.740. Likewise, if the game being played is the team’s third game in
four days, the model estimates their game net rating to decrease by -1.290.
However, it is important to acknowledge these results had adjusted R-squared
indicating they explain less than .1 percent of the variance in the response.
Meaning they offer little information about why a game turned out the way it
did.

It was perhaps more insightful, however, to observe that many of the collected
travel metrics had insignificant effects on game net rating. Contrary to
conventional wisdom, the observed game being the second leg of a back to back
did not have a significant negative effect. Neither did the majority of our
time zone shift metrics, save for three hours back. In addition, all of the
distance metrics- team flight distance, flight distance difference, and
windowed versions of these metrics were also insignificant. This addresses the
aforementioned concern that having to travel and cross more time zones
inherently puts teams at a disadvantage. The results offer that this is not
the case. Moreover, this provides incredibly useful insight to the NBA league
office because it suggests that scheduling teams to fly more miles does not
necessarily put them at a structural disadvantage, but that scheduling them to
have 3 games in 4 days appears to be more consequential to performance. Such
insight can help the league office effectively design schedules that are more
fair.

##  Relevance to Stakeholders

Fatigue analysis provides valuable insight to NBA teams and the league office,
from strategic, schedule fairness, as well as viewer and profit maximization
standpoints.

###  NBA teams

NBA teams can utilize fatigue analysis to understand how to best allocate
their resources and important areas of focus for the team. When teams,
trainers, and athletes understand how fatigue is impacting each game, they can
develop strategies to minimize fatigue. A team can use this knowledge to
reallocate time and effort, such as changing team travel plans, changing how
and when the team gets to the stadium, or adjusting player’s rest to
practice/warm-up ratio in order to decrease fatigue’s impact on game net
rating. Our work in particular suggests that teams’ existing strategies are
highly effective at mitigating the impacts of fatigue. Across all analyses the
impact of player fatigue on game outcomes was ambiguous with the exception of
playing 3 games in 4 days, and the rest difference between teams at a schedule
level. These rare scheduling events may represent the last holdout of
improvements that can be achieved by NBA teams. Therefore, from an optimizing
player performance perspective there is little incentive for raised
investments from teams into player recovery when existing procedures appear to
be eliminating any regressions in performance.

The case study of Kawhi Leonard shows that the impact of load management was
somewhat ambiguous, at least in the short run. An NBA team would want to
perform a similar analysis with other key players. This is because there
exists an opportunity cost every time a player does not suit up for the game.
They will be more rested and have less potential for injury, but there is a
higher chance of losing that game, and other players will have to exert
themselves more. Therefore it would be important for NBA teams to carry out
some form of a cost-benefit analysis to explore whether or not the payoffs of
load management outweigh the deficits.

###  NBA league office

The structural fatigue analysis identified certain features of the schedule
that detrimentally impacted team performance, albeit with a moderately small
overall impact on game outcome. From a league office perspective, they want to
strategically minimize occurrences of these events, but if they are
unavoidable, the league will want to make sure that these events are
distributed fairly so that certain teams are not put at structural
disadvantages.

Additionally, structural fatigue analysis is important when considering
viewership and profit maximization. Using strategies like starting the season
a week earlier as a blueprint, the league office can continue to take measures
that adjust the schedule so that players have adequate time to rest and
recover between games. Schedule analysis drives insight to inform what kind of
changes need to be made in order to do so. This kills two birds with one stone
as it also serves as a means of discouraging teams from load managing their
players. Organizing the schedule in such a way that coaches like Popovich or
Nurse do not even have the need to rest their stars could effectively ensure
that marquee players are suited up on the court. Doing so would very well
increase the entertainment value of NBA games, which in turn boosts
broadcasting viewership and revenue generated.

##  Future Work

###  Exertion and performance

In the section Results on Exertion and Athlete Performance, we are unable to
establish a string association between player exertion and performance. There
is almost certainly a more pronounced relationship between player performance
and exertion. In particular, there is significant room to improve this
analysis with more granular tracking data that better characterizes in-game
events. For example, off the ball movement requires less energy than on ball
movement and more accurately reflects the game. Achieving better mappings
between the muscular properties of athletes and game outcomes will create a
more nuanced picture of how exertion influences performance. At a team level,
this picture could be improved by incorporating practice workload into the
picture to build the most accurate measure of cumulative workload.

###  Long run load management

The case study suggests that Kawhi Leonard’s player impact estimate does not
necessarily improve in games immediately following additional rest compared to
those with no load management. This calls into question the effectiveness of
load management as a strategy in the short run, however it would ultimately be
more important to figure out how load management pays off in the long run.

##  References

  1. 1.

Holmes, B. “Schedule alert! Every game tired teams should lose this month.”
ESPN Enterprises, Inc, 1 Nov. 2017, https://www.espn.com/nba/story/_/id/
21236405/nba-schedule-alert-20-games-tired-teams-lose-november.

  2. 2.

Holmes, B. “Which games will your team lose because of the NBA schedule?.”
ESPN Enterprises, Inc, 30 Oct. 2018, https://www.espn.com/nba/story/_/id/
25117649/nba-schedule-alert-games-your-team-lose-2018-19.

  3. 3.

Magness, S. “Fatigue and the NBA Finals: How Players Raise Their Game When It
Matters Most.” Medium, 3 June 2017, https://medium.com/@stevemagness/how-
player-fatigue-will-impact-the-nba-finals-74d8003d84de.

  4. 4.

“How Much Does Fatigue Matter in NBA Betting?.” OddsShark,
https://www.oddsshark.com/
nba/how-much-does-fatigue-matter-nba-betting.

  5. 5.

McMahan, I. “How sleep and jet-lag influences success in the travel-crazy
NBA.” Guardian News & Media Limited, 26 Oct. 2018,
https://www.theguardian.com/sport/2018/oct/
26/sleep-nba-effects-basketball.

  6. 6.

“Professional teams are waking up to new methods of managing athlete fatigue.”
Fatigue Science, 10 Oct. 2013, https://www.fatiguescience.com/blog/
professional-teams-are-waking-up-to-new-methods-of-managing-athlete-fatigue.

  7. 7.

Schuhmann, J. “Schedule Breakdown: How rest, road trips and other factors will
impact 2021-22.” 2021 NBA Media Ventures, LLC, 22 Aug. 2021,
https://www.nba.com/news/how-rest-road-trips-and-other-factors-played-out-
in-2021-22-schedule.

  8. 8.

Mullin, B. “Injury-Plagued NBA Draws Fewer TV Viewers.” 2021 Dow Jones &
Company, Inc., 27 Dec. 2019, https://www.wsj.com/articles/injury-plagued-nba-
draws-fewer-tv-viewers-11577487612.

  9. 9.

Pickman, B. “Report: NBA Could Issue Teams Fines of $100K For Resting Players
Throughout Season.” 2021 ABG-SI LLC. Sports Illustrated, 7 Dec. 2020,
https://www.si.com/nba/2020/12/07/nba-resting-policy-load-management.

  10. 10.

Edwards T., Spiteri T., Piggott B., Bonhotal J., Haff G.G., and Joyce C.
“Monitoring and managing fatigue in Basketball.” Sports. 6: 19, 2018.

  11. 11.

Seamon, B.A., and Harris-Love, M.O. “Clinical Assessment of Fatigability in
Multiple Sclerosis: A Shift from Perception to Performance.” Front. Neurol. 7:
194, 2016.

  12. 12.

Stojanović, E., Stojiljković, N., Scanlan, A.T. et al. “The Activity Demands
and Physiological Responses Encountered During Basketball Match-Play: A
Systematic Review.” Sports Med 48: 111–135, 2018.

  13. 13.

Fox, J., Stanton, R., O’Grady, C., Teramoto, M., Sargent, C., and Scanlan, A.
“Are acute player workloads associated with in-game performance in
basketball?.” Biology of Sport, 2020.

  14. 14.

García, A.C. “Training Load and Performance in Basketball: The More The
Better?.” FC Barcelona Innovation Hub, 22 Oct. 2020,
https://barcainnovationhub.com/training-load-and-performance-in-basketball-
the-more-the-better/.

  15. 15.

Holmes, B. “How fatigue shaped the season, and what it means for the
playoffs.” ESPN Enterprises, Inc, 10 Apr. 2018,
https://www.espn.com/nba/story/_/ id/23094298/how-fatigue-shaped-nba -season-
means-playoffs.

  16. 16.

Huyghe T., Scanlan A., Dalbo V., Calleja-González J. “The negative influence
of air travel on health and performance in the National Basketball
Association: A narrative review.” Sports. 6: 89, 2018.

  17. 17.

Drew, M.K., Cook, J., and Finch, C.F. “Sports-related workload and injury
risk:simply knowing the risks will not prevent injuries: Narrative review.”
British Journal of Sports Medicine 50: 1306-1308, 2016.

  18. 18.

Caparrós T., Casals M., Solana Á., and Pena J. “Low External Workloads Are
Related to Higher Injury Risk in Professional Male Basketball Games.” J.
Sports Sci. Med. 17: 289–297, 2018.

  19. 19.

Lewis, M. “It’s a hard-knock life: game load, fatigue, and injury risk in the
National Basketball Association.” J. Athletic Train. 53, 503–509, 2018.

  20. 20.

Bresler, A. (2021). Nbastatr. R package version 0.1.1504.
http://asbcllc.com/nbastatR/.

  21. 21.

Fernandez, J. (2020). airball: Schedule & Travel Related Metrics in
Basketball. R package version 0.4.0. https://github.com/josedv82/airball.

  22. 22.

Chodowski, P. (2021). NBAr. R package version 4.0.4. https://github.com/
PatrickChodowski/NBAr.

  23. 23.

Abdelkrim, B.N., Castagna, C.; El Fazaa, S., and El Ati, J. “The Effect of
Players’ Standard and Tactical Strategy on Game Demands in Men’s Basketball.”
Journal of Strength and Conditioning Research 24 - Issue 10: 2652-2662, 2010.

  24. 24.

Abdelkrim, B.N., El Fazaa, S., and El Ati J. “Time–motion analysis and
physiological data of elite under-19-year-old basketball players during
competition.” British Journal of Sports Medicine 41:69-75, 2007.

  25. 25.

Lin, Mingfeng & Lucas, Henry & Shmueli, Galit. (2013). Too Big to Fail: Large
Samples and the p-Value Problem. Information Systems Research. 24. 906-917.
10.1287/isre.2013.0480.


# Chaos in Coupled Heteroclinic Cycles and its Piecewise-Constant
Representation

Arkady Pikovsky<EMAIL_ADDRESS>Alexander Nepomnyashchy
<EMAIL_ADDRESS>Institute for Physics and Astronomy, University of
Potsdam, Karl-Liebknecht Str. 24/25, 14476 Potsdam, Germany  National Research
University Higher School of Economics, Nizhny Novgorod, Russia  Department of
Mathematics, Technion - Israel Institute of Technology, Haifa, Israel

###### Abstract

We consider two stable heteroclinic cycles rotating in opposite directions,
coupled via diffusive terms. A complete synchronization in this system is
impossible, and numerical exploration shows that chaos is abundant at low
levels of coupling. With increase of coupling strength, several symmetry-
changing transitions are observed, and finally a stable periodic orbit appears
via an inverse period-doubling cascade. To reveal the behavior at extremely
small couplings, a piecewise-constant model for the dynamics is suggested.
Within this model we construct a Poincaré map for a chaotic state numerically,
it appears to be an expanding non-invertable circle map thus confirming
abundance of chaos in the small coupling limit. We also show that within the
piecewise-constant description, there is a set of periodic solutions with
different phase shifts between subsystems, due to dead zones in the coupling.

##  1 Introduction

A heteroclinic cycle is an interesting object in nonlinear dynamics, attracted
large of attention recently. It was probably first described by May and
Leonard [1] in the context of a generalization of the Lotka-Volterra model of
species competition. One of the early examples is also the dynamics of three
unstable modes in a convection problem in a rotating fluid layer heated from
below [2]. From the mathematical viewpoint, a stable heteroclinic cycle is a
robust object, demonstrating oscillations the period of which grows in time,
and tends to infinity [3, 4]. Among recent applications we mention the concept
of winnerless competition in neurosciences, based on a representation of
certain dynamical patterns as vicinities of a heteroclinic cycle [5, 6];
studies of the dynamics of spin-torque nano-spintronic oscillators [7];
synchronization dynamics of coupled oscillators [8]. In the simplest setup, a
heteroclinic cycle includes three interacting “species”, but recently
generalizations to larger networks have been also investigated [9, 10, 11,
12].

In this paper, we consider a simple system of two coupled subsystems, each of
which possesses a stable heteroclinic cycle. For brevity, we will speak on
“coupled heteroclinic cycles”. A variant with a symmetric (or nearly
symmetric) coupling has been studied in [13]. There, several synchronous and
quasiperiodic regimes have been reported. Here we consider a strongly
antisymmetric coupling of two cycles. We assume, that we have two cycles
oscillating in opposite directions. Namely, while in one cycle we have a
sequence of states $1\to 2\to 3\to 1\to 2\ldots$, in another cycle we have a
sequence of states $1\to 3\to 2\to 1\to 3\ldots$. The coupling is an
attraction of the corresponding states in two cycles. Because of opposite
rotations, a completely synchronous state, where the states of the coupled
variables in two subsystems coincide (and thus the coupling diffusion terms
vanish), is impossible (even if all parameters of the cycles are the same, as
we assume below). We will however show, that partially synchronized periodic
regimes, where non-interacting variables coincide, can appear and are stable
for large coupling. Our main finding is that chaos dominates the dynamics at
small coupling. We first show this numerically for the full system of
equations. Then, we develop an approximate approach where the dynamics is
approximated by piecewise-constant functions. This approach captures the
properties of coupled heteroclinic cycles in the limit of small coupling.

The paper is organized as follows. We discuss the basic heteroclinic cycle in
Section 2. Coupled heteroclinic cycles are introduced in Section 3. The
coupled system possesses several symmetries, which we discuss in Section 4.
Numerical exploration of coupled cycles is performed in Section 5. In Section
6 we introduce the piecewise-constant model of the cycle dynamics, which is
then extended to coupled units in Section 7. Chaotic and periodic regimes
within the piecewise-constant model are explored in Section 8, and compared
with the regimes in the full system in Section 9. We conclude with discussion
of the results in Section 10.

##  2 Heteroclinic cycle

We discuss here briefly the canonical model of a heteroclinic cycle,
introduced by May and Leonard [1]:

$\displaystyle\dot{u}$ $\displaystyle=u(1-u-\alpha v-\beta w)\;,$ (1)
$\displaystyle\dot{v}$ $\displaystyle=v(1-v-\alpha w-\beta u)\;,$
$\displaystyle\dot{w}$ $\displaystyle=w(1-w-\alpha u-\beta v)\;.$

We assume below that $\alpha>1$ and $0<\beta<1$. The crucial parameter here is
the sum $\alpha+\beta$. For $\alpha+\beta>2$, a stable heteroclinic cycle
exists, in which a trajectory consequently approaches steady states
$(1,0,0)\to(0,1,0)\to(0,0,1)\to(1,0,0)\to\ldots$ as time grows. The time
intervals that a trajectory spends in a vicinity of these states grow
exponentially.

For an illustration and for further analysis below, it is convenient to
introduce new variables allowing for a better resolution of vicinities of the
steady states, stable and unstable manifolds of which constitute the
heteroclinic cycle. These variables are also suitable for numerical
integration and for a piecewise-constant approximation (see Section 6 below).
These new variables read

$x=\ln\frac{u}{1-u}\;,\quad y=\ln\frac{v}{1-v}\;,\quad z=\ln\frac{w}{1-w}\;.$
(2)

The equations (1) in these variables take the form

$\displaystyle\dot{x}$
$\displaystyle=1-\alpha\frac{1+e^{x}}{1+e^{-y}}-\beta\frac{1+e^{x}}{1+e^{-z}}\;,$
(3) $\displaystyle\dot{y}$
$\displaystyle=1-\alpha\frac{1+e^{y}}{1+e^{-z}}-\beta\frac{1+e^{y}}{1+e^{-x}}\;,$
$\displaystyle\dot{z}$
$\displaystyle=1-\alpha\frac{1+e^{z}}{1+e^{-x}}-\beta\frac{1+e^{z}}{1+e^{-y}}\;.$

We illustrate a trajectory approaching the heteroclinic cycle in Fig. 1.

Figure 1: A trajectory in systems (1),(3). Colors: Red: $x(t)$, green: $y(t)$,
blue: $z(t)$. Large positive values of $x$ correspond to $u\approx 1$, and
large negative values of $x$ correspond to $u\approx 0$; similarly for other
variables.

##  3 Coupled heteroclinic cycles

Here we introduce our basic model of two interacting heteroclinic cycles. In
contradistinction to the previously explored case of a nearly symmetric
coupling [13], we consider an asymmetric coupling: we assume that two
heteroclinic cycles of type (1) rotate in opposite directions. Namely, while
in cycle 1 with variables $u_{1},v_{1},w_{1}$ the excitations are according to
the rule $u_{1}\to v_{1}\to w_{1}\to u_{1}\ldots$, in cycle 2 the excitations
go as $u_{2}\to w_{2}\to v_{2}\to u_{2}\to\ldots$. The coupled variables are
pairs $(u_{1},u_{2})$, $(v_{1},v_{2})$, and $(w_{1},w_{2})$. The equations of
our model read

$\displaystyle\dot{u}_{1}$ $\displaystyle=u_{1}(1-u_{1}-\alpha v_{1}-\beta
w_{1})+D(u_{2}-u_{1})\;,$ (4) $\displaystyle\dot{v}_{1}$
$\displaystyle=v_{1}(1-v_{1}-\alpha w_{1}-\beta u_{1})+D(v_{2}-v_{1})\;,$
$\displaystyle\dot{w}_{1}$ $\displaystyle=w_{1}(1-w_{1}-\alpha u_{1}-\beta
v_{1})+D(w_{2}-w_{1})\;,$ $\displaystyle\dot{u}_{2}$
$\displaystyle=u_{2}(1-u_{2}-\beta v_{2}-\alpha w_{2})+D(u_{1}-u_{2})\;,$
$\displaystyle\dot{v}_{2}$ $\displaystyle=v_{1}(1-v_{2}-\beta w_{2}-\alpha
u_{2})+D(v_{1}-v_{2})\;,$ $\displaystyle\dot{w}_{2}$
$\displaystyle=w_{1}(1-w_{2}-\beta u_{2}-\alpha v_{2})+D(w_{1}-w_{2})\;.$

One can see that because of exchange $\alpha\leftrightarrow\beta$, the cycles
in subsystems 1 and 2 rotate in opposite directions. Parameter $D$ describes
coupling between subsystems. We expect that already a small coupling will
significantly influence the dynamics, so $D$ will be mostly small. Therefore,
it is convenient to represent this parameter as $D=\exp[-Q]$, so that large
values of $Q$ correspond to a weak coupling. Then, adopting the variables (2),
we can rewrite system (4) as

$\displaystyle\dot{x}_{1}$
$\displaystyle=1-\alpha\frac{1+e^{x_{1}}}{1+e^{-{y_{1}}}}-\beta\frac{1+e^{x_{1}}}{1+e^{-z_{1}}}+e^{-Q}\left[\frac{(1+e^{-x_{1}})(1+e^{x_{1}})}{1+e^{-x_{2}}}-(1+e^{x_{1}})\right]\;,$
(5) $\displaystyle\dot{y}_{1}$
$\displaystyle=1-\alpha\frac{1+e^{y_{1}}}{1+e^{-z_{1}}}-\beta\frac{1+e^{y_{1}}}{1+e^{-x_{1}}}+e^{-Q}\left[\frac{(1+e^{-y_{1}})(1+e^{y_{1}})}{1+e^{-y_{2}}}-(1+e^{y_{1}})\right]\;,$
$\displaystyle\dot{z}_{1}$
$\displaystyle=1-\alpha\frac{1+e^{z_{1}}}{1+e^{-x_{1}}}-\beta\frac{1+e^{z_{1}}}{1+e^{-{y_{1}}}}+e^{-Q}\left[\frac{(1+e^{-z_{1}})(1+e^{z_{1}})}{1+e^{-z_{2}}}-(1+e^{z_{1}})\right]\;,$
$\displaystyle\dot{x}_{2}$
$\displaystyle=1-\beta\frac{1+e^{x_{2}}}{1+e^{-y_{2}}}-\alpha\frac{1+e^{x_{2}}}{1+e^{-z_{2}}}+e^{-Q}\left[\frac{(1+e^{-x_{2}})(1+e^{x_{2}})}{1+e^{-x_{1}}}-(1+e^{x_{2}})\right]\;,$
$\displaystyle\dot{y}_{2}$
$\displaystyle=1-\beta\frac{1+e^{y_{2}}}{1+e^{-z_{2}}}-\alpha\frac{1+e^{y_{2}}}{1+e^{-x_{2}}}+e^{-Q}\left[\frac{(1+e^{-y_{2}})(1+e^{y_{2}})}{1+e^{-y_{1}}}-(1+e^{y_{2}})\right]\;,$
$\displaystyle\dot{z}_{2}$
$\displaystyle=1-\beta\frac{1+e^{z_{2}}}{1+e^{-x_{2}}}-\alpha\frac{1+e^{z_{2}}}{1+e^{-y_{2}}}+e^{-Q}\left[\frac{(1+e^{-z_{2}})(1+e^{z_{2}})}{1+e^{-z_{1}}}-(1+e^{z_{2}})\right]\;.$

While Eqs. (4) are appropriate for an analytical consideration of regimes
where the co-existence state is weakly unstable, Eqs. (5) are suitable for
numerical explorations of a deeply heteroclinic dynamics.

##  4 Invariant synchronous manifolds

###  4.1 Symmetry properties of systems (4), (5)

System (4) is invariant with respect to renamings of variables inside each
subsystem:

$\begin{gathered}R_{1}:(u_{1},v_{1},w_{1},u_{2},v_{2},w_{2})\to(v_{1},w_{1},u_{1},v_{2},w_{2},u_{2})\\\
\text{and
}R_{2}:(u_{1},v_{1},w_{1},u_{2},v_{2},w_{2})\to(w_{1},u_{1},v_{1},w_{2},u_{2},v_{2}).\end{gathered}$
(6)

Obviously, $R_{2}=R_{1}^{-1}=R_{1}^{2}$, $R_{1}=R_{2}^{-1}=R_{2}^{2}$. Also,
system (4) is symmetric with respect to three exchange transformations,

$\begin{gathered}T_{1}:(u_{1},v_{1},w_{1},u_{2},v_{2},w_{2})\to(u_{2},w_{2},v_{2},u_{1},w_{1},v_{1})\;,\\\
T_{2}:(u_{1},v_{1},w_{1},u_{2},v_{2},w_{2})\to(w_{2},v_{2},u_{2},w_{1},v_{1},u_{1})\;,\\\
\text{and
}T_{3}:(u_{1},v_{1},w_{1},u_{2},v_{2},w_{2})\to(v_{2},u_{2},w_{2},v_{1},u_{1},w_{1})\;.\end{gathered}$
(7)

One can see that $T_{1}^{2}=T_{2}^{2}=T_{3}^{2}=I$, $T_{2}=R_{2}T_{1}$,
$T_{3}=R_{1}T_{1}$.

System (5) has the same set of symmetries with replacement of
$(u_{i},v_{i},w_{i})$ with $(x_{i},y_{i},z_{i})$, $i=1,2$.

The subspaces invariant with respect to transformations $T_{1}$, $T_{2}$ and
$T_{3}$ are invariant manifolds for system (4). Let us consider, for instance,
i.e., the transformation

$T_{1}:(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2})\to(x_{2},z_{2},y_{2},x_{1},z_{1},y_{1}).$

System (5) can be written as

$\dot{x}_{1}=f(x_{1},y_{1},z_{1},x_{2}),\;\dot{y}_{1}=f(y_{1},z_{1},x_{1},y_{2}),\;\dot{z}_{1}=f(z_{1},x_{1},y_{1},z_{2}),$
$\dot{x}_{2}=f(x_{2},z_{2},y_{2},x_{1}),\;\dot{y}_{2}=f(y_{2},x_{2},z_{2},y_{1}),\;\dot{z}_{2}=f(z_{2},y_{2},x_{2},z_{1}),$

where

$f(x,y,z,u)=1-\alpha\frac{1+e^{x}}{1+e^{-y}}-\beta\frac{1+e^{x}}{1+e^{-z}}+e^{-Q}(1+e^{x})\frac{e^{-x}-e^{-u}}{1+e^{-u}}\;.$

Let us define

$x=\frac{x_{1}+x_{2}}{2},\;y=\frac{y_{1}+z_{2}}{2},\;z=\frac{z_{1}+y_{2}}{2},\;X=\frac{x_{2}-x_{1}}{2},\;Y=\frac{z_{2}-y_{1}}{2},\;Z=\frac{y_{2}-z_{1}}{2}.$

Then we obtain:

$\dot{X}=\frac{1}{2}[f(x+X,y+Y,z+Z,x-X)-f(x-X,y-Y,z-Z,x+X)]\;,$
$\dot{Y}=\frac{1}{2}[f(y+Y,z+Z,x+X,z-Z)-f(y-Y,z-Z,x-X,z+Z)]\;,$
$\dot{Z}=\frac{1}{2}[f(z+Z,x+X,y+Y,y-Y)-f(z-Z,x-X,y-Y,y+Y)]\;,$
$\dot{x}=\frac{1}{2}[f(x+X,y+Y,z+Z,x-X)+f(x-X,y-Y,z-Z,x+X)]\;,$
$\dot{y}=\frac{1}{2}[f(y+Y,z+Z,x+X,z-Z)+f(y-Y,z-Z,x-X,z+Z)]\;,$
$\dot{z}=\frac{1}{2}[f(z+Z,x+X,y+Y,y-Y)+f(z-Z,x-X,y-Y,y+Y)]\;.$

One can see that the system has a 3-dimensional invariant manifold $X=Y=Z=0$.
The dynamics on that manifold is considered in the next subsection 4.2. The
invariance to transformation $T_{1}$ corresponds to the reversion symmetry
with respect to that manifold:

$(x,y,z,X,Y,Z)\to(x,y,z,-X,-Y,-Z).$

Therefore, the system can have a symmetric attractor that includes
simultaneously both points $(x,y,z,X,Y,Z)$ and $(x,y,z,-X,-Y,-Z)$. Otherwise,
it has two asymmetric attractors connected by $T_{1}$.

###  4.2 Invariant manifolds

System (4) has a 3-dimensional invariant manifold

$M1:\;u_{1}=u_{2}\equiv u,\;v_{1}=w_{2}\equiv v,\;w_{1}=v_{2}\equiv w,$ (8)

on which the dynamics reduces to a system of 3 equations:

$\displaystyle\frac{du}{dt}$ $\displaystyle=u(1-u-\alpha v-\beta w),$ (9)
$\displaystyle\frac{dv}{dt}$ $\displaystyle=v(1-v-\alpha w-\beta u)+D(w-v),$
$\displaystyle\frac{dw}{dt}$ $\displaystyle=w(1-w-\alpha u-\beta v)+D(v-w).$

Such a regime can be described as a partially synchronous one. Here, as one
can see from expressions (8), one pair of interacting variables are identical,
while two other pairs are “cross-identical”. Due to this, the coupling terms
are essential for the dynamics.

Due to the renaming symmetry, system (4) has two other invariant manifolds,

$M2:\;v_{1}=v_{2}\equiv u,\;w_{1}=u_{2}\equiv v,\;u_{1}=w_{2}\equiv w,\qquad
M3:\;w_{1}=w_{2}\equiv u,\;u_{1}=v_{2}\equiv v,\;v_{1}=u_{2}\equiv w,$ (10)

where the dynamics is described by exactly the same system (9).

###  4.3 Stationary points and their stability

For $D=0$, system (9) coincides with (1) and has five steady states:

$(0,0,0),\quad(1,0,0),\quad(a,a,a),\quad(0,1,0),\quad(0,0,1)$

where $a=(\alpha+\beta+1)^{-1}$. The first three of these points do not depend
on parameter $D$. For small values of $D$, we can approximately find the last
two steady states up to the first order in $D$:

$(0,1+D\frac{\alpha+\beta-1}{1-\beta},D\frac{-1}{1-\beta}),\quad(0,D\frac{1}{\alpha-1},1+D\frac{-\alpha+1-\beta}{\alpha-1})$

Because $\beta<1$, the first of these steady states lies in a “forbidden
domain” $w<0$, and only the second of these steady states exists. Thus, a
heteroclinic cycle disappears for any small $D$.

For $D=0$, the “coexistence point” $(a,a,a)$ loses its stability at
$\alpha+\beta=2$, and the heteroclinic cycle appears. In system (4) the
relevant eigenvalues of the coexistence point are

$\sigma_{2,3}=\frac{\alpha+\beta-2}{2(\alpha+\beta+1)}-D\pm
i\left[\frac{3(\alpha-\beta)^{2}}{4(\alpha+\beta+1)^{2}}-D^{2}\right]^{1/2}.$
(11)

Thus, the “coexistence point” is oscillatory unstable if
$D<\frac{\alpha+\beta-2}{2(\alpha+\beta+1)}$. For small $D$, the shift of the
stability threshold is $\alpha+\beta-2=6D+O(D^{2})$. Furthermore, as is shown
in A, the transition is now a standard Hopf bifurcation at which a stable
limit cycle appears. In this paper, however, we will not explore the case
where parameters $\alpha$ and $\beta$ are close to this bifurcation point,
rather we focus on a situation where the trajectory in coupled systems is
close to former well-developed heteroclinic cycles (i.e. with $\alpha+\beta-2$
not small).

Figure 2: Panel (a): period of the cycle (magenta dashed line, right axis) and
maximal and minimal values of the normalized by $Q$ variables $x$ (red), $y$
(green) and $z$ (blue, left axis). Panel (b): absolute values of the
multipliers of the cycle. Stability is lost when one real multiplier crosses
the border $-1$.

Below in numerical simulations in Sections 4,5 we use

$\begin{gathered}\alpha=1+\Omega\approx 2.324,\quad\beta=1-\Omega^{-1}\approx
0.2451,\\\ \text{where}\quad\Omega=1.32471595724475\ldots\quad\text{is the
spiral mean }\Omega^{3}=\Omega+1\;.\end{gathered}$ (12)

Furthermore, for better comparison with theory of Section 7, we normalize the
depicted variables $x,y,z$ by $Q$.

Figure 3: Bifurcation diagram obtained at the Poincare section $x_{1}=0$,
$\dot{x}_{1}>0$; depicted are values of the variable $x_{2}$ at these cross-
sections.

###  4.4 Periodic dynamics and its transversal stability

Figure 2 shows the bifurcation diagram of the periodic orbit on the symmetric
manifold M1, and its stability (including transversal one). Figure 3 shows the
bifurcation diagram of the full 6-dimensional system. It has been constructed
by following the cycle on manifold M1 starting from small values of $Q$. This
cycle is stable for $Q\lesssim 5.43$ (as can be seen from the plot of
multipliers Fig. 2), at this value it experiences period-doubling. A next
transition is at $Q\approx 7.36$, where through a pitchfork bifurcation, two
mutually symmetric cycles arise. These cycles are shown in Fig. 3 with two
different colors (red and blue). Further increase of values of $Q$ leads to a
transition to chaos. We will discuss this part of the diagram in the next
section 5, starting from large values of $Q$ and following decrease of these
values.

##  5 Numerical exploration of chaos in coupled heteroclinic cycles

###  5.1 Bifurcation diagram

For fixed values $\alpha,\beta$ (12), the only bifurcation parameter is the
strength of coupling $Q$. The diagram of the states in dependence on $Q$ is
presented in Fig. 3. Here we show values of the variable $x_{2}$ at the
Poincaré surface of section $x_{1}=0,\dot{x}_{1}>0$. One can see that in a
large range of values of $Q$ the dynamics is chaotic, with some periodic
windows present.

###  5.2 Symmetric chaotic attractor for weak coupling (large $Q$)

Figure 4: Panel (a): Points at the Poincaré map projected on plane
$(x_{2},y_{2})$. Although the figure looks like a curve, in fact it is a dense
set of points. Colors of these points indicate the closest invariant
synchronous manifolds (blue: M1; green: M2; red: M3). Panel (b): a one-
dimensional map constructed as described in text.

We start with chaotic states at weak coupling (large values of $Q$, bottom
panel at Fig. 3). We take $Q=35$ and explore the Poincaré section $x_{1}=0$,
$\dot{x}_{1}>0$ like in Fig. 3. At this section, we plot the values of
variables $(x_{2},y_{2})$ in Fig. 4(a). This image looks like a closed curve,
but in fact it is a set of points. This indicates that, with a good accuracy,
the system’s dynamics is described by a one-dimensional map of a circle onto
itself. To construct this map, we rescale the variables and attribute to each
point in Fig. 4(a) an angle $\phi_{n}=\arctan(\eta_{n}/\xi_{n})$, where
$\xi=(x_{2}+15)/56$ and $\eta=(y_{2}+20)/117$. Figure 4(b) shows the
transformation $\phi_{n}\to\phi_{n+1}$. One can see that it looks like a non-
invertible circle map. This map has many expanding regions, producing chaos.
However, the map is continuous with small regions where its derivative vanish;
such regions give rise to periodic windows.

It is instructive to explore the role of symmetric manifolds M1,M2,M3 (cf.
Eqs. (8),(10)). We define the closeness to these manifolds via the observables

$\begin{gathered}d_{M1}=(x_{1}-x_{2})^{2}+(y_{1}-z_{2})^{2}+(z_{1}-y_{2})^{2}\;,\quad
d_{M2}=(x_{1}-z_{2})^{2}+(y_{1}-y_{2})^{2}+(z_{1}-x_{2})^{2}\;,\\\
d_{M3}=(x_{1}-y_{2})^{2}+(y_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}\;,\end{gathered}$
(13)

which vanish exactly on the manifolds.

Figure 5: Characterization of closeness of a trajectory to the submanifolds
M1,M2,M3, at $Q=35$. The lines show quantities $d_{M1}$ (blue); $d_{M2}$
(green); $d_{M3}$ (red).

The time evolution of these observables at $Q=35$ is shown in Fig. 5. One can
see that most of the time one of these variables is much smaller than another
ones (notice the logarithmic scale of these distances). This allows for a
representation of the dynamics as an irregular sequence of nearly synchronous
patches. In Fig. 4(a) we use the same colors as in Fig. 5, to indicate
positions on this Poincaré map, corresponding to the vicinities of different
synchronous manifolds.

Figure 6: Panel (a) Three Poincare maps with the same section for three
symmetric attractors at $Q=14$. Panel (b) A first return map $x_{2}(n)\to
x_{2}(n+1)$ for the green attractor in panel (a). This piece is chosen because
its projection on the variable $x_{2}$ is well-defined. Figure 7: Distances to
the manifolds M1-M3 for the green attractor of Fig. 6 at $Q=14$. For other
attractors, the picture is the same, with the corresponding rotation of
colors.

###  5.3 Symmetry breaking transitions

Here, we describe several symmetry-breaking transitions that occur as the
coupling becomes stronger (parameter $Q$ decreases).

The first transition, at which the attractor of Fig. 4 breaks into three
symmetric attractors, happens at $Q\approx 14.5$. This transition is clearly
seen as a jump in the bifurcation diagram (Fig. 3). The Poincaré map for
$Q=14$ is shown in Fig. 6. One can see three pieces, each correspond to one of
the three attractors. Because of the symmetry, we present the effective one-
dimensional map just for one of these attractors in panel (b). One can see
that it looks like a distorted $\sin$-map.

Figure 8: Panel (a) Three Poincare maps with the same section for three
symmetric attractors at $Q=9.5$. Panel (b) A first return map $x_{2}(n)\to
x_{2}(n+1)$ for the green attractor in panel (a). This piece is chosen because
its projection on the variable $x_{2}$ is well-defined. Figure 9: Panel (a)
Six Poincare maps with the same section for three symmetric attractors at
$Q=8.5$. We use the same colors as in Fig. 8, but the two attractors are shown
with thin and bold markers and using a dark and a light version of the color.
Panel (b) A first return map $x_{2}(n)\to x_{2}(n+1)$ for the green attractors
in panel (a). Here again two attractors are shown with different size of the
points and with dark and light versions of the color.

We show the time evolution of the distances to the submanifolds M1-M3 for the
green attractor of Fig. 6 in Fig 7. One can see that now the dynamics is
persistently close to the symmetric manifold M2.

At the next transition (at $Q\approx 10$), each of the symmetric attractors
splits into two period-two bands. The attractors and the corresponding first-
return map are shown in Fig. 8. At the next transition (at $Q\approx 8.55$)
each of three attractors splits into two, as demonstrated in Fig. 9. With
further decrease of value of $Q$, these attractors undergo an inverse cascade
of period-doublings. At $Q\approx 8.08$ stable period-2 asymmetric cycle at
the end of this cascade appears. At $Q\approx 7.36$, through an inverse
pitchfork bifurcation, two asymmetric stable period-2 cycles merge into a
stable symmetric period-2 cycle, as described in Section 4.4 above.

##  6 Piecewise-constant model of heteroclinic cycle

In this section, we construct an approximate piecewise-constant (PC) model of
the heteroclinic cycle. It is based on its representation in variables (2)
within system (3). The functions appearing in (3) resemble an expression for a
Fermi distribution, which has a simple form in the zero-temperature limit
$T\to 0$

$\frac{1}{1+e^{-x}}\quad\Leftrightarrow\quad\frac{1}{1+e^{\epsilon/T}}=\begin{cases}0&\epsilon>0\;,\\\
1&\epsilon<0\;.\end{cases}$ (14)

In our case, there is no small parameter $T$, but one can consider the range
of variations of the variable $x$ (which goes to large positive and large
negative values) as a large parameter. Then, a smooth transition in a region
$x\approx 0$ can be approximated with a step function, similarly to the zero-
temperature limit in (14).

Figure 10: Piecewise-constant representation of function $F(x,y)$,
contributing to the dynamics of variable $x$. The white area is a no-go domain
where formally $F=\infty$. Arrows show the direction of the evolution of $x$
according to this term.

Let us consider, as an example, the second term in the equation for $x$ in
(3). We can approximate it according to the values of $(x,y)$ as follows:

$\frac{1+e^{x}}{1+e^{-y}}\approx F(x,y)=\begin{cases}\infty&\text{for
}x>0,y>0\text{ and }y<0,x>-y\;,\\\ 0&\text{for }y<0,x<-y\;,\\\ 1&\text{for
}y>0,x<0\;.\end{cases}$ (15)

The domain where the value of this term is $\infty$ is in fact a no-go area,
it is never entered during the dynamics. In the rest available area, this term
takes values $1$ and $0$. We illustrate this in Fig. 10. The available domain
$(y<0,x<-y)\cup(y>0,x<0)$ has two borders: $(x=-y;y<0)$ and $(x=0;y>0)$. The
latter border is not reachable, because there $\dot{x}<0$. However, at the
former border the time derivative of $x$ according to the term (15) vanishes.
Thus, a piece of trajectory along this border $(x=-y;y<0)$ is possible, if the
variable $y$ grows.

Regime |  Conditions |  Equations |  Next regime
---|---|---|---
A |  $x>0$
$|y|>|x|$ |  $\dot{x}=1$
$\dot{y}=1-\beta$
$\dot{z}=1-\alpha$ |  $\to$ Ab
Ab (boundary of A) |  $x>0$
$|y|=|x|$ |  $x=-y$
$\dot{y}=1-\beta$
$\dot{z}=1-\alpha$ |  $\to$ B
B |  $y>0$
$|z|>|y|$ |  $\dot{x}=1-\alpha$
$\dot{y}=1$
$\dot{z}=1-\beta$ |  $\to$ Bb
Bb (boundary of B) |  $y>0$
$|z|=|y|$ |  $\dot{x}=1-\alpha$
$y=-z$
$\dot{z}=1-\beta$ |  $\to$ C
C |  $z>0$
$|x|>|z|$ |  $\dot{x}=1-\beta$
$\dot{y}=1-\alpha$
$\dot{z}=1$ |  $\to$ Cb
Cb (boundary of C) |  $z>0$
$|x|=|z|$ |  $\dot{x}=1-\beta$
$\dot{y}=1-\alpha$
$z=-x$ |  $\to$ A
Table 1: Piecewise-constant dynamics, approximating the heteroclinic cycle in
Eqs. (3) according to (15).

The approximation (15) can be applied to all the terms on the r.h.s. of (3).
It can be summarized in a sequence of possible constant-velocity pieces listed
in Table 1 (there, for brevity, we list only pieces along a trajectory close
to the heteroclinic cycle for $\alpha>1$, $\beta<1$). Matching these pieces
$A\to Ab\to B\to Bb\to C\to Cb\to A\to\ldots$ yields a full trajectory
approaching the heteroclinic cycle (see Fig. 11).

Figure 11: Heteroclinic cycle in Eqs. (3) (panel (a)) and in the PC
representation Table 1 (panel (b)). Depicted are the time derivatives
$\dot{x}$ (red), $\dot{y}$ (green), and $\dot{z}$ (blue) as functions of time
along a trajectory. One can see a good correspondence, except for an initial
transient in panel (a), during which a trajectory approaches the heteroclinic
cycle. Dashed lines in panel (b) show the switchings between the piecewise-
constant patches; these patches are shown with the same coding as in Table 1.

##  7 Piecewise-constant model of coupled heteroclinic cycles

###  7.1 Effect of coupling in the PC approximation

We now extend the analysis of Section 6 above to the case of coupled
heteroclinic cycles (5). Let us consider the equation for $x_{1}$, where the
additional term is

$C(x_{1},x_{2})=e^{-Q}\frac{(1+e^{-x_{1}})(1+e^{x_{1}})}{1+e^{-x_{2}}}-(1+e^{x_{1}})\;.$

Similarly to (15) we obtain in the PC approximation

$C(x_{1},x_{2})=\begin{cases}-\infty&\text{for }(x_{2}>0,x_{1}>x_{2}+Q)\text{
and }(x_{2}<0,x_{1}>Q)\;,\\\ 0&\text{for }(x_{2}>0,Q<x_{1}<x_{2}+Q)\text{ and
}(x_{2}<0,-Q+x_{2}<x_{1}<Q)\;,\\\ \infty&\text{for }(x_{2}>0,x_{1}<-Q)\text{
and }(x_{2}<0,x_{1}<-Q+x_{2})\;.\end{cases}$ (16)

(a)  (b)

Figure 12: Panel (a): Accessible domain for $x_{1}$ on plane $x_{1},x_{2}$
according to Eqs. (16). Panel (b): Full accessible for $x_{1}$ and $x_{2}$.

The available domain where $C=0$ is depicted in Fig. 12(a). If we combine this
domain with the symmetric one resulting from the equation for $x_{2}$, we
obtain the effect of coupling as an accessible “corridor” in Fig. 12(b).
Similar corridors restrict the dynamics of variables $(y_{1},y_{2})$ and
$(z_{1},z_{2})$.

###  7.2 Scaling

Let us consider scaling of functions $F(x,y)$ (Eq. (15)) and $C(x_{1},x_{2})$
(Eq. (16)). The function $F(x,y)$ is scale-free. The coupling function
$C(x_{1},x_{2})$ can be written as a universal function of the normalized
variables $C(x_{1},x_{2})=c\left(\frac{x_{1}}{Q},\frac{x_{2}}{Q}\right)$.
Thus, in the PC approximation the equations for the variables
$x_{1},y_{1},z_{1},x_{2},y_{2},z_{2}$ can be renormalized

$x_{k}=Q\tilde{x}_{k},\;y_{k}=Q\tilde{y}_{k},\;,z_{k}=Q\tilde{z}_{k},\qquad
t=Q\tilde{t},\qquad(k=1,2)\;.$ (17)

In the new equations, parameter $Q$ is replaced by one, and the only remaining
parameters are $\alpha,\beta$. This is the reason why in the figures above, we
have depicted variables $x,y,z$ scaled by factor $Q$.

The scaling property of the PC model means that there is no “bifurcation
diagram” in dependence on the coupling strength $Q$ in this approximation. The
dynamics for fixed $(\alpha,\beta)$ is universal (there can be still
multistability, i.e. dependence on initial conditions). This is in a
contradistinction to observation in the original system, where the bifurcation
diagram is nontrivial (Fig. 3). However, the part of Fig. 3 at large $Q$
demonstrates predominantly chaos, and the maximal and minimal values of
variable $x_{2}$ scale $x_{max,min}\sim Q$ as the scaling (17) suggests.
According to (17), also the characteristic timescale (in particular, the
period of a periodic orbit) scales $\sim Q$. Below in this section, we use the
PC approximation with $Q=1$.

###  7.3 Synchronous cycle

Figure 13: Periodic cycle in the PC approximation. Panel (a): variables
$x(t),y(t),z(t)$. Panel (b): their derivatives
$\dot{x}(t),\dot{y}(t),\dot{z}(t)$. Markers of the time axis indicate the
instants of switches of the dynamics.

Figure 14: Projections of the trajectory Fig. 13 on phase planes
$(x,y),(y,z),(z,x)$ (red line). The numbers indicate the same instants of
switches as in Fig. 13. The black dashed lines shows borders of admissible
regions according to Fig. 10 (panels (b),(c)) and to Fig. 12(b) (panel (a)).
The blue dotted line shows the rescaled cycle for $Q=30$ in the full Eqs. (9).

As the first example of application of the PC approximation, we consider a
symmetric cycle on the manifold M1 (cf. Section 4). This cycle is described by
Eqs. (9), which in variables $(x,y,z)$ read

$\displaystyle\dot{x}$
$\displaystyle=1-\alpha\frac{1+e^{x}}{1+e^{-y}}-\beta\frac{1+e^{x}}{1+e^{-z}}\;,$
(18) $\displaystyle\dot{y}$
$\displaystyle=1-\alpha\frac{1+e^{y}}{1+e^{-z}}-\beta\frac{1+e^{y}}{1+e^{-x}}+e^{-Q}\left[\frac{(1+e^{-y})(1+e^{y})}{1+e^{-z}}-(1+e^{y})\right]\;,$
$\displaystyle\dot{z}$
$\displaystyle=1-\alpha\frac{1+e^{z}}{1+e^{-x}}-\beta\frac{1+e^{z}}{1+e^{-y}}+e^{-Q}\left[\frac{(1+e^{-z})(1+e^{z})}{1+e^{-y}}-(1+e^{z})\right]\;.$

The PC approximation of Eqs. (18) is formulated as

$\displaystyle\dot{x}$ $\displaystyle=1-\alpha F(x,y)-\beta F(x,z)\;,$ (19)
$\displaystyle\dot{y}$ $\displaystyle=1-\alpha F(y,z)-\beta F(y,x)+C(y,z)\;,$
$\displaystyle\dot{z}$ $\displaystyle=1-\alpha F(z,x)-\beta F(z,y)+C(z,y)\;.$

The solution of these equations trivially scales with parameter $Q$, but its
form depends on parameters $\alpha,\beta$. For the values of these parameters
explored above (see Eq. (12)), the cycle is presented in Figs. 13,14. Points
$0,\ldots,7$ indicate instants at which the dynamics changes; altogether the
cycle consists of 8 pieces. As Fig. 14 shows, coupling term $\sim C$ is
“working” only at stages $6\to 7\to 0$ and $2\to 3$; at other stages it is
“switched off”. Stages $1\to 2\to 3$ and $7\to 0$ follow the borders of the
affordable domains according to function $F$. In the same figures we show the
cycle in the full Eqs. (9) for $Q=30$ with blue dots; this gives an impression
of accuracy of the PC approximation.

##  8 Chaotic and periodic regimes in the PC approximation.

Figure 15: Projections of many periodic orbits of type Aa at
$\beta=1/6,\alpha=1.9$ on different subplanes. Dashed lines show available
domains according to Figs. 10,12. Figure 16: The range of possible phase
shifts $\Delta\varphi_{y}$ (red domain) and $\Delta\varphi_{z}$ (green domain)
between two subsystems for periodic orbits of type Aa. The phase shifts were
calculated as follows. For a periodic orbit of period $T$, we calculated the
times $t_{x_{1},y_{1},z_{1}}$ and $t_{x_{2},y_{2},z_{2}}$, at which the
corresponding variable attains a maximum (there is one such point on the
trajectory). After that, a phase shift in the pairs of variables was
calculated according to $\Delta\varphi_{x}=2\pi(t_{x_{2}}-t_{x_{1}})/T$,
$\Delta\varphi_{y}=2\pi(t_{y_{2}}-t_{y_{1}})/T$,
$\Delta\varphi_{z}=2\pi(t_{z_{2}}-t_{z_{1}})/T$.

In this section we report on stable periodic and chaotic regimes in the PC
model of coupled heteroclinic cycles

$\displaystyle\dot{x}_{1}$ $\displaystyle=1-\alpha F(x_{1},y_{1})-\beta
F(x_{1},z_{1})+C(x_{1},x_{2})\;,$ (20) $\displaystyle\dot{y}_{1}$
$\displaystyle=1-\alpha F(y_{1},z_{1})-\beta F(y_{1},x_{1})+C(y_{1},y_{2})\;,$
$\displaystyle\dot{z}_{1}$ $\displaystyle=1-\alpha F(z_{1},x_{1})-\beta
F(z_{1},y_{1})+C(z_{1},z_{2})\;,$ $\displaystyle\dot{x}_{2}$
$\displaystyle=1-\beta F(x_{1},y_{1})-\alpha F(x_{1},z_{1})+C(x_{2},x_{1})\;,$
$\displaystyle\dot{y}_{2}$ $\displaystyle=1-\beta F(y_{1},z_{1})-\alpha
F(y_{1},x_{1})+C(y_{2},y_{1})\;,$ $\displaystyle\dot{z}_{2}$
$\displaystyle=1-\beta F(z_{1},x_{1})-\alpha F(z_{1},y_{1})+C(z_{2},z_{1})\;.$

As mentioned above, the only relevant parameters are $\alpha$ and $\beta$.
Below we fix $\beta=1/6$ and consider the range of other parameter
$1.85\leq\alpha\leq 2.85$. The lower boundary is slightly larger than the
value $\alpha_{c}=2-\beta=1.833$ at which the heteroclinic cycle in a single
system appears. We observed 2 types of stable periodic orbits and two types of
chaos.

###  8.1 Stable periodic orbits of type A

These periodic orbits build 3 families, which are transformed to each other by
the renaming transformation (6), we denote these families as Aa, Ab, Ac. A
remarkable property of these orbits is that the phase shift between two
periodic trajectories in systems 1 and 2 is to a large extent arbitrary. We
illustrate this in Fig. 15. Here we overlap 160 periodic orbits of this type,
appearing at different random initial conditions. One can see that the
trajectories projected on the planes of variables of subsystem 1 and of
subsystem 2 coincide. However, projections on the planes showing variables of
both subsystems do not coincide, because the phase shifts are different.

Periodic orbits of type A exist for small values of $\alpha\lesssim 2.0025$.
In Fig. 16 we show the range of possible phase shifts between the two
subsystems.

We attribute the existence of synchronous periodic orbits with phase shifts in
some interval to the property of existence of dead zones in coupling, recently
introduced in [14, 15]. Indeed, in the PC model coupling between the systems 1
and 2 is only at the borders of the corridors depicted in Fig. 12(b); inside
these corridors the systems do not interact. This is exactly the property of
dead zones in the coupling terms, which, as has been shown in [14, 15], leads
to an effective coupling function with neutral intervals. Thus, phase shifts
in some range become possible.

Figure 17: Projection of a cycle of type B on different subplanes, for
$\beta=1/6,\alpha=2.4$. Dashed lines show available domains according to Figs.
10,12.

###  8.2 Stable periodic orbits of type B

Stable periodic orbit of type B exists in the whole range of explored values
of the parameter $\alpha$. This orbit is symmetric with respect to renamings,
but asymmetric with respect to transformation $1\leftrightarrow 2$. We show
projections of this trajectory on different subplanes in Fig. 17. This orbit
can be characterized as a 2:1 synchronous one, because here in system 1 there
is one maximum and one minimum per period, while in system 2 there are two
maxima and two minima. Correspondingly, the amplitude in system 2 is smaller.

###  8.3 Chaos of type A

Two types of chaotic regimes are illustrated in Fig. 18. Chaos of type A (Fig.
18(A)) inherits symmetry of the periodic cycle of type A: there are three
different attractors Aa, Ab, Ac according to renaming symmetry (6). Chaotic
regime of type B (Fig. 18(B)) is fully symmetric, both to renamings and to
exchange.

Figure 18: Projections of a trajectory on different subplanes, for
$\beta=1/6,\alpha=2.0$ (chaos of type Ax, panel (A)) and
$\beta=1/6,\alpha=2.2$ (chaos of type B, panel (B)). Dashed lines show
available domains according to Figs. 10,12.

Below we study chaos of type A at $\alpha=2$ in more details.

The starting point is a construction of a Poincaré map. As a hyperplane of
section, we have chosen $y_{1}+0.5=0$, $\dot{y}_{1}<0$. At this section, the
value of the variable $z_{1}$ is also constant, because it is slaved by
$y_{1}$. The values of the other four variables $x_{1},x_{2},y_{2},z_{2}$ are
in general different, but there are segments where one of the variables is
constant. Therefore, to resolve all distinct pieces of the attractor, we show
the obtained sequence of dots is shown in Fig. 19(a) in “skew” coordinates
$x_{2}+0.2x_{1},y_{2}-0.1x_{1},z_{2}+0.3x_{1}$. One can see that the attractor
consists of five straight intervals (the borders of these linear segments are
shown with markers and are marked with numbers). This suggests that the
attractor in the Poincaré section is one-dimensional and can be described with
a one-dimensional map.

Figure 19: Panel (a): Poincaré map and the borders of the segments. Panel (b):
Poincaré map in the one-dimensional parametrization of the curve. Dashed lines
are borders of the linear segments (marked at the axes with the same numbers
as in panel (a)).

To construct this map, we parametrize the curve according to the Euclidean
distance in the $6$-dimensional space. The total length of five segments is
$L=0.4523$. We denote the coordinate along the curve as $0\leq s\leq L$. In
Fig. 19(b) we plot the one-dimensional Poincaré map $s_{n}\to s_{n+1}$. The
points on the axes and the grid are the borders of the linear segments.
Remarkably, this map has no flat intervals, i.e. except for a finite set of
points $|\frac{ds_{n+1}}{s_{n}}|>1$. This confirms chaoticity of the dynamics.

###  8.4 Chaos of type B

Here our approach is the same as in the characterization of chaos of type A
above, the parameter now is set to $\alpha=2.2$.

Figure 20: Panel (a): Poincaré map in transformed coordinates. The borders of
linear segments are denoted with numbers from $1$ to $17$. Panel (b): Poincaré
map in the one-dimensional parametrization of the curve. Dashed lines are
borders of the linear segments.

We use the same condition $y_{1}=-0.5$, $\dot{y}_{1}<0$ for the Poincaré
section, and plot the obtained points in a three-dimensional space with
coordinates $(y_{2}-0.1x_{1},z_{2}+0.3x_{1},x_{2}+0.2x_{1})$ (the value of
$z_{1}$ at this section is constant). Now the points form not a one-
dimensional piecewise-linear curve, but a closed curve, consisting of 17
linear segments, see Fig. 20(a). This means that the attractor in the
5-dimensional Poincaré map is one-dimensional, and topologically a circle.

Next, we parametrize the attractor according to the Euclidean distance in the
$6$-dimensional space. The total length is $L=9.885$. The coordinate along the
curve is $s\in[0,L)$. Starting a trajectory on the attractor, and looking for
a first return, we obtain a one-dimensional map $s_{n}\to s_{n+1}$ plotted in
Fig. 20(b) we plot the one-dimensional Poincaré map $s_{n}\to s_{n+1}$. The
points on the axes and the grid are the borders of the linear segments (the
same numbers as in Fig. 20, point “18” is the same as “1”). This map is non-
invertable and has 6 intervals of monotonicity. On all these intervals
$|ds_{n+1}/ds_{n}|>1$, what means that chaos is robust. However, there are
regions where this derivative is close to one (e.g., in segment 9), and
presumably at the border of existence domain of chaos of type B, a stable
cycle appears when at some segments the derivative is less than one. We
cannot, however, construct the corresponding one-dimensional map for such a
range of parameters, because we cannot obtain enough iteration points to form
a closed curve like in Fig. 20(a).

##  9 Exploration of the continuous system close to the PC limit

In the previous sections, we explored the original system of coupled
heteroclinic cycles (4) and its PC model (20). It is difficult to compare
these findings because the PC model is valid in the limit $Q\to\infty$ and the
explored in Section 5 values $Q\lesssim 40$ are far from this limit.
Therefore, in this section we present simulations at rather large $Q=200$
(what corresponds to the coupling $D=1.38\cdot 10^{-87}$, at which numerical
simulations are still reliable). Below in this section we fix $\beta=1/6$ and
vary $\alpha$, to be able to compare the observed attractors with those of
Section 8.

Figure 21: Poincaré maps at section $y_{1}=-0.5;\;\dot{y}_{1}<0$. Points at
this section are shown in coordinates $x_{2},y_{2},z_{2}$. (a) $\alpha=1.86$.
Here, only stable symmetric orbits (green squares) and neutral periodic orbits
(red dots) are present. (b) $\alpha=1.93$. Here, stable symmetric orbits
(green squares), neutral periodic orbits (red dots), and small-scale chaos
(blue dots) are present. (c) $\alpha=2.08$. Here, the asymmetric large-scale
attractor (blue dots) is present. (d) $\alpha=2.70$. Here, the symmetric
large-scale attractor is present.

Numerical exploration at each value of $\alpha$ was performed as follows: for
a large set of randomly chosen initial conditions, the final state
(“attractor”) after a long transient $10^{6}$ was analyzed. The following
attractors are observed (we depict them in the Poincaré plots in Fig. 21):

Stable synchronous periodic orbits on manifolds M1,M2,M3. These attractors are
observed in the range $1.85\leq\alpha\leq 2.07$. In Figure 21(a,b) these
periodic orbits are depicted with green squares.

“Neutral periodic orbits”. These orbits, which are similar to the periodic
orbits of type A in the PC model, appear to build three families; they are
characterized by a continuous parameter. In Figure 21(a,b) these orbits are
depicted with red dots. These attractors are observed in the range
$1.85\leq\alpha\leq 1.98$. Remarkably, for larger values of parameter
$\alpha$, the range of these orbits becomes smaller. Numerical limitations do
not allow concluding, whether these orbits are true neutral ones, like the
theory of dead zones predicts [14, 15], or they just evolve so slowly that
within the explored time interval they look like stationary. To illustrate
this, we plot in Fig. 22 the evolution of these orbits together with
neighboring ones that converge eventually to real attractors - stable
symmetric periodic orbits.

Figure 22: Evolution of neutral periodic orbits (red lines) and of
trajectories close to them (green lines). The horizontal axis shows discrete
times for the Poincaré map, in real time this corresponds to $\approx 1.2\cdot
10^{6}$. One can see practically no evolution of the red trajectories in the
bulk, and convergence of green trajectories starting at the boundaries to the
stable synchronous orbits. This picture is very similar to one reported for
the oscillators coupled with a dead zone coupling (Fig 3(b1) in [15]).

Large-scale asymmetric chaos and large-scale symmetric chaos These attractors
are similar to chaos of types A and B, respectively, in the PC model, they are
observed for $\alpha>2.07$, and the transition from asymmetric to symmetric
one occurs at $\alpha\approx 2.11$. In Figure 21(c,d) these orbits are
depicted with blue dots. Figure 23(a) show the bifurcation diagram of these
attractors, obtained by slowly decreasing parameter $\alpha$ from large
values. Inside these chaotic states, there are periodic windows.

Figure 23: Panel (a): Bifurcation diagram (Poincaré maps at section
$y_{1}=-0.5;\;\dot{y}_{1}<0$) for the large-scale chaotic attractor. Panel
(b): Bifurcation diagram for the small-scale attractor. Because this branch is
isolated, to obtain it we started within chaos at $\alpha=1.97$ and varied
$\alpha$ with small steps in both directions until this attractor disappeared.
On both panels the fixed point at $y_{2}=-0.5$ is a symmetric stable periodic
orbit.

Small-scale chaos These three symmetric branches of solutions appears to be
isolated from other attractors, they are observed for $1.87\leq\alpha\leq
1.97$ and are shown in Figures 21(b) with blue dots. with blue dots.
Additionally, a bifurcation diagram for one of such attractors is shown in
Fig. 23(b).

To conclude this section, we compare the regimes in the continuous system at
$Q=200$ with those in the PC model:

  * 1.

Large scale chaos (asymmetric and symmetric) for $Q=200$ corresponds to
chaotic regimes of types A and B in the PC model. However, the ranges of
existence in dependence on parameter $\alpha$ are different.

  * 2.

The property of having a continuous set of periodic orbits in the PC model
corresponds to the observation of such periodic orbits also for $Q=200$,
although here possibly the orbits are slightly non-neutral with an extremely
slow evolution toward stable symmetric orbits. However, this slow evolution
can be hardly followed for the available length scale and for a limited
precision of calculations.

  * 3.

In the PC model the symmetric solutions on manifolds M1, M2, M3 appear
unstable, while in the continuous model for $Q=200$ they are stable in a
certain range of the parameter $\alpha$.

  * 4.

Small-scale chaos in the continuous model at $Q=200$ does not have a
counterpart in the PC model.

  * 5.

In the PC model for large $\alpha$, a stable periodic orbit is observed, while
such an orbit appears unstable (if it exists) at $Q=200$, where chaos is
observed at large values of parameter $\alpha$.

##  10 Conclusion

In this paper, we explored coupled heteroclinic cycles, rotating in opposite
directions. This makes their dynamics completely different to that of coupled
heteroclinic cycles rotating in the same direction [13], or to coupling of a
heteroclinic and a usual limit cycle [9]. In the present setup, the two
systems necessarily come to states where there is an interaction between one
large and one small variable. This makes the effect of even a very small
coupling enormous. In our numerical simulations of the full system, we
observed that chaos persist for very small coupling strengths. An increase of
coupling strength results in a series of transitions related to change of
symmetry of the chaotic attractors; finally, an inverse period-doubling
transition to a stable periodic cycle occurs.

Existence of chaos in the limit of small coupling allows for a construction of
a proper analytical model valid in this limit. The model constructed is based
on a piecewise-constant representation of r.h.s. of the equations in
transformed (logarithmic) variables. In many previous studies, such an
approximation of dynamical equations led to a possibility of an analytic
representation of the dynamics via explicit construction of a Poincaré map
[16, 17, 18, 19]. In other cases, already the basic model is formulated as a
piecewise-linear one, what enabled for an explicit construction of chaotic
solutions there [20, 21, 22]. Mostly close to our consideration is the
analysis of chaos in complex heteroclinic connections for a Bianchi IX model
in cosmology [23, 24], where an asymptotic one-dimensional map has been
constructed for a complex heteroclinic cycle. Remarkably, the piecewise-
constant model is scale-invariant: the dynamics does not depend on the
coupling parameter, what suggests that it really captures the singular limit
of zero coupling.

Unfortunately, the dynamics in the full 6-dimensional phase space in this
model appears too complex to represent it completely analytically (like in
other systems with the piecewise-constant dynamics, demonstrating chaos [25,
26, 27]). Therefore, we followed a numerical approach, where conditions and
borders of available domains were determined numerically, while the pieces of
trajectories between these turning points were drawn from the exact equations.
In this way, we constructed Poincaré maps, which appeared to be expanding non-
invertable one-dimensional maps.

A peculiar property of the interaction in the PC limit is that there are
domains of variables where the coupling vanishes exactly, so that the coupling
works only at the boundaries of these domains. Oscillatory systems with such a
property have been recently introduced by using a concept of dead zones (zones
without coupling) [14, 15]. The coupled heteroclinic cycles in the PC limit
appear to be an example of a dead zone system. This explains a rather unusual
property of existence of a family of locked cycles with phase shifts from some
interval. We further demonstrated that this property is also observed in the
continuous system, although we hypothesize that this observation may be
limited to finite although very large time intervals, because in reality the
coupling does not vanish but is exponentially small (cf. [28]).

Above, we restricted ourselves to a symmetric situation, where the parameters
inside the cycles are the same and the coupling is also symmetric. We expect
that the constructed chaotic attractor is robust under breaking of these
conditions, but this hypothesis has to be checked in future studies.

In this paper, we considered two coupled heteroclinic cycles, it appears
natural to extend this study to heteroclinic networks. We stress here that
“heteroclinic networks” can be understood in different senses. One can
consider a lattice or a network, where on each cite there is a heteroclinic
cycle, and these cycles are coupled via links with other cycles like in Eqs.
(4), see [10]. Our approach could be generalized for such a setup, if one
assumes different rotation directions at different nodes. In another
interpretation, a network consists of many nodes
$\xi_{1},\xi_{2},\xi_{3},\ldots,\xi_{n}$ with saddle connections
$\xi_{k}\to\xi_{m}$ between them [29]. In the context of our study above, one
could look at a weak coupling of two such networks, where at least some
connections are reverted.

## Acknowledgements

The authors acknowledge the partial financial support by the Israel Science
Foundation (grant No. 843/18). AP acknowledges support by the Laboratory of
Dynamical Systems and Applications NRU HSE of the Russian Ministry of Science
and Higher Education (Grant No. 075-15-2019- 1931). We thank Michael Zaks for
fruitful discussions.

## References

  * [1] R. M. May, W. J. Leonard, Nonlinear aspects of competition between three species, SIAM J. Appl. Math. 29 (2) (1975) 243–253.
  * [2] F. H. Busse, R. M. Clever, Heteroclinic cycles and phase turbulence, in: Pattern formation in continuous and coupled systems (Minneapolis, MN, 1998), Vol. 115 of IMA Vol. Math. Appl., Springer, New York, 1999, pp. 25–32.
  * [3] J. Guckenheimer, P. Holmes, Structurally stable heteroclinic cycles, Math. Proc. Cambridge Philos. Soc. 103 (1) (1988) 189–192.
  * [4] M. Krupa, Robust heteroclinic cycles, J. Nonlinear Sci. 7 (2) (1997) 129–176.
  * [5] V. S. Afraimovich, M. I. Rabinovich, P. Varona, Heteroclinic contours in neural ensembles and the winnerless competition principle, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 14 (4) (2004) 1195–1208.
  * [6] M. A. Komarov, G. V. Osipov, J. A. K. Suykens, M. I. Rabinovich, Numerical studies of slow rhythms emergence in neural microcircuits: bifurcations and stability, Chaos 19 (1) (2009) 015107, 8.
  * [7] D. Li, Y. Zhou, B. Hu, C. Zhou, Coupled perturbed heteroclinic cycles: Synchronization and dynamical behaviors of spin-torque oscillators, Phys. Rev. B 84 (2011) 104414.
  * [8] C. Bick, Heteroclinic dynamics of localized frequency synchrony: heteroclinic cycles for small populations, J. Nonlinear Sci. 29 (6) (2019) 2547–2570.
  * [9] M. Tachikawa, Specific locking in populations dynamics: Symmetry analysis for coupled heteroclinic cycles, Journal of Computational and Applied Mathematics 201 (2) (2007) 374–380, special Issue: Dynamical Systems Theory and Its Applications to Biology and Environmental Sciences.
  * [10] M. Tachikawa, Multiplicity of Limit Cycle Attractors in Coupled Heteroclinic Cycles, Progress of Theoretical Physics 109 (1) (2003) 133–138.
  * [11] M. Voit, H. Meyer-Ortmanns, Dynamics of nested, self-similar winnerless competition in time and space, Physical Review Research 1 (2) (2019) 023008.
  * [12] M. Voit, S. Veneziale, H. Meyer-Ortmanns, Coupled heteroclinic networks in disguise, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (8) (2020) 083113.
  * [13] D. Li, M. C. Cross, C. Zhou, Z. Zheng, Quasiperiodic, periodic, and slowing-down states of coupled heteroclinic cycles, Phys. Rev. E 85 (2012) 016215\.
  * [14] P. Ashwin, C. Bick, C. Poignard, State-dependent effective interactions in oscillator networks through coupling functions with dead zones, Philosophical Transactions of the Royal Society A 377 (2160) (2019) 20190042.
  * [15] P. Ashwin, C. Bick, C. Poignard, Dead zones and phase reduction of coupled oscillators, Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (9) (2021) 093132.
  * [16] A. S. Pikovsky, M. I. Rabinovich, A simple self–sustained generator with stochasic behavior, Sov. Phys. Doklady 23 (3) (1978) 183–185.
  * [17] S. V. Kijashko, A. S. Pikovsky, M. I. Rabinovich, A radio–frequency generator with stochastic behavior, Sov. J. Commun. Technol. Electronics 25 (2).
  * [18] A. S. Pikovsky, M. I. Rabinovich, Stochastic oscillations in dissipative systems, Physica D 2 (1981) 8–24.
  * [19] J. Guckenheimer, M. Wechselberger, L.-S. Young, Chaotic attractors of relaxation oscillators, Nonlinearity 19 (3) (2006) 701.
  * [20] V. N. Belykh, N. V. Barabash, I. V. Belykh, A lorenz-type attractor in a piecewise-smooth system: Rigorous results, Chaos: An Interdisciplinary Journal of Nonlinear Science 29 (10) (2019) 103108.
  * [21] V. N. Belykh, N. V. Barabash, I. V. Belykh, Bifurcations of chaotic attractors in a piecewise smooth lorenz-type system, Automation and Remote Control 81 (8) (2020) 1385–1393.
  * [22] V. N. Belykh, N. V. Barabash, I. V. Belykh, Sliding homoclinic bifurcations in a lorenz-type system: Analytic proofs, Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (4) (2021) 043117.
  * [23] O. Bogoyavlenskii, S. Novikov, Singularities of the cosmological model of the bianchi ix type according to the qualitative theory of differential equations, Zh. Eksp. Teor. Fiz 64 (1973) 1475–1494.
  * [24] O. I. Bogoyavlenskii, S. P. Novikov, Homogeneous models in general relativity and gas dynamics, Russian Mathematical Surveys 31 (5) (1976) 31.
  * [25] T. Schurmann, I. Hoffmann, The entropy of’strange’billiards inside n-simplexes, Journal of Physics A: Mathematical and General 28 (17) (1995) 5033.
  * [26] K. Peters, U. Parlitz, Hybrid systems forming strange billiards, International Journal of Bifurcation and Chaos 13 (09) (2003) 2575–2588.
  * [27] M. Blank, L. Bunimovich, Switched flow systems: pseudo billiard dynamics, Dynamical Systems 19 (4) (2004) 359–370.
  * [28] M. I. Bolotov, G. V. Osipov, A. Pikovsky, Marginal chimera state at cross-frequency locking of pulse-coupled neural networks, Phys. Rev. E 93 (2016) 032202.
  * [29] P. Ashwin, S. B. Castro, A. Lohse, Almost complete and equable heteroclinic networks, Journal of Nonlinear Science 30 (1) (2020) 1–22.

##  Appendix A Hopf bifurcation of a coexistence point on a synchronous
manifold

Below we consider the bifurcation of the limit cycle on the invariant manifold
M1 defined by relations (8), which appears due to the oscillatory instability
of the “coexistence point” $(a,a,a)$, $a=(\alpha+\beta+1)^{-1}$ (see Section
4.2).

It is convenient to transform system (9) using the following rescaling of
variables:

$t=a^{-1}\tau,\;D=ad,\;u=a(1+U),\;v=a(1+V),\;w=a(1+W).$

We obtain:

$\frac{dU}{d\tau}+U+\alpha V+\beta W=-U(U+\alpha V+\beta W),$
$\frac{dV}{d\tau}+\beta U+(1+d)V+(\alpha-d)W=-V(V+\alpha W+\beta U)$ (21)
$\frac{dW}{d\tau}+\alpha U+(\beta-d)V+(1+d)W=-W(W+\alpha U+\beta V).$

In order to derive the amplitude equation governing the evolution of small-
amplitude disturbances near the coexistence point in the vicinity of the
instability threshold point,

$d=d_{0}+\epsilon^{2}d_{2},\;\epsilon\ll 1,$ (22)

we apply the multi-scale approach. We present the variables $U$, $V$ and $W$
as functions of two variables, $\tau_{0}=\tau$ and
$\tau_{2}=\epsilon^{2}\tau$,

$U=U(\tau_{0},\tau_{2}),\;V=V(\tau_{0},\tau_{2}),\;W=W(\tau_{0},\tau_{2}),$

and replace the differentiation operator with

$\frac{d}{d\tau}=\frac{d}{d\tau_{0}}+\epsilon^{2}\frac{d}{d\tau_{2}}.$ (23)

Here $\tau_{0}$ corresponds to fast oscillations of the variables, and
$\tau_{2}$ describes the “slow” evolution of the amplitude and the phase of
oscillations near the instability threshold. The variables are presented as
the series,

$(U,V,W)=\epsilon(U^{(1)},V^{(1)},W^{(1)})+\epsilon^{2}(U^{(2)},V^{(2})+\epsilon^{3}(U^{(3)},V^{(3)},W^{(3)})+\ldots$
(24)

At the leading order we obtain the linear system

$\frac{dU^{(1)}}{d\tau_{0}}+U^{(1)}+\alpha V^{(1)}+\beta W^{(1)}=0,$
$\frac{dV^{(1)}}{d\tau_{0}}+\beta U^{(1)}+(1+d_{0})V^{(1)}+(\alpha-
d_{0})W^{(1)}=0,$ (25) $\frac{dW^{(1)}}{d\tau_{0}}+\alpha U^{(1)}+(\beta-
d_{0})V^{(1)}+(1+d_{0})W^{(1)}=0.$

The expression for the growth rates,

$\sigma_{\pm}=\frac{\alpha+\beta-2}{2}-d_{0}\pm
i\sqrt{\frac{3(\alpha-\beta)^{2}}{4}-d_{0}^{2}},$

which is equivalent to (11) up to rescaling of variables, determines the
instability threshold

$d_{0}=\frac{\alpha+\beta-2}{2}$

and the frequency of oscillations,

$\omega_{0}=\mbox{Im}\sigma=\frac{1}{2}\sqrt{3(\alpha-\beta)^{2}-(\alpha+\beta-2)^{2}}.$

The solution of (25) is

$(U^{(1)},V^{(1)},W^{1})=A(\tau_{2})e^{i\omega_{0}\tau_{0}}(U_{1}^{(1)},V_{1}^{(1)},W_{1}^{(1)})+A(\tau_{2})^{*}e^{-i\omega_{0}\tau_{0}}(U_{1}^{(1)*},V_{1}^{(1)*},W_{1}^{(1)*}),$

where $A(\tau_{2})$ is a complex amplitude function, which is still unknown.
Choosing $U_{1}^{(1)}\equiv 1$, we obtain:

$U_{1}^{(1)}=1,\;V_{1}^{(1)}=-\frac{1-\beta+i\omega_{0}}{\alpha-\beta},\;W_{1}^{(1)}=-\frac{\alpha-1-i\omega_{0}}{\alpha-\beta}.$

At the second order, we find

$(U^{(2)},V^{(2)},W^{(2)}=|A(\tau_{2})|^{2}(U_{0}^{(2)},V_{0}^{(2)},W_{0}^{(2)})+$
$A^{2}(\tau_{2})e^{2i\omega_{0}\tau_{0}}(U_{2}^{(2)},V_{2}^{(2)},W_{2}^{(2)})+A^{*2}(\tau_{2})e^{-2i\omega_{0}\tau_{0}}(U_{2}^{(2)*},V_{2}^{(2)*},W_{2}^{(2)*}),$

where

$U_{0}^{(2)}=\frac{\alpha+\beta-2}{\alpha+\beta+1}.\;V_{0}^{(2)}=-\frac{(\alpha+\beta-2)(2\beta-1)}{(\alpha+\beta+1)(\alpha-\beta)},\;W_{0}^{(2)}=\frac{(\alpha+\beta-2)(2\alpha+1)}{(\alpha+\beta+1)(\alpha-\beta)};$
$U_{2}^{(2)}=\frac{1}{12\omega_{0}^{2}}[\alpha^{2}-10\alpha\beta+\beta^{2}+8\alpha+8\beta-8-2i\omega_{0}(\alpha+\beta-2)]+\frac{1}{6\omega_{0}^{2}}\left(\frac{\alpha+\beta-2}{\alpha-\beta)}\right)^{2}[2\alpha\beta-\alpha-\beta+2i\omega_{0}(\alpha+\beta)],$
$V_{2}^{(2)}=-\frac{1}{12\omega_{0}^{2}}[2\alpha^{2}-5\alpha\beta-\beta^{2}+\alpha+7\beta-4+2i\omega_{0}(-2\alpha+\beta+1)]+$
$\frac{1}{6\omega_{0}^{2}}\left(\frac{\alpha+\beta-2}{\alpha-\beta}\right)^{2}\frac{\alpha+\beta-2-i\omega_{0})(2\alpha^{2}-7\alpha\beta+\beta^{2}+4\alpha+3\beta-3-2i\omega_{0}\beta)}{\alpha+\beta+1+2i\omega_{0}},$
$W_{2}^{(2)}=\frac{1}{12\omega_{0}^{2}}[\alpha^{2}+5\alpha\beta-2\beta^{2}-\beta-7\alpha+4-2i\omega_{0}(\alpha-2\beta+1)]+$
$\frac{1}{6\omega_{0}^{2}}\left(\frac{\alpha+\beta-2}{\alpha-\beta}\right)^{2}\frac{(\alpha+\beta-2-i\omega_{0})(\alpha^{2}-7\alpha\beta+2\beta^{2}+3\alpha+4\beta+3-2i\omega_{0}\alpha)}{\alpha+\beta+1+2i\omega_{0}}.$

At the third order in $\epsilon$, following the idea of the method of multiple
scales, we demand the solvability of the obtained system in the class of
functions periodic in $\tau_{0}$ (i.e., we eliminate the secular terms). The
solvability condition gives the amplitude equation

$\frac{dA}{d\tau_{2}}=sd_{2}A-\kappa|A|^{2}A$

that describes the slow change of the amplitude function. Here

$s=\left(\frac{\partial\sigma_{+}}{\partial
d}\right)_{d=d_{0}}=-1-i\frac{\alpha+\beta-2}{2\omega_{0}}.$

The expression for $\kappa(\alpha,\beta)$ is rather cumbersome, and we do not
write it here. Let us present that expression in the limit of small
$\delta\equiv d_{0}=(\alpha+\beta-2)/2$:

$\kappa=2i\sqrt{3}(\alpha-1)+2\delta(4-i\sqrt{3}\alpha)+O(\delta^{2}).$

One can see that for arbitrary small $\delta>0$, Re$\kappa>0$. Thus, there is
a direct Hopf bifurcation: the small stable solution with
$|A|^{2}=d_{2}\mbox{Re}s/\mbox{Re}\kappa$ exists in the region $d_{2}<0$,
where the coexistence point is unstable.

Note that in the case $\alpha+\beta=2$, $d_{0}=0$ (i.e., for system (1),
Re$s=0$, hence we do not observe the standard Hopf bifurcation. Indeed, it is
known that for $\alpha+\beta=2$, system (1) has an infinite number of periodic
solutions corresponding to oscillations with different amplitudes.