# A general method for estimating zonal transmission interface limits from
nodal network data

Patrick R. Brown<EMAIL_ADDRESS>National Renewable Energy Laboratory
Clayton P. Barrows<EMAIL_ADDRESS>National Renewable Energy
Laboratory  Jarrad G. Wright  National Renewable Energy Laboratory  Gregory L.
Brinkman National Renewable Energy Laboratory  Sourabh Dalvi  National
Renewable Energy Laboratory  Jiazi Zhang  National Renewable Energy Laboratory
Trieu Mai  National Renewable Energy Laboratory

###### Abstract

Capacity expansion models for the electric power system often employ zonal
(rather than nodal) resolution, necessitating estimates of aggregate power
transfer limits across the interfaces between model zones. Interface limits
between planning areas are sometimes published, but they are not generalizable
to arbitrary zone shapes. There is thus a need for a reproducible method for
estimating interface transfer limits (ITLs) between user-defined zones
directly from nodal transmission system data. Here, we present a simple method
for estimating ITLs using a DC power flow approximation via the power transfer
distribution factor (PTDF) matrix. Linear optimization is performed to
identify the distribution of power flows that maximizes the total flow on
interface-crossing lines, subject to individual line ratings, limits on bus
injection/withdrawal, and the relationships among flows, injections, and
withdrawals imposed by the PTDF matrix. We demonstrate the application of the
method on a 134-zone $\scriptstyle\sim$$65\,000$-bus system, and we explore
the influence of flow direction, contingency level, and zone size on the
estimated ITLs. There is significant heterogeneity in the ratio of the ITL to
the sum of interface-crossing line ratings, which highlights the importance of
accounting for the physical constraints on power flows imposed by Kirchhoff’s
laws when estimating zonal ITLs.

## Nomenclature

### Indices and sets

$z\in Z$ |  Zones
---|---
$i\in I$ |  Interfaces, where $i_{z^{\prime},z^{\prime\prime}}$ is the interface between zones $z^{\prime}$ and $z^{\prime\prime}$ (sometimes indicated by $z^{\prime}||z^{\prime\prime}$)
$b\in B$ |  Buses
$b_{v}\in B_{v}\subset B$ |  Generator buses
$b_{q}\in B_{q}\subset B$ |  Load buses
$b_{t}\in B_{t}\subset B$ |  Transmission buses
$l\in L$ |  Transmission lines and transformers (both typically referred to as “lines”), where $l^{b^{\prime},b^{\prime\prime}}$ is a line defined from starting bus $b^{\prime}$ to ending bus $b^{\prime\prime}$
$l_{i}\in L_{i}\subset L$ |  Lines crossing interface $i_{z^{\prime},z^{\prime\prime}}$

### Parameters

$p(l,b)$ |  PTDF matrix value for line $l$ and bus $b$
---|---
$r(l)$ |  Rating for line $l$ [$\mathrm{MW}\text{/}$]

### Variables

$F(l)$ |  Power flow on line $l$ [$\mathrm{MW}\text{/}$] (positive when power flows from $b^{\prime}$ to $b^{\prime\prime}$ on line $l^{b^{\prime},b^{\prime\prime}}$; negative when power flows from $b^{\prime\prime}$ to $b^{\prime}$)
---|---
$G(b)$ |  Net injection at bus $b$ [$\mathrm{MW}\text{/}$] (positive for generation, negative for load)

##  1 Introduction

Multidecadal capacity expansion models (CEMs) for interregional power system
planning typically must sacrifice spatial resolution to achieve tolerable
problem sizes and solution times [1, 2, 3]. Rather than model individual
generators, loads, and transmission lines, modelers group power systems into
geographic areas (here referred to as “zones”). The zones themselves are
treated as “copper plates”, without constraints on intrazonal power flows;
power flows across zonal interfaces are constrained by interface transfer
limits (here abbreviated as ITLs), which are meant to approximate the combined
carrying capacity of interface-crossing transmission lines under standard
operating conditions.

The selection of model zone boundaries is often constrained by the
availability of underlying input data; hourly load profiles may only be
available for balancing areas or independent system operators (ISOs),
emissions policies may apply within state boundaries, and installed generation
capacity may be available by county or state. ITLs may be published for
individual transmission corridors, as in the case of intertie paths in the
Western Electricity Coordinating Council (WECC) [4], or for the interfaces
between balancing areas or regional transmission organizations (RTOs), but
they are not generalizable to the zone shapes that may be needed to represent
other model features. Policymakers are also interested in understanding ITLs
between different planning regions and their impact on resource adequacy and
reliability [5]. Zonal ITLs necessarily represent a simplification compared to
nodal power flow and should be treated as approximate estimates rather than
precise limits applicable under all operating conditions; nevertheless, their
utility across a range of models necessitates a reproducible estimation
method.

In general, all transmission lines crossing an interface may not necessarily
be used at their rated capacities simultaneously. In large meshed networks,
the physics of electrical power flow tend to result in some lines reaching
their rated capacity before others, thus preventing the others from
contributing their full rated capacity to the ITL. Maximum power flows may
also be direction-dependent; if an outlying zone has only load buses (or
generator buses), it may only be able to accept power from (or send power to)
connected zones. Thus, ITLs often cannot be accurately approximated by the sum
of steady-state ratings of all transmission lines crossing an interface.

The convenience of quantifying interzonal transfer limits has been recognized
in planning and operations procedures for decades, and there is an associated
rich body of literature on ITL definitions [6] and calculation methods [7, 8,
9, 10]. A comprehensive review [10] compares ITL calculation methods and the
generation, load, and transmission details considered by each one. However,
the discussed methods focus on parametric evaluation of transfer capability
based on incrementally increasing power injected and withdrawn at
corresponding generation and load buses. Exceptions are presented in [8, 9],
where maximum flow problems are formulated to calculate the ITL. Nonetheless,
the maximum flow problems are formulated for individual interfaces between bus
or zone pairs and little discussion is devoted to the scalability and
generality of the presented methods. ITLs are used as the reference from which
available transfer capability (ATC) is calculated as the ITL minus the sum of
the transfer reserve margin and existing transmission commitments. The
prevalent use of ATC in electric industry planning and operations procedures
[11] suggests accurate ITL calculation is critical. Furthermore, in an
industry that increasingly relies on a distributed set of resources and is
subject to a range of policies applied over a variety of geographic regions,
ITL calculation should be both scalable and general so that ITLs can be
frequently updated and are not overly sensitive to system dispatch conditions.
Thus, there is a need for a reproducible method for estimating ITLs, obeying
the physical constraints on power flows imposed by Kirchhoff’s laws.

Here, we present a simple method for calculating ITLs between arbitrarily-
defined zones, combining individual line ratings, bus types, and a
representation of DC power flow by way of the power transfer distribution
factor (PTDF) matrix derived from the individual branch impedances and network
topology (Section 2.1). The method produces ITLs for both power flow
directions across a given interface, and it can be used at the $n-0$, $n-1$,
or higher contingency levels. We demonstrate its use on a simple 5-bus 3-zone
test system [12] (Section 2.2), and on the $\scriptstyle\sim$$65\,000$-bus
2024 system used in the North American Renewable Integration Study (NARIS)
[13] (Section 2.3) combined with the 134 default zones of the Regional Energy
Deployment System (ReEDS) model [1] (Section 3.1). Computer code for
implementing the PTDF-based method is made available open source [14], and it
can be used on both empirical (typically proprietary) network data and on
open-source synthetic test systems such as the “TAMU” dataset [15, 16].

##  2 Methods and input data

###  2.1 Calculation of interface transfer limits

####  2.1.1 Conventions

A line segment is said to cross interface $i_{z^{\prime},z^{\prime\prime}}$ if
one of its endpoints is in zone $z^{\prime}$ and its other endpoint is in zone
$z^{\prime\prime}$. Each line $l^{b^{\prime},b^{\prime\prime}}$ has a starting
bus $b^{\prime}$ and an ending bus $b^{\prime\prime}$; flow on line
$l^{b^{\prime},b^{\prime\prime}}$ is positive if power flows from bus
$b^{\prime}$ to bus $b^{\prime\prime}$ and negative if power flows from
$b^{\prime\prime}$ to $b^{\prime}$. (We refer to all network branches as
lines, even though some branches represent transformers.) For most lines, we
use the polarity of starting and ending buses defined by the underlying nodal
dataset, but interface-crossing lines are always redefined such that the
starting bus is in the first alphanumerically sorted zone and the ending bus
is in the second alphanumerically sorted zone.

####  2.1.2 Objective function

For each interface $i_{z^{\prime},z^{\prime\prime}}$ we perform independent
optimizations to determine the ITL in the “forward” direction (from
$z^{\prime}$ to $z^{\prime\prime}$) and “reverse” direction (from
$z^{\prime\prime}$ to $z^{\prime}$). In the forward direction the objective is
to maximize the sum of flows on lines crossing interface
$i_{z^{\prime},z^{\prime\prime}}$:

$\displaystyle\max_{F,G}\sum_{l_{i}}^{L_{i}}F(l_{i})$ (1a) In the reverse
direction the objective is to minimize the sum of flows on lines crossing
interface $i_{z^{\prime},z^{\prime\prime}}$, recalling that negative flows
represent flows from $z^{\prime\prime}$ to $z^{\prime}$:
$\displaystyle\min_{F,G}\sum_{l_{i}}^{L_{i}}F(l_{i})$ (1b)

These optimizations are performed independently for each interface $i$ in the
set of all interfaces $I$.

####  2.1.3 Constraints

Line ratings are enforced for all lines included in the simulation. Power can
flow along an individual line in either direction, subject to its rating:

$\displaystyle-r(l)$ $\displaystyle\leq F(l)\leq r(l)$ $\displaystyle\forall
l$ (2)

The relationship between line flows and bus injection/withdrawal is governed
by the PTDF matrix:

$\displaystyle F(l)$ $\displaystyle=\sum_{b}^{B}p(l,b)G(b)$
$\displaystyle\forall l$ (3)

The PTDF matrix provides a simple linear representation of DC power flow [17];
calculated from the reactance matrix as described in textbooks (e.g. [18]), it
relates the flows on all lines $l$ to the injections/withdrawals on all buses
$b$ in an AC network. The PTDF matrix element $p(l^{\prime},b^{\prime})$ gives
the change in power flow along line $l^{\prime}$ that accompanies a 1-unit
change in power injection from bus $b^{\prime}$.

As we are interested in the physical capacity of the transmission network
itself rather than typical operating conditions, the level of power
injection/withdrawal at each bus is not constrained by the installed
generation capacity or load participation factor. But whenever data on the bus
type is available, the bus is constrained to either inject power (for
generator buses), withdraw power (for load buses), or to neither inject nor
withdraw (for transmission buses):

$\displaystyle G(b_{v})$ $\displaystyle\geq 0$ $\displaystyle\forall b_{v}\in
B_{v}$ (4a) $\displaystyle G(b_{q})$ $\displaystyle\leq 0$
$\displaystyle\forall b_{q}\in B_{q}$ (4b) $\displaystyle G(b_{t})$
$\displaystyle=0$ $\displaystyle\forall b_{t}\in B_{t}$ (4c)

If no information on bus type is available for a given bus, or if a bus
qualifies as both a generator and load bus, constraint 4 is not applied for
that bus.

Together, constraints 3 and 4 tend to limit the ITL to a value less than the
simple sum of interface-crossing line ratings.

###  2.2 Example using 5-bus test system

We first illustrate the method for a simple 5-bus test system [12], with
topology and line parameters (reactance and power rating) defined in Table 1.

Table 1: 5-bus test system line parameters. Line |  From bus |  To bus |  From zone |  To zone |  Reactance |  Rating [MW]
---|---|---|---|---|---|---
A$|$B |  A |  B |  1 |  1 |  0.0281 |  400
B$|$C |  B |  C |  1 |  1 |  0.0108 |  400
C$|$D |  C |  D |  1 |  2 |  0.0297 |  400
D$|$E |  D |  E |  2 |  3 |  0.0297 |  240
A$|$E |  A |  E |  1 |  3 |  0.0064 |  400
A$|$D |  A |  D |  1 |  2 |  0.0304 |  400

The corresponding PTDF matrix for this network is provided in Table 2, using
bus A as the slack bus (the results are independent of the choice of slack
bus). This example does not include limitations on injection/withdrawal by bus
type; each bus is allowed to act as either generation or load, so constraint 4
is inactive and the ITL for each interface is the same in both directions.

Table 2: PTDF matrix for the 5-bus test system. | A | B | C | D | E
---|---|---|---|---|---
A$|$B | 0 | -0.6698 | -0.5429 | -0.1939 | -0.0344
B$|$C | 0 | 0.3302 | -0.5429 | -0.1939 | -0.0344
C$|$D | 0 | 0.3302 | 0.4571 | -0.1939 | -0.0344
D$|$E | 0 | 0.1509 | 0.2090 | 0.3685 | -0.1120
A$|$E | 0 | -0.1509 | -0.2090 | -0.3685 | -0.8880
A$|$D | 0 | -0.1792 | -0.2481 | -0.4376 | -0.0776

The resulting ITLs are given in Table 3. Each ITL is associated with an
independent set of line flows, with flows corresponding to the ITL of the
$1||2$ interface shown in Figure 1.

Table 3: 5-bus test system transfer limits. Interface |  Interface-crossing lines |  Sum of ratings [MW] |  ITL [MW]
---|---|---|---
1$||$2 |  A$|$D, C$|$D |  800 |  719
2$||$3 |  D$|$E |  240 |  400
1$||$3 |  A$|$E |  400 |  240
Figure 1:  Line flows in the 5-bus test system corresponding with the ITL of
the $\bm{1||2}$ interface. The first number in each black or green pair gives
the flow on the indicated line; the second number gives the line rating. Lines
crossing the $1||2$ interface are shown in green; lines not crossing the
$1||2$ interface are shown in black. The width of the gray and light green
bars indicates the line rating; the width and color of the overlaid arrow
indicate the flow from source (wide end of arrow) to sink (narrow end of
arrow). Blue numbers indicate power injection (+) and red numbers indicate
power withdrawal (–) at each bus.

Note that while the total flow from Zone 1 to Zone 2 (including flow through
Zone 3) is $959\text{\,}\mathrm{MW}\text{/}$, the ITL only includes the
contribution from interface-crossing lines A$|$D and C$|$D. In this case, line
C$|$D is fully utilized under interface-flow-maximizing conditions, but line
A$|$D is underutilized, contributing $319\text{\,}\mathrm{MW}\text{/}$ of a
potential $400\text{\,}\mathrm{MW}\text{/}$ to the ITL. As a result of the
constraints imposed by the PTDF matrix (i.e. the impact of Kirchhoff’s voltage
law on power flows), increasing the flow on line A$|$D would require
increasing the flow on lines A$|$E (in the positive direction) and D$|$E (in
the negative direction, from bus E to bus D), both of which are already at
their rated limit. Increasing the ITL between Zones 1 and 2 could thus be
accomplished by increasing the rating of line C$|$D; changing the reactance of
lines A$|$D, A$|$E, and D$|$E; or increasing the ratings of lines A$|$E and
D$|$E; but not by increasing the rating of interface-crossing line A$|$D.

###  2.3 Demonstration data for contiguous U.S.

The lines and buses included in the NARIS nodal dataset [13] are shown in
Figure 2 and Figure 3 respectively.

Figure 2:  Transmission line dataset for the contiguous U.S. The rating of
each transmission line is indicated by the line width; the voltage is
indicated by the color. The data shown here were assembled as part of the
North American Renewable Integration Study (NARIS) [13].  Figure 3:  Bus
dataset for the contiguous U.S. [13]

The following data-cleaning steps were employed:

  * •

Buses without location data were assigned to the location of the bus separated
by the smallest number of “hops” through the network. If a bus without
location data is connected to two or more buses by the same number of network
hops (e.g. if bus A, without location data, is connected to buses B and C,
which do have location data), the missing-location bus was assigned to a point
located at the centroid of the connected buses with location data.

  * •

Buses located more than $100\text{\,}\mathrm{km}\text{/}$ outside the
contiguous U.S. were dropped from the dataset.

  * •

Network reduction was performed to remove radial buses connected to only a
single other bus, reducing problem size. For the purposes of constraint 4, if
a generator bus was removed, the bus it was attached to is treated as a
generator; if a load bus was removed, the bus it was attached to is treated as
a load.

  * •

DC lines and back-to-back interties were dropped, as the method described here
only applies to AC transmission. For application in a CEM or production cost
model (PCM), DC lines and interties would be modeled separately, with flow
limits defined by their rated capacities.

  * •

For lines without known rating data, missing ratings were inferred using
practical capability limits for transmission lines based on the loadability
concepts introduced by St. Clair [19] (Figure S1 in the Supplemental
Information). Normalized line characteristics at typical subtransmission and
transmission voltages were based on [20] while calculations for thermal and
surge impedance loading (SIL) ratings were based on [21]. Line lengths were
calculated as the geodesic distance between the two line end points, likely
underestimating the actual length and hence slightly overestimating actual
ratings.

  * •

For lines without reactance data, the missing reactance was copied from the
line of the same voltage with the length nearest to the length of the line
without reactance data. For transformers without reactance data (which have
zero length), the missing reactance was taken as the median across
transformers of the same voltage that include reactance data.

  * •

The three interconnections (the Eastern Interconnection, Western
Interconnection, and Texas Interconnection) were simulated independently. When
assessing interfaces within the Eastern Interconnection, only buses within
$800\text{\,}\mathrm{km}\text{/}$ of the specified interface were included in
the simulation to reduce problem size. For interfaces within the Western
Interconnection and Texas Interconnection, all buses within the
interconnection were included.

The ITL calculation described here can be applied to arbitrary zone shapes
ranging from political boundaries (counties or states) to actual balancing
areas. For demonstration purposes, we apply it to the 134 default zones used
in the ReEDS model, shown in Figure 4 [1].

Figure 4:  Default model zones used in ReEDS. Adjacent zones within the same
interconnection are connected by red lines.

For each pair of zones within the same interconnection connected by at least
one line, we first calculate the ITL in both directions (i.e. from the
alphanumerically first zone to the second zone, and vice versa), including all
lines in the network (henceforth referred to as “$n-0$” conditions). For each
interface and direction, we then identify the interface-crossing line with the
highest flow, remove it from the network, and recalculate the ITL. These
(typically lower) ratings are referred to as “$n-1$” ratings. Because
intrazonal, non-interface-crossing lines could also be binding, a true
contingency analysis would repeat the calculation for every line in the
network, not just interface-crossing lines. As this approach would increase
computation times by many orders of magnitude, we do not employ it here;
instead we focus our contingency analysis on interface-crossing lines alone.

###  2.4 Caveats and limitations

Before discussing the resulting ITLs for the ReEDS interfaces, we first note
some limitations of the method described here.

  * •

This method does not consider the use of phase-shifting transformers or other
flexible AC transmission system (FACTS) devices that can modify power flow
over individual branches or groups of branches. Increased use of FACTS devices
would relax constraint 3 and could increase ITLs.

  * •

As we are concerned with the physical capabilities of the transmission
network, we do not place limits on the level of injections/withdrawals to or
from individual buses beyond constraint 4. For the existing system, injections
would be limited by the amount of installed generation capacity (combined with
weather conditions for variable renewable energy resources), and withdrawals
would likely be distributed according to historically observed load
participation factors rather than varying arbitrarily in order to enable
maximum interface flows. Including limits based on installed generation
capacity and observed withdrawals (load participation factors) would be
expected to reduce the calculated ITLs.

  * •

The method used to infer missing branch ratings (described in the Supplemental
Information) considers the combination of individual branch thermal and
angular stability limits but not voltage stability limits. This approach
provides an upper limit on individual branch ratings as a function of line
length. Branch ratings reported in the NARIS dataset are used as-is, assuming
the rating includes applicable de-rating beyond thermal limits. It is expected
that a majority of individual branch ratings would be based on thermal limits,
with exceptions for long lines that could include appropriate de-ratings
relative to thermal limits. If a reported branch rating only includes the
thermal limit and neglects applicable de-rating, the resulting ITL could be
overestimated.

  * •

Zonal capacity-expansion or production-cost simulations using interzonal ITLs
offer lower resolution than a full nodal model since zonal models treat the
area within zones as congestion-free “copper plates.” When possible, it is
preferable to locate zone boundaries such that they intersect important
congestion-limited corridors; otherwise, the ITL could give an inflated
estimate of the actual ability to move power across a zone with the existing
network. For example, if a zone were defined such that it encompassed the
boundary between two asynchronous interconnections, the ability to flow power
across that zone between interconnections would be limited not only by the
ITLs on either side of the zone but also by the capacity of any back-to-back
interties located within the zone. Thus, while the method described here can
be used for any specified zone boundaries, the utility of the estimated ITL
will depend on how well the specified zones intersect with congestion patterns
of interest.

##  3 Results

###  3.1 ITLs resolved by flow direction and contingency level

Figure 5 shows the difference between the sum of line ratings, the ITL at
$n-0$ contingency level, and the ITL at $n-1$ contingency level for the 293
synchronous zonal interfaces associated with the 134 ReEDS model zones.

Figure 5:  ITLs for ReEDS zones as a function of direction and contingency
level. a, Map showing the sum of ratings for interface-crossing lines (gray),
the ITL at $n-0$ contingency level (red), and the ITL at $n-1$ contingency
level (black). Each interface is represented by a group of bars drawn between
the centroids of the two zones it connects, with the thickness indicating the
rating; ITLs are shown for each direction, with the ITL flowing into a zone
indicated by the thickness of the half of the bar pointing into that zone. b,
Distribution of the ratios of the ITL between the higher-limit and lower-limit
directions for each interface, for both $n-0$ (red) and $n-1$ contingency
levels (black). c, Distribution of the ratios between the ITL and the sum of
ratings of interface-crossing lines for that interface. Red lines indicate
$n-0$ conditions; black lines indicate $n-1$ conditions. Solid lines indicate
ratios calculated using the higher-capacity direction; dashed lines indicate
ratios calculated using the lower-capacity direction. The distribution on the
$x$ axis includes 293 interfaces between the 134 ReEDS model zones that are
connected by at least one AC transmission line.

There is substantial heterogeneity in the impact of both the direction of flow
and the contingency level on ITLs across the zonal interfaces assessed here.
Under $n-0$ conditions, $31\text{\,}\mathrm{\char 37\relax}\text{/}$ of
interfaces have ITLs that vary between the two directions of flow by
$\leq$5\text{\,}\mathrm{\char 37\relax}\text{/}$$ , while
$6\text{\,}\mathrm{\char 37\relax}\text{/}$ of interfaces have ITLs that vary
between the two directions by more than a factor of 2 (red trace in Figure
5b). In the higher-flow direction, $16\text{\,}\mathrm{\char
37\relax}\text{/}$ of interfaces have an ITL equal to the sum of ratings for
interface-crossing lines, while $12\text{\,}\mathrm{\char 37\relax}\text{/}$
of interfaces have an ITL equivalent to less than half of the sum of ratings
for interface-crossing lines (red solid trace in Figure 5c). The median ratio
of the ITL to the sum of interface-crossing line ratings at the $n-0$
contingency level is $84\text{\,}\mathrm{\char 37\relax}\text{/}$ for the
higher-ITL direction and $69\text{\,}\mathrm{\char 37\relax}\text{/}$ for the
lower-ITL direction.

Under $n-1$ conditions, the ITL is typically substantially smaller than the
sum of interface-crossing line ratings: the median ratio is
$44\text{\,}\mathrm{\char 37\relax}\text{/}$ for the higher-ITL direction and
$33\text{\,}\mathrm{\char 37\relax}\text{/}$ for the lower-ITL direction.
$12\text{\,}\mathrm{\char 37\relax}\text{/}$ of interfaces have an ITL of zero
in at least one of the two directions (black trace in Figure 5b).

Figure 6 shows the impact of applying $n-1$ contingency conditions relative to
the $n-0$ ITL.

Figure 6:  Distribution of the impacts of subtracting the line with largest
flow from each interface, shown for the direction with lower $\bm{n-0}$ rating
(blue) and the direction with higher $\bm{n-0}$ rating (orange). a, Difference
between $n-1$ and $n-0$ ratings. b, Ratio between $n-1$ and $n-0$ ratings. The
distribution includes 293 interfaces between the 134 ReEDS model zones that
are connected by at least one AC transmission line.

The median reduction in ITL for $n-1$ conditions is
0.9–$1.1\text{\,}\mathrm{GW}\text{/}$, or 39–$42\text{\,}\mathrm{\char
37\relax}\text{/}$, depending on the direction of flow. As noted above, the
reduction is much larger for a small collection of interfaces; 9–10 interfaces
(depending on the flow direction) lose $\geq$4\text{\,}\mathrm{GW}\text{/}$$
of ITL capacity at the $n-1$ contingency level.

###  3.2 Aggregating transfer limits of smaller zones to obtain transfer
limits for larger zones

Accurate ITLs are likely to be useful at a variety of geographic resolutions.
ITLs for smaller zones can be useful for assessing the trade-off between
installing renewable energy capacity in high-quality resource regions far from
load (incurring transmission expansion costs) and installing capacity in
lower-quality resource regions nearer to load; ITLs between larger planning
regions (such as ISOs or RTOs) may be useful for longer-term interregional
resource planning. Having calculated the ITLs between the 134 ReEDS regions,
it is natural to ask whether these ITLs may simply be summed to obtain ITLs
between larger regions, or whether the ITLs should be recalculated directly
including all interface-crossing lines for the larger regions.

Here we consider 11 “planning regions” that approximate the FERC Order 1000
planning regions [22] and ERCOT. These planning regions consist of groupings
of the 134 ReEDS model zones (Figure 7a).

Figure 7:  Dependence of planning-region ITLs on calculation method. a, Map of
11 planning regions, approximating FERC Order 1000 planning regions [22] and
ERCOT. Black lines indicate the boundaries between synchronous
interconnections. b, ITLs between planning regions as a function of
calculation method, direction, and contingency level. Green values were
directly calculated from flows on lines that cross the interfaces between
planning regions; purple values were obtained by summing the ITLs of the
underlying ReEDS zones; gray bars indicate the sum of ratings of interface-
crossing lines. $n-1$ capacities are shown by lighter bars; additional
capacity under $n-0$ conditions is shown by darker bars.

Figure 7b shows the ITLs between these planning regions calculated using two
different methods. Values indicated by green bars were calculated directly
from the flows on all lines crossing the interface between a given pair of
planning regions. Values indicated by purple bars are sums over the relevant
ITLs between ReEDS zones shown in Figures 4 and 5; for example, the
“ISONE$\rightarrow$NYISO” ITL under $n-1$ conditions is the sum of the
p129$\rightarrow$p127, p131$\rightarrow$p127, and p132$\rightarrow$p127 ITLs
under $n-1$ conditions.

Two consistent trends are apparent in Figure 7b. First, $n-0$ ITLs are
uniformly higher when calculated by aggregating the ITLs of constituent zones
than when calculated directly. This trend is directly analogous to the trend
observed in Figure 5c; just as the sum of ratings of interface-crossing lines
tends to overestimate the ITL, the sum of ITLs of constituent regions tends to
overestimate the ITL between larger regions. The difference is larger for
planning region interfaces that encompass more constituent-zone interfaces;
for example, the interfaces with MISO have a larger discrepancy between the
two methods than the interfaces with NYISO.

Second, the $n-1$ ITLs are consistently lower when obtained from the sum of
$n-1$ ITLs for constituent zones than when calculated directly. The
MISO$||$SPP interface encompasses 18 zonal interfaces, so the sum of those 18
$n-1$ capacities corresponds to a situation where 18 interface-crossing lines
are not contributing to the ITL. When the MISO$||$SPP $n-1$ ITL is calculated
directly, only a single interface-crossing line is dropped. As such, if an
$n-1$ ITL is desired between larger aggregated regions, it should be
calculated directly (as represented by the green bars in Figure 7) rather than
as a sum of the $n-1$ ITLs of constituent interfaces (as represented by the
purple bars in Figure 7).

##  4 Conclusions

In summary, we have demonstrated a scalable method for calculating zonal
interface transfer limits from nodal network transmission data that is
applicable to arbitrary zone shapes and is compatible with different
contingency levels. Results for the $\scriptstyle\sim$$65\,000$-bus NARIS
network (Figure 5) illustrate that while summing the ratings of interface-
crossing lines provides a reasonable estimate of $n-0$ ITLs for some
interfaces, it overestimates the ITL for the majority of interfaces. In
general, accounting for the limitations imposed by line reactances and bus
injection/withdrawal capabilities helps improve the accuracy of estimated
interface capacities and reduces the divergence between zonal models and more
detailed nodal models.

The method presented here could be extended in several ways. Comparing
estimated ITLs to published path or interface ratings could assist in cleaning
of the underlying network data and improve the quality of the estimated ITLs.
The dual values of constraint 2 provide the marginal increase in ITL that
would result from increasing the rating of a particular line; if combined with
estimates of upgrade costs (e.g. via reconductoring), one could construct an
“upgrade supply curve” that identifies the most cost-effective interventions
for achieving a desired increase in ITL. Modifying network parameters via
FACTS could relax constraint 3, and adding energy storage or distributed
energy resources (DERs) at buses without existing generation assets could
alleviate constraint 4; both strategies could accordingly increase ITLs.
Adopting dynamic line ratings (DLR) could also increase ITLs but would
increase complexity, given the need to repeat the ITL calculation for
different weather conditions.

The greatest limitation of this method is its reliance on detailed nodal
network data, which typically are not available publicly. Necessary data
fields for the application of this method include line reactance, line power
rating, bus location, and connectivity between lines and buses; bus type, if
available, can improve the estimate by allowing the use of constraint 4. Some
open-source synthetic network datasets are available (e.g. [15, 16]);
comparing ITLs calculated from actual network data with those calculated from
synthetic networks could provide a useful validation metric for improving the
quality of synthetic datasets.

## Acknowledgments

The authors acknowledge Luke Lavin, Hamody Hindi, Josh Novacheck, Leonardo
Rese, and Mike Meshek for helpful discussions and feedback. This work was
authored by the National Renewable Energy Laboratory, operated by the Alliance
for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under
contract no. DE-AC36-08GO28308. Funding was provided by the U.S. Department of
Energy Office of Energy Efficiency and Renewable Energy Wind Energy
Technologies Office. The views expressed in the article do not necessarily
represent the views of the DOE or the U.S. Government.

## References

  * [1] Ho, J., Becker, J., Brown, M., Brown, P., Chernyakhovskiy, I., Cohen, S., Cole, W., Corcoran, S., Eurek, K., Frazier, W., et al. (2021).  Regional Energy Deployment System (ReEDS) Model Documentation: Version 2020.  National Renewable Energy Laboratory (NREL) https://www.nrel.gov/docs/fy21osti/78195.pdf.
  * [2] Young, D. (2020).  US-REGEN Model Documentation.  Electric Power Research Institute.  https://www.epri.com/research/products/000000003002016601.
  * [3] United States Environmental Protection Agency (2021).  Documentation for EPA’s Power Sector Modeling Platform v6 Using the Integrated Planning Model.  United States Environmental Protection Agency.  https://www.epa.gov/system/files/documents/2021-09/epa-platform-v6-summer-2021-reference-case-09-11-21-v6.pdf.
  * [4] Western Electricity Coordinating Council (2022).  2022 Path Rating Catalog - Public Version.  https://www.wecc.org/Reliability/2022PathRatingCatalogPublic.pdf.
  * [5] U.S. Federal Energy Regulatory Commission (2022).  Staff-Led Workshop on Establishing Interregional Transfer Capability Transmission Planning and Cost Allocation Requirements.  https://www.ferc.gov/news-events/events/staff-led-workshop-establishing-interregional-transfer-capability-transmission.
  * [6] Ilić, M., Galiana, F., Fink, L., Bose, A., Mallet, P., and Othman, H. (1998).  Transmission capacity in power networks.  International Journal of Electrical Power & Energy Systems _20_ , 99–110.  http://dx.doi.org/10.1016/S0142-0615(97)00074-4.
  * [7] Austria, R., Chao, X., Reppen, N., and Welsh, D. (1995).  Integrated approach to transfer limit calculations.  IEEE Computer Applications in Power _8_ , 48–52.  http://dx.doi.org/10.1109/67.363302.  Conference Name: IEEE Computer Applications in Power.
  * [8] Gan, D., Thomas, R., and Zimmerman, R. (2000).  Stability-constrained optimal power flow.  IEEE Transactions on Power Systems _15_ , 535–540.  http://dx.doi.org/10.1109/59.867137.  Conference Name: IEEE Transactions on Power Systems.
  * [9] Gan, D., Luo, X., Bourcier, D.V., and Thomas, R.J. (2003).  Min–max transfer capabilities of transmission interfaces.  International Journal of Electrical Power & Energy Systems _25_ , 347–353.  http://dx.doi.org/10.1016/S0142-0615(02)00099-6.
  * [10] Mohammed, O.O., Mustafa, M.W., Mohammed, D.S.S., and Otuoze, A.O. (2019).  Available transfer capability calculation methods: A comprehensive review.  International Transactions on Electrical Energy Systems _29_ , e2846.  http://dx.doi.org/10.1002/2050-7038.2846.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/2050-7038.2846.
  * [11] PJM - Interregional Planning.  https://www.pjm.com/planning/interregional-planning.
  * [12] Li, F. and Bo, R. (2010).  Small test systems for power system economic studies.  IEEE PES General Meeting, PES 2010 , 16–19http://dx.doi.org/10.1109/PES.2010.5589973.
  * [13] Brinkman, G., Bain, D., Buster, G., Draxl, C., Das, P., Ho, J., Ibanez, E., Jones, R., Koebrich, S., Murphy, S., et al. (2021).  The North American Renewable Integration Study (NARIS): A U.S. Perspective.  NREL.  https://www.nrel.gov/docs/fy21osti/79224.pdf.
  * [14] Barrows, C.P. (2023).  InterfaceLimits.jl.  https://github.com/NREL-Sienna/InterfaceLimits.jl.
  * [15] University, T.A. (2022).  Texas A&M University Electric Grid Datasets.  https://electricgrids.engr.tamu.edu/.
  * [16] Xu, Y., Myhrvold, N., Sivam, D., Mueller, K., Olsen, D.J., Xia, B., Livengood, D., Hunt, V., D’Orfeuil, B.R., Muldrew, D., et al. (2020).  U.S. test system with high spatial and temporal resolution for renewable integration studies.  IEEE Power and Energy Society General Meeting _2020-Augus_.  http://dx.doi.org/10.1109/PESGM41954.2020.9281850.
  * [17] Christie, R.D. and Wollenberg, B.F. (2000).  Transmission Management in the Deregulated Environment.  Proceedings of the IEEE _88_.  http://dx.doi.org/10.1109/5.823997.
  * [18] Wood, A.J., Wollenberg, B.F., and Sheble’, G.B. (2014).  Power generation, operation, and control (John Wiley & Sons, Inc., Hoboken, NJ), 3 edition.  https://handoutset.com/wp-content/uploads/2022/03/Power-Generation-Operation-and-Control-3rd-edition-by-Allen-J.-Wood-Bruce-F.-Wollenberg-Gerald-B.-Shebl.pdf.
  * [19] St. Clair, H.P. (1953).  Practical concepts in capability and performance of transmission lines.  Transactions of the American Institute of Electrical Engineers. Part III: Power Apparatus and Systems _72_ , 1152–1157.  http://dx.doi.org/10.1109/AIEEPAS.1953.4498751.
  * [20] Hao, J. and Xu, W. (2008).  Extended transmission line loadability curve by including voltage stability constrains.  In 2008 IEEE Canada Electric Power Conference. pp. 1–5.  http://dx.doi.org/10.1109/EPC.2008.4763379.
  * [21] Saadat, H. (1999).  Power System Analysis (McGraw-Hill), 1st edition.
  * [22] U.S. Federal Energy Regulatory Commission (2021).  Order No. 1000 - Transmission planning and cost allocation.  https://www.ferc.gov/electric-transmission/order-no-1000-transmission-planning-and-cost-allocation.

## Supplemental Information

Figure S1 shows the line loadability curves used to fill missing line ratings
as discussed in Section 2.3. Shorter lines are limited by their thermal rating
(current carrying capability) while longer lines are limited by their
stability limits. For the purposes of this analysis, selected simplifications
in the representation of lines were made to characterize these loadability
curves and include:

  * •

Lossless lines (simplified representation of receiving-end to sending-end
voltages as well as phase-angles).

  * •

A maximum voltage drop of $5\text{\,}\mathrm{\char 37\relax}\text{/}$ from
receiving end to sending end.

  * •

A maximum power angle of $45\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$
between the receiving end and sending end (to ensure a reasonable stability
margin). Although the absolute theoretical power angle limit under steady-
state conditions is $90\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$, ranges
of $35\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$–
$45\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$ are typical to ensure
transient stability during large system disturbances.

  * •

Voltage specific lumped model parameters from [20].

  * •

The exclusion of voltage stability limits (only angular stability limits are
included).

Figure S1:  Line loadability versus distance used to fill missing line
ratings.


# Labyrinth Chaos: Revisiting the Elegant, Chaotic and Hyperchaotic Walks

Vasileios Basios  Service de Physique des Systèmes Complexes et Mécanique
Statistique, Université Libre de Bruxelles, Belgium  Interdisciplinary Center
for Nonlinear Phenomena and Complex Systems (CeNoLi), Belgium††thanks: E-mail:
<EMAIL_ADDRESS>Chris G. Antonopoulos  Department of Mathematical Sciences,
University of Essex, UK††thanks: E-mail<EMAIL_ADDRESS>Anouchah Latifi
Department of Physics, Faculty of Sciences, Qom University of Technology,
Iran††thanks: E-mail<EMAIL_ADDRESS>
###### Abstract

Labyrinth chaos was discovered by Otto Rössler and René Thomas in their
endeavour to identify the necessary mathematical conditions for the appearance
of chaotic and hyperchaotic motion in continuous flows. Here, we celebrate
their discovery by considering a single labyrinth walks system and an array of
coupled labyrinth chaos systems that exhibit complex, chaotic behaviour,
reminiscent of chimera-like states, a peculiar synchronisation phenomenon. We
discuss the properties of the single labyrinth walks system and review the
ability of coupled labyrinth chaos systems to exhibit chimera-like states due
to the unique properties of their space-filling, chaotic trajectories, what
amounts to elegant, hyperchaotic walks. Finally, we discuss further
implications in relation to the labyrinth walks system by showing that even
though it is volume-preserving, it is not force-conservative.

Keywords: Thomas-Rössler systems, Labyrinth chaos, Labyrinth walks, Chaos,
Hyperchaos, Force-conservative dynamics, Volume-preserving systems

In the course of their pioneering work on the properties of feedback circuits
and their relation to chaos and hyperchaos, Otto Rössler and René Thomas
proposed a minimal model of a dynamical system they termed Labyrinth Chaos. It
turned out that even though it is simple, it is full of surprising properties
that no-one could thought of. Simple and elegant as it is, it still holds
great promise in elucidating aspects of chaotic dynamics that are not evident
in other systems. Our paper revisits their work and highlights the incredible
riches of this system in its disconcerting simplicity and importance in the
context of dynamical systems and in other fields.

##  1 Introduction

Undoubtedly, two of the most celebrated, low-dimensional, chaotic systems in
the theory of dynamical systems are the Rössler and Lorenz systems with their
iconic Rössler and Lorenz attractors [1, 2]. Rössler’s work to elucidate the
fundamentals of the topology of chaos, led him to investigate expanding and
contracting circuits for his famous folding mechanism. His construction
remains one of the simplest and most elegant chaotic dynamical systems to
date, the Rössler system [3].

Rössler’s imperative to focus on the role of feedback, bore fruits quickly:
For the first time, the Rössler system provided a clear demonstration of one
of the two scenaria for the Shilnikov-route to chaos, the period doubling
cascade [4]. In his work, in bringing forward the concept of feedback circuits
in dynamical systems theory, he found a great ally and coworker, Thomas.
Thomas was already a renown biologist and deeply familiar with the idea of
feedback circuits and their Boolean-logic aspects, present in phenomena
associated with gene-expression, homeostasis, multi-stationarity and memory.
He joined forces with his friend Rössler and turned him to an “amateur
mathematician”, as he would put it! Rössler, wrote in [5] about his attractor
and the role of feedback circuits: “Later, René Thomas saw much deeper into
the topology of this feedback circuit”. However, the truth is that together,
they managed to see so deep that they discovered many things which are
fundamental in chaotic dynamics. Among them, was the importance of the
coexistence of positive and negative circuits as the necessary mathematical
conditions for the appearance of chaos and hyperchaos, in the sense of the
presence of more than one positive Lyapunov exponents. In this context, the
term “circuits” was used by Rössler and Thomas to refer to those terms in the
Jacobian matrix of the dynamical system, whose row and column indices are in
circular permutation. Therefore, circuits can be positive or negative
according to the sign of the product of their terms as they appear in the
determinant of the Jacobian matrix [6, 7] of the system.

In their venture [6, 7, 8], Rössler and Thomas proposed a class of continuous-
flow dynamical systems, called Thomas-Rössler (TR) systems [9]. Their motive
was to investigate what kind of nonlinearity is necessary to support chaos and
hyperchaos in the simplest possible way. This led them to focus on the
necessary mathematical conditions for the appearance of chaos, in terms of
positive and negative circuits.

Their pioneering work on TR systems provided the ground for a plethora of
published works of biological relevance (e.g. the special issue in memory of
Thomas [10], and contributions therein) and of a more general scope in complex
systems [11, 10], especially related to the emergence of complex behaviour via
simple-circuit feedback structures in complex-systems science at large.

The main results of their work can be summarised in the following three
statements: (i) a positive circuit is necessary for the system to have stable
states, (ii) a negative circuit is necessary for the system to exhibit,
robust, sustained oscillations, and (iii) a necessary condition for chaos is
the presence of both a positive and a negative circuit in the system (for a
more complete analysis, one may have a look at [10], where the authors discuss
the subject in more details). However, concerning hyperchaos, we can add a
fourth statement: (iv) a necessary condition for hyperchaos is the presence of
more than one positive and negative, circuits.

What Rössler and Thomas called labyrinth chaos and “Arabesques” [7, 12, 9]
drew wide attention as examples of elegant chaos, also discussed in [3]. Here,
the term “elegant” is understood in the sense of Poincaré’s demand for beauty
and simplicity in expressing mathematical ideas, something that Rössler and
Thomas’ work surely adheres to.

Interestingly, the class of TR systems exhibit typical properties, such as
complex periodicity, multi-stationarity, symmetries, etc. Notably, for
appropriate parameter ranges, there typically appear coexisting, symmetric,
chaotic attractors. However, one of their most spectacular findings is that of
the special case of a totally new type of chaos in the absence of attractors:
a chaotic, volume-preserving dynamical system with no attractors, the so-
called labyrinth walks system. The existence of chaos in the absence of
chaotic attractors is due to the existence of only unstable fixed points,
located on an infinitely large, periodic lattice in a volume-preserving state-
space. Hence, the trajectories wander through the lattice of unstable fixed
points, exhibiting one or more positive Lyapunov exponents (i.e. hyperchaotic
behaviour), depending on the dimensionality of the system. They termed these
wandering, space filling trajectories, labyrinth walks [6, 7, 8], in the sense
they exhibit super-diffusive, fractional Brownian properties, which can also
be described by complex, symbolic dynamics [13].

What has also been attributed to Rössler since the early dates of the
development of chaos theory, is the term “hyperchaos” [14], drawing attention
to the topological idiosyncrasies of chaos, what refers to dynamics with more
than one exponentially expanding directions, i.e. to more than one positive
Lyapunov exponents. Ever since, its special significance has been studied in a
plethora of works, ranging from fractional-order dynamical systems [15] to
coupled quantum systems [16]. His contribution was and still remains seminal
as certain fundamental questions, especially about the topology of chaos, are
still open today [17, 18, 2]. Rössler and co-workers have also demonstrated
that labyrinth chaos and labyrinth walks are hyperchaotic systems in more than
three dimensions, showing the truth of condition (iv) as the necessary
condition for the existence of hyperchaos [7].

Amazingly, these systems still continue to surprise us with their elegance,
beauty and combination of properties: As it has been reported recently in [9],
arrays of coupled labyrinth chaos systems exhibit stereotypical, chimera-like
states, reminiscent of chimera states observed in coupled Kuramoto oscillators
[19, 20, 21, 9] and in other systems [22, 23, 24, 25, 26, 27]. These results
emphasise further its significance as an example of an elegant dynamical
system with a repertoire of rich properties, among which, hyperchaoticity is
prominent.

The paper is organised as follows: In Sec. 2, we review the concept of
circuits introduced by Rössler and Thomas and their role in the appearance of
chaos and hyperchaos in dynamical systems. In Sec. 3, we explore the dynamics
in the neighbourhood of the fixed points of labyrinth walks and discuss how
there can exist chaos and hyperchaos in the absence of attractors.
Subsequently, in Sec. 4, we turn our attention to a more complex setting;
namely, to that of an array of coupled TR systems and the role of labyrinth
walks in the emergence of chimera-like states. In Sec. 5, we discuss further
implications in relation to the labyrinth walks system by showing that even
though it is volume-preserving, it is not force-conservative. Finally, we
conclude our work in Sec. 6, where we discuss possible future directions of
work.

##  2 Feedback Circuits and Chaos

The concepts of circuits and feedback stem from electromechanical systems and
biology. However, they were redefined by Rössler and Thomas to describe
certain properties of dynamical systems and their equations of motion. We own
them that “complex dynamics can be thought of in terms of feedback circuits”
[7]. According to them, a feedback circuit (in short, a circuit) is the
influence of a variable to its own evolution. This can happen directly or
indirectly by the coupling of the variable with other variables in the system.
Importantly, the coupling can be linear or nolinear and their realisation was
pivotal in connecting graph theory with multistationarity and chaos [28, 29].

In particular, what proved fruitful was that the feedback circuits are
directly identifiable through the Jacobian matrix $J$ of a dynamical system,

$\frac{dX}{dt}=F(X),$

defined by

$J\equiv J(X)=\left[\begin{array}[]{cccc}J_{11}&J_{12}&\cdots&J_{1n}\\\
J_{21}&J_{22}&\cdots&J_{2n}\\\ \vdots&\vdots&\ddots&\vdots\\\
J_{n1}&J_{n2}&\cdots&J_{nn}\end{array}\right],$

where $t$ is the time, $X=(X_{1},\dots,X_{n})\in\mathbb{R}^{n}$ the solution-
vector and $F\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ a differentiable function
with $F(X)=(F_{1}(X),\dots,F_{n}(X))$. The entries in the Jacobian matrix are
the partial derivatives

$J_{ij}\equiv J_{ij}(X)=\frac{\partial F_{i}}{\partial X_{j}},$

where $i,j=1,\ldots,n$.

Consequently, if $J_{ij}=0$, variable $X_{j}$ does not influence variable
$X_{i}$. However, if $J_{ij}$ is different than zero, then variable $X_{j}$
influences variable $X_{i}$. This is in the sense that if
$J_{ij}=\frac{\partial F_{i}}{\partial X_{j}}$ is positive, then there is a
positive feedback from variable $X_{j}$ to variable $X_{i}$, whereas if
$J_{ij}=\frac{\partial F_{i}}{\partial X_{j}}$ is negative, it results in a
negative feedback from variable $X_{j}$ to variable $X_{i}$. Of course, the
values of the entries in $J$ might depend on the solution $X(t)$.

Since a non-zero entry in $J$ implies that variable $X_{j}$ influences
variable $X_{i}$, one can define a vertex from $X_{j}$ to $X_{i}$ in a graph
of incidences (see Figs. 1 and 2 for two examples). Thus, circuits are the
non-zero entries in $J$, with their indices forming cyclic permutations with
each other.

Figure 1: Graph of incidences of the circuits in the Jacobian matrix $J$ of
the system of Eq. (1), where $n=3$ and $b\neq 0$. The curved arrows indicate
self-loops (one-element circuits), where entries in $J$ influence themselves
as they are its diagonal elements.

For example, for the labyrinth chaos system of Eq. (1) where $n=3$, the graph
of incidences is shown in Fig. 1 for $b\neq 0$, which shows its three element
circuits, the closed path $X_{1}\rightarrow X_{3}\rightarrow X_{2}\rightarrow
X_{1}$ and the three, one-element, circuits which are self-loops, due to the
diagonal elements in $J$.

In contrast, in Fig. 2, the graph of incidences of the labyrinth walks system
(2), where $n=3$ and $b=0$, has only one, three-element, circuit, which is the
same as the three-element circuit in system (1), shown in Fig. 1. However,
there is no one-element circuit, since the diagonal elements in $J$ are zero,
which implies there are no self-loops.

Figure 2: Graph of incidences of the circuits in the Jacobian matrix $J$ of
the system of Eq. (2), where $n=3$. Note the absence of self-loops as the
diagonal elements in $J$ are equal to 0, thus there are no one-element
circuits, in contrast to the graph of incidences of system (1) in Fig. 1.

The circuits that can be identified in the Jacobian matrix are directly
implicated in its characteristic equation, eigenvalues, eigenstates and
determinant. Importantly, they are also implicated in the stability of fixed
points and chaoticity, quantified by the spectrum of Lyapunov exponents.

Furthermore, circuits can be positive or negative depending on the sign of the
product of their terms, i.e. on the signs of the entries in $J$. In this
sense, the parity of the number of negative interactions of a circuit
determines its sign: if it is even, the circuit is positive and if it is odd,
it is negative [30]. If its sign depends on the values of the solution $X$,
i.e. on the position in the state-space, then the circuit is called an
ambiguous circuit.

The locus of points that the sign of a non-ambiguous circuit changes from
positive to negative is called a frontier. Depending of the character of
change of the eigenvalues of $J$, the frontiers have been categorised in [30],
where computer codes are also offered for their calculation.

The principal frontier [30], F1, is the locus of points where the sign of the
determinant of the Jacobian matrix, $\det J$, changes and partitions the phase
space according to the sign of the product of its eigenvalues. This is so as
$\prod\lambda_{k}=\det J$, where $\lambda_{k}$ are the eigenvalues of $J$.

Consequently, positive and negative circuits play contrasting roles in the
dynamics of a system. The necessary condition for multistationarity is the
presence of a positive circuit, whereas the presence of a negative circuit is
a necessary condition for the existence of a stable attractor. Coexisting
positive and negative circuits are necessary conditions for oscillations,
chaos and hyperchaos (see conditions (i)-(iv)). Based on their work on TR
systems and labyrinth walks, Rössler, Thomas and co-workers, proposed and
demonstrated a necessary condition for the existence of hyperchaos: namely
that hyperchaos of order $m$ (i.e. with $m$ positive Lyapunov exponents) can
be generated by a single feedback circuit in $n=2m+1$ dimensions [7].

The idea of feedback circuits is a common theme in engineering, biology,
network science, nonlinear dynamics and complexity. The formulation of
feedback circuits for the analysis of dynamical systems, as discussed by
Rössler, Thomas, Rosen, Shilnikov, Nicolis and other pioneers in nonlinear
dynamics and chaos theory, has become an integral part of the vocabulary, not
only in relation to biological systems [10] per se, but also, in complexity
theory at large [11, 31, 32].

We shall now turn to the role of feedback circuits in the dynamics of the
surprisingly elegant 3D labyrinth walks system.

##  3 Chaos and Hyperchaos in the Absence of Attractors

We focus now on the system [7]

$\frac{dX_{i}}{dt}=\sin(X_{i+1})-bX_{i}$ (1)

where $i=1,\dots,n$, $t$ is the time and $b$ a real constant. In this context,
$X_{n+1}=X_{1}$, signifying the cyclic permutation in variables $X_{i}$. The
system has $n$ one-element circuits, which makes it ideal for checking
graphically, analytically and computationally the ideas developed by Rössler
and Thomas. The particular case $b=0$ corresponds to the labyrinth walks
system that we study next.

In particular, we study the case $n=3$ and $b=0$, what amounts to the 3D
labyrinth walks system

$\displaystyle\frac{dX_{1}}{dt}$ $\displaystyle=\sin(X_{2})=F_{1},$ (2)
$\displaystyle\frac{dX_{2}}{dt}$ $\displaystyle=\sin(X_{3})=F_{2},$
$\displaystyle\frac{dX_{3}}{dt}$ $\displaystyle=\sin(X_{1})=F_{3},$

where $F=(F_{1},F_{2},F_{3})$. System (2) can be written in the compact form

$\dot{X}=P_{M}f(X),$

where $X=(X_{1},X_{2},X_{3})$,
$\dot{X}=(\dot{X}_{1},\dot{X}_{2},\dot{X}_{3})=\Bigl{(}\frac{dX_{1}}{dt},\frac{dX_{2}}{dt},\frac{dX_{3}}{dt}\Bigr{)}$,

$f(X)=\left[\begin{array}[]{c}\sin\left(X_{1}\right)\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\cr\sin\left(X_{2}\right)\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\sin\left(X_{3}\right)\end{array}\right],$

and

$P_{M}=\left[\begin{array}[]{ccc}0&1&0\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr 0&0&1\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
1&0&0\end{array}\right]$

is the permutation matrix.

The Jacobian matrix of system (2) is given by

$\displaystyle J$
$\displaystyle=\left[\begin{array}[]{ccc}0&\cos\left(X_{2}\right)&0\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\cr 0&0&\cos\left(X_{3}\right)\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr\cos\left(X_{1}\right)&0&0\end{array}\right]$ (3)
$\displaystyle=P_{M}\left[\begin{array}[]{ccc}\dot{F_{3}}&0&0\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr 0&\dot{F_{1}}&0\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr 0&0&\dot{F_{2}}\end{array}\right].$

Thus, the signs of circuit-changes occur in an alternating, nonlinear fashion
depending on the cosines in $\dot{F}=(\dot{F_{1}},\dot{F_{2}},\dot{F_{3}})$,
due to the cyclic permutations encoded in $P_{M}$. The determinant of $J$ is
given by

$\det J(X)=\cos(X_{1})\cos(X_{2})\cos(X_{3}).$

In the general case of $n$ variables, one can easily see that $J$ has the form
of the permutation matrix $P_{M}$ in $n$ dimensions, where its entries are
given by the generalisation of Eq. (3) in $n$ dimensions. It follows that in
the $n$D case, its determinant is given by the product of the cosines of all
variables, i.e. by

$\det J(X)=\prod_{i=1}^{n}\cos(X_{i}),$

since the determinant of $P_{M}$ is 1 and the determinant of the corresponding
diagonal matrix (see Eq. (3) for $n=3$) is
$\prod_{i=1}^{n}\dot{F}_{i}=\prod_{i=1}^{n}\cos(X_{i})$. This shows that the
principal frontier F1 of system (2), where the circuits change sign [30], has
zero dimension. This is because F1 is an $nD$-lattice of points defined by the
solution to the equation $\prod_{i=1}^{n}\cos(X_{i})=0$. This is a set of
distinct points and not a surface or hyper-surface, what is common in other
systems [30, 12] studied by Rössler, Thomas and co-workers.

(a)

(b)

Figure 3: Unit-cell of fixed points in panel (a) and part of the 3D fixed-
points lattice with a labyrinth walk trajectory in panel (b). (a) The black
and grey sets of unstable fixed points $X_{i}=k\pi,\;i=1,2,3$ of the labyrinth
walks system (2) in the unit cell for $k=0$ and 1. We note that the unit-
length of the cell is $\pi$ in all 3 directions and that the state-space of
system (2) in panel (b) is composed of collated versions of the cell in (a) in
all 3 directions. (b) Part of the lattice of fixed points in the 3D state-
space of system (2) (black points) and a labyrinth walks trajectory (grey
curve).

The stability analysis of system (2) shows there is a simple, yet elegant
structure endowed with an intricate symmetry. The fixed points are the
solutions to the equation $\frac{dX}{dt}=0$, therefore solutions to the system
of equations

$\displaystyle\sin(X_{1})$ $\displaystyle=0,$ $\displaystyle\sin(X_{2})$
$\displaystyle=0,$ $\displaystyle\sin(X_{3})$ $\displaystyle=0.$

This leads to a countably-infinite set of fixed points, arranged in the 3D
lattice where $\sin(X_{i})=0,\;i=1,2,3$, i.e. the fixed points are given by
$X_{i}=k\pi,\;k\in\mathbb{N}$ (see Fig. 3(b) for a subset of them). We note
that these properties define an infinitely-large lattice of fixed points in
the 3D state-space of the system composed of collated, identical cubic cells
of unit-length $\pi$ (see Fig. 3(a) for an example of such a cubic cell and
(b) for a part of the infinitely-large 3D lattice of fixed points, what
amounts to a number of collated, cubic cells of unit-length $\pi$). These
properties can be carried out straightforwardly for any $n$. For example, for
$k=0$ and 1, the set of fixed points of system (2) is given by $\\{(0,0,0)$,
$(0,0,\pi)$, $(0,\pi,0)$, $(\pi,0,0)$, $(\pi,\pi,0)$, $(\pi,0,\pi)$,
$(0,\pi,\pi)$, $(\pi,\pi,\pi)\\}$ and is shown in Fig. 3(a). These are the
only fixed points of the system and they all turn out to be unstable, what
amounts to the absence of attractors. In particular, the eigenvalues of $J$
for each of the eight fixed points evaluate either to the cubic roots of
unity, $\omega_{1,2,3}=\Big{\\{}1,\frac{-1\pm i\sqrt{3}}{2}\Big{\\}}$, or to
their opposites, $\omega_{1,2,3}=\Big{\\{}-1,\frac{1\mp
i\sqrt{3}}{2}\Big{\\}}$. Thus, the fixed points belong to two sets: The first
set (black points in Fig. 3(a)) consists of fixed points with one positive
real eigenvalue and two complex, with negative real parts, eigenvalues, what
signifies an outgoing unstable, spiral motion. The second set (grey points in
Fig. 3(a)) consists of fixed points with one negative real eigenvalue and two
complex, with positive real parts, eigenvalues, what signifies an incoming,
unstable, spiral motion. All fixed points, in both sets, have complex
eigenvectors attributed to their complex eigenvalues, which adds a very
interesting aspect for further investigations in terms of spinor properties
[33]. All fixed points are thus unstable and give rise to a structure that can
sustain complex unstable periodic orbits, wandering, space-filling
trajectories (the so-called, labyrinth walks) and thus, to chaotic behaviour.

(a)

(b)

Figure 4: Time evolution of $X_{1}$, $X_{2}$ and $X_{3}$ of the grey curve in
Fig. 3(b) and its projections in the $(X_{3},X_{1})$, $(X_{2},X_{1})$ and
$(X_{3},X_{2})$ planes. (a) Evolution of $X_{1}$, $X_{2}$ and $X_{3}$ of the
grey curve shown in Fig. 3(b) over time, what constitutes an example of a
labyrinth walks trajectory. (b) Projections of the same trajectory in the
$(X_{3},X_{1})$, $(X_{2},X_{1})$ and $(X_{3},X_{2})$ planes.

In Fig. 4(a), we present the time-evolution of $X_{1}$, $X_{2}$ and $X_{3}$ of
system (2) for the labyrinth walks, grey trajectory shown in Fig. 3(b) and in
Fig. 4(b), its projection in the planes $(X_{3},X_{1})$, $(X_{2},X_{1})$ and
$(X_{3},X_{2})$, where some very interesting patterns can be seen. This
trajectory is an example of a wandering, space-filling trajectory in the 3D
state-space of system (2), part of which is shown in panel (b) in Fig. 3, what
amounts to the fractional, super-diffusive character [13, 34] of labyrinth
walks.

The results for the 3D case can be readily generalised in $n$ dimensions,
where the unit cells will be $n$D hypercubes, building up an $n$D lattice of
unstable fixed points, located at their vertices. In this case, the
eigenvalues of the $n\times n$ Jacobian matrix $J$ will be the $n^{th}$ roots
of unity, resulting in a structure of higher-order symmetries that can sustain
higher-order hyperchaos without the presence of attractors. In [7], it has
been reported that systems with even and odd dimensionality have been examined
and it has been predicted by Rössler and Thomas that they exhibit hyperchaos
of order $m$, i.e. $m$ positive Lyapunov exponents, for $n=2m+1$. One-element
circuits have been verified for $n=5,7,25,99$.

It is the existence of infinitely-many, unstable fixed points arranged on an
infinite lattice, where there are no attractors that makes the labyrinth walks
system so interesting. Even though it is deterministic, these properties lead
to fractional-like chaotic motion reported in [13, 34], reminiscent of Levy-
flights in stochastic systems! The coexistence of channels of free, ballistic
flights along with unstable periodic orbits brings in mind trajectories and
motion in infinite-horizon Lorenz-gas type billiards [35]. This combination of
properties makes labyrinth walks not only an elegant system to study but also,
a good example to elucidate the concept of volume preservation, a topic we
will touch upon in Sec. 5.

Before that, we turn our attention in the next section to coupled TR systems
(1) and revisit their surprising role in sustaining chimera-like states, a
peculiar synchronisation phenomenon.

##  4 Chimera-like States in arrays of Coupled Labyrinth Walks and TR systems

The fractional, super-diffusive character [13, 34] of labyrinth walks provides
the mechanism for the emergence of chimera-like states in arrays of coupled 3D
TR systems [9] (see Eqs. (1)). In this section, we revisit the recent results
in [9], in light of the ability of labyrinth walks and its lattice of unstable
fixed points to facilitate complex, chimera-like states, even with very strong
coupling.

To this end, the authors in [9] considered the following system of coupled 3D
labyrinth walks and TR systems arranged, in a simple ring topology with non-
local coupling

$\displaystyle\frac{dX_{1}^{k}}{dt}$
$\displaystyle=-b_{k}X_{1}^{k}+\sin\big{(}X_{2}^{k}\big{)}+\frac{d}{2P}\sum_{j=k-P}^{k+P}\big{(}X_{1}^{k}-X_{1}^{j}\big{)},$
(4) $\displaystyle\frac{dX_{2}^{k}}{dt}$
$\displaystyle=-b_{k}X_{2}^{k}+\sin\big{(}X_{3}^{k}\big{)},$
$\displaystyle\frac{dX_{3}^{k}}{dt}$
$\displaystyle=-b_{k}X_{3}^{k}+\sin\big{(}X_{1}^{k}\big{)},$

where $k=1,\dots,n$, $2P\leq n$ is the number of coupled neighbours, $d\geq 0$
is the strength of the linear coupling and $b_{k}\geq 0$ for all $k$. The term
$\frac{1}{2P}$ is a normalisation constant that distributes equally the linear
coupling among the $n$ coupled systems.

(a)

(b)

Figure 5: Snapshots of temporal phenomena of coherent and incoherent patterns,
reminiscent of chimera-like states in system (4) (see text for the details) as
depicted in their spatio-temporal evolution in panel (b) in Fig. 6. (a)
Snapshots at time $t\approx 13401$ and (b) at a later time $t\approx 13404$.
Note that here, $d=0.6$, $n=40$, $b_{k}=0$ for $k=1,\ldots,20$ (labyrinth
walks) and $b_{k}=0.17$ for $k=21,\ldots,40$ (labyrinth chaos). Note also that
$X_{1}^{k}$ values seen at almost the same height (values at the vertical
axes), imply they are locked (see text for more details). The group of locked
$X_{1}^{k}$ values is the coherent group and the group of non-locked,
$X_{1}^{k}$ values, the incoherent group.

(a)

(b)

Figure 6: Hyperchaotic and, spatio-temporal behaviour and emergence of
alternating patterns of chimera-like states in system (4). (a) Evolution of
the first four largest Lyapunov exponents $\lambda_{i}$, $i=1,\ldots,4$, seen
to converge to positive, non-zero values (hyperchaoticity) and (b) chimera-
like behaviour in $n=40$ TR and labyrinth walks systems (4) (see text for the
details). Note that here $d=0.6$, $b_{k}=0$ for $k=1,\ldots,20$ (labyrinth
walks), $b_{k}=0.17$ for $k=21,\ldots,40$ (labyrinth chaos) and that the whole
system is hyperchaotic with more than one positive Lyapunov exponents, as
shown in panel (a).

The system of Eqs. (4) was introduced and studied in [9], where $n=40$ TR
systems were considered, with the first half ($k=1,\ldots,20$) performing
labyrinth walks ($b_{k}=0$) and the other half set in the labyrinth chaos
regime with $b_{k}=0.18$ for $k=21,\ldots,40$. Each system was coupled with
its $P=20$ nearest-neighbour systems in both sides (10 + 10), with the
coupling strength set at $d=0.6$. For these parameters and setting, system (4)
is hyperchaotic, with more than 1 positive Lyapunov exponents.

For such a strong coupling ($d=0.6$) and without the presence of labyrinth
walks in the second part of the system, one would expect the prevalence of
chaotic synchronisation, what would result in an almost, fully synchronised
system over time; yet, due to the presence of labyrinth walks, the emergence
of partial, spatio-temporal synchronisation, reminiscent of chimera-like
states [20, 21, 22, 24, 23, 26, 27] was observed [9]. In this context, partial
is meant in the sense that the absolute differences of the $x_{1}^{k}$,
$k=1,\ldots,n$ solutions were less than a threshold of the order of $10^{-4}$
to $10^{-3}$, i.e. they were locked in time, but not continuously over time.
Here, we show in Figs. 5 and 6, the results of a similar study where the only
difference with the study in [9] is that $b_{k}=0.17$ for $k=21,\ldots,40$,
which corresponds again to labyrinth chaos. In particular, Fig. 5(a) shows an
example of locking at two different times as depicted in the two panels for
$t\approx 13401$ and $t\approx 13404$. Even though, these differences are
locked (but not continuously) over time, the system remains hyperchaotic as a
whole, as shown in Fig. 6(a), where the first four largest Lyapunov exponents
$\lambda_{i}$, $i=1,\ldots,4$, can be seen to converge to positive, non-zero,
values.

Apparently, it is the presence of chaotic walks that gives rise to the
emergence of chimera-like states, what offers a novel phenomenon of chaotic
synchronisation that is worth analysing further, both in terms of labyrinth
walks and in terms of the appearance of chimera-like states in arrays of
hyperchaotic dynamics. Moreover, these results are promising as they also
offer the opportunity to study the emergence of chimera-like states in scale-
free, small-world and Erdős-Rényi networks of coupled TR systems.

##  5 Labyrinth Walks: A Volume Preserving, non Force-conservative System

We now turn to two very important properties of the labyrinth walks system
(2).

The first one is that it is volume-preserving [36, 37], in the sense it is
preserving volumes in its 3D state-space. Indeed, if we write system (2) in
the form

$\dot{X}=F(X),$

where, we recall, $X=(X_{1},X_{2},X_{3})$ and
$F=(F_{1},F_{2},F_{3})=(\sin(X_{2}),\sin(X_{3}),\sin(X_{1}))$, the
incompressibility condition (divergence) gives

$\nabla\dot{X}=\nabla F=\frac{\partial F_{1}}{\partial X_{1}}+\frac{\partial
F_{2}}{\partial X_{2}}+\frac{\partial F_{3}}{\partial X_{3}}=0.$ (5)

This means the flow $F$ of system (2) is divergence-free, i.e. it is
incompressible. This, in turn, implies that the solution to Liouville’s
equation [38]

$\frac{dV}{dt}=\int_{V}\big{(}\nabla\dot{X}\big{)}dV$

is

$V(t)=e^{(\nabla\dot{X})t}V(0)=V(0),$

where $V$ is the volume of a set of initial conditions in the 3D state-space
of system (2). The last result shows that volumes in its state-space are
preserved in time, and thus it is a volume-preserving dynamical system. We
note here that Eq. (5) is essentially the sum of the diagonal entries in the
Jacobian matrix $J$ (see Eq. (3)) of system (2), thus it is equal to its
trace. In other words, the divergence of the flow is the trace of the Jacobian
matrix of the dynamical system.

As we shall show next, the second property is that system (2) is not force-
conservative, in the sense that

$\oint\overrightarrow{F}.\hskip 2.84526pt\overrightarrow{dr}\neq 0.$ (6)

The “acceleration” field $\ddot{X}$, that would be proportional to a force
$\mathcal{F}$ (from Newton’s second law), can be obtained by computing the
time-derivative of Eq. (2). This yields

$\displaystyle\frac{d^{2}X_{1}}{dt^{2}}$
$\displaystyle=\sin(X_{3})\cos(X_{2}),$
$\displaystyle\frac{d^{2}X_{2}}{dt^{2}}$
$\displaystyle=\sin(X_{1})\cos(X_{3}),$
$\displaystyle\frac{d^{2}X_{3}}{dt^{2}}$
$\displaystyle=\sin(X_{2})\cos(X_{1}),$

which can be written in the compact form

$\ddot{X}=\mathcal{F},$

with $\ddot{X}=(\ddot{X_{1}},\ddot{X_{2}},\ddot{X_{3}})$ and

$\mathcal{F}=(\mathcal{F}_{1},\mathcal{F}_{2},\mathcal{F}_{3}),$

where

$\displaystyle\mathcal{F}_{1}$ $\displaystyle=\sin(X_{3})\cos(X_{2}),$ (7)
$\displaystyle\mathcal{F}_{2}$ $\displaystyle=\sin(X_{1})\cos(X_{3}),$
$\displaystyle\mathcal{F}_{3}$ $\displaystyle=\sin(X_{2})\cos(X_{1}).$

The potential $U$, if it existed, would be related to the force $\mathcal{F}$
by

$\mathcal{F}=-\nabla{U},$

it would be path independent and would represent the “mechanical energy” of
the system, when added to its “kinetic energy”. In this context, $\nabla$ is
the gradient of the scalar-valued differentiable function
$U:\mathbb{R}^{3}\rightarrow\mathbb{R}$ that we want to check whether it is
the potential energy of system (2).

As we show here, this is not the case as we can find a counter-example of two
different paths (solid and dashed in Fig. 7) along which $U$ is not the same
function: To this end, let us integrate $\mathcal{F}$ along the two paths
shown in Fig. 7. The solid-line path starts at the origin $O(0,0,0)$ and ends
at $M(X_{1},X_{2},X_{3})$, following the three straight lines from $O(0,0,0)$
to $P(X_{1},0,0)$ to $Q(X_{1},X_{2},0)$ to $M(X_{1},X_{2},X_{3})$ and the
dashed-line path, following the straight lines from $O(0,0,0)$ to
$H(0,0,X_{3})$ to $M(X_{1},X_{2},X_{3})$.

Figure 7: An example of two paths (depicted by solid and dashed lines) along
which $U$ is different, implying that $\mathcal{F}$ is not a conservative
force and thus, the labyrinth walks system (2) is not force-conservative.

Following the dashed-line path, $U$ reads

$\displaystyle U=$
$\displaystyle-\int_{0}^{X_{1}}\mathcal{F}_{1}(X^{\prime}_{1},X^{\prime}_{2},X^{\prime}_{3})\Bigg{|}_{\begin{subarray}{c}X_{2}^{\prime}{}=0\\\
X_{3}^{\prime}{}=0\end{subarray}}dX_{1}^{\prime}{}$
$\displaystyle-\int_{0}^{X_{2}}\mathcal{F}_{2}(X^{\prime}_{1},X^{\prime}_{2},X^{\prime}_{3})\Bigg{|}_{\begin{subarray}{c}X_{1}^{\prime}{}=X_{1}\\\
X_{3}^{\prime}{}=0\end{subarray}}dX_{2}^{\prime}{}$
$\displaystyle-\int_{0}^{X_{3}}\mathcal{F}_{3}(X^{\prime}_{1},X^{\prime}_{2},X^{\prime}_{3})\Bigg{|}_{\begin{subarray}{c}X_{1}^{\prime}{}=X_{1}\\\
X_{2}^{\prime}{}=X_{2}\end{subarray}}dX_{3}^{\prime}{}.$

Substituting $\mathcal{F}_{1}$, $\mathcal{F}_{2}$ and $\mathcal{F}_{3}$ from
Eq. (7) into the last equation, $U$ yields

$U=-X_{2}\sin(X_{1})-X_{3}\sin(X_{2})\cos(X_{1}).$ (8)

Now, following the dashed-line path, $U$ reads

$\displaystyle U=$
$\displaystyle-\int_{0}^{X_{3}}\mathcal{F}_{1}(X^{\prime}_{1},X^{\prime}_{2},X^{\prime}_{3})\Bigg{|}_{\begin{subarray}{c}X_{1}^{\prime}{}=0\\\
X_{2}^{\prime}{}=0\end{subarray}}dX_{3}^{\prime}{}$
$\displaystyle-\int_{H}^{M}\mathcal{F}\hskip
2.84526pt(dX_{1}^{\prime}{}+dX_{2}^{\prime}{}),$

which can be written as

$\displaystyle U=$
$\displaystyle-\int_{0}^{X_{3}}\mathcal{F}_{3}(X^{\prime}_{1},X^{\prime}_{2},X^{\prime}_{3})\Bigg{|}_{\begin{subarray}{c}X_{1}^{\prime}{}=0\\\
X_{2}^{\prime}{}=0\end{subarray}}dX_{3}^{\prime}{}$
$\displaystyle-\int_{0}^{X_{1}}\mathcal{F}_{1}(X^{\prime}_{1},X^{\prime}_{2},X^{\prime}_{3})\Bigg{|}_{\begin{subarray}{c}X_{3}^{\prime}{}=X_{3}\\\
X_{2}^{\prime}{}=(\frac{X_{2}}{X_{1}})X_{1}^{\prime}{}\end{subarray}}dX_{1}^{\prime}{}$
$\displaystyle-\int_{0}^{X_{2}}\mathcal{F}_{2}(X^{\prime}_{1},X^{\prime}_{2},X^{\prime}_{3})\Bigg{|}_{\begin{subarray}{c}X_{3}^{\prime}{}=X_{3}\\\
X_{1}^{\prime}{}=(\frac{X_{1}}{X_{2}})X_{2}^{\prime}{}\end{subarray}}dX_{2}^{\prime}{}.$

Again, substituting Eq. (7) into the last equation, $U$ becomes

$U=-\frac{{X_{1}}}{X_{2}}\sin(X_{2})\sin(X_{3})+\frac{{X_{2}}}{X_{1}}\cos(X_{1})\cos(X_{3}).$
(9)

Comparing, $U$ from Eqs. (8) and (9) computed along the two paths, we deduce
they are not the same function. Therefore, $U$ depends on the path, a
potential function cannot be defined and system (2) is not force-conservative,
since Eq. (6) holds.

This discussion is motivated further by the fact that when one studies the
volume-preservation properties of system (2), given that the trace of its
Jacobian matrix is zero (what corresponds to a divergence-free, incompressible
flow) and in view of Liouville’s theorem for Hamiltonian systems [36, 37, 39],
one is tempted to consider whether it falls in the realm of Kolmogorov-Arnold-
Moser (KAM) theory [40]. This might be driven further by the computation of
Poincaré surfaces of section, such as those in [12, 13], where islands of
stability can be seen within the chaotic sea and where the diffusive aspect of
wandering trajectories is similar in appearance to those of chaotic
Hamiltonian systems [41, 42, 43]. Yet, this is where similarities stop as
labyrinth walks (2), even though is not a force-conservative system, preserves
volumes in its 3D state-space. However, the fact it preserves volumes does not
mean there is an integral of motion. The existence of a global conserved
quantity in the 3D state-space is prohibited since that would render the
system effectively two dimensional and hence chaos could not appear.

##  6 Conclusions & Outlook

As all Rössler’s pioneering contributions, labyrinth chaos still holds promise
for very interesting further developments. Its simplicity and elegance, both
in terms of symmetries, topology and feedback-circuit structure, makes it a
good candidate to compare it with other nonlinear, cyclically coupled systems,
such as the Arabesques [8, 12], the Lotka-Voltera system and its variants
[44], and the Arnold-Beltrami-Childress (1:1:1 ABC) model [45, 46].

In connection to arrays of coupled TR systems, where chimera-like states have
been shown to emerge even under strong coupling settings, the investigation of
the role of inhomogeneities in parameters and in network topology (e.g. scale-
free, small-world, Erdős-Rényi) are worth studying further as well as in
connection to other systems [27, 47].

Here, we discussed further implications pertaining to the labyrinth walks
system by showing that, even though it is volume-preserving, it is not force-
conservative. This makes it a fine and elegantly simple example of a low-
dimensional volume-preserving, chaotic system with no attractors and with only
unstable fixed points. We note that in Statistical Mechanics, where phase
spaces are high-dimensional, there are also examples of systems that are
volume-preserving, chaotic, with no attractors and with only unstable fixed
points, such as extensions of the Nose-Hoover thermostat known as the Martyna-
Klein-Tuckerman, chain thermostats [48].

The investigation of the appropriate Hamiltonian (non-holonomic) or Lagrangian
(with constrains) formalism that might encompass the labyrinth walks system,
could reveal deeper connections between chaos, volume preservation, symplectic
structure, diffusion and KAM theory. It would be interesting to seek a
comparison of the labyrinth walks system with the class of systems pertaining
to the construction of Hoover thermostats in Statistical Mechanics [49, 50,
48]. This is because this class of systems provides: (a) a framework for the
required analysis by determining phase-space distributions assuming ergodicity
and (b) expresses the balance between the kinetic, dissipative and potential
parts, what makes their phase-space incompressible.

Finally, information transfer in networks of coupled labyrinth walks systems
as “carriers of information” [51, 52, 53] emerges as another exciting field of
future research, since it is related to the volume-preservation property of
individual labyrinth walks systems.

##  7 Acknowledgements

This research was supported by our respective Institutes: Department of
Physics, Qom University of Technology for AL, Service de Physique des Systèmes
Complexes et Mécanique Statistique and CeNoLi, Université Libre de Bruxelles
for VB and Department of Mathematical Sciences, University of Essex for CGA.

## Availability of Data

The data that support the findings of this study are available within the
article.

## References

  * [1] R. Gilmore and M. Lefranc.  The Topology of Chaos.  Wiley, 2002.
  * [2] R. Gilmore, C. Letellier, and M. Lefranc.  Chaos topology.  Scholarpedia, 3(7):4592, 2008.
  * [3] J. C. Sprott.  Elegant Chaos: Algebraically Simple Chaotic Flows.  World Scientific, 2010.
  * [4] P. Gaspard and G. Nicolis.  What can we learn from homoclinic orbits in chaotic dynamics?  Journal of Statistical Physics, 31:499–518, 1983.
  * [5] C. Letellier and O. Rössler.  Rössler attractor.  Scholarpedia, 1(10):1721, 2006.
  * [6] R. Thomas.  Deterministic chaos seen in terms of feedback circuits: Analysis, synthesis, ‘Labyrinth Chaos’.  Int. J. Bifurcation and Chaos, 9(10):1889–1905, 1999.
  * [7] R. Thomas, V. Basios, M. Eiswirth, T. Kruel, and O. Rössler.  Hyperchaos of arbitrary order generated by a single feedback circuit, and the emergence of chaotic walks.  Chaos, 14(3):669–674, 2004.
  * [8] R. Thomas and O. Rössler.  Genèse de Formes.  Revue des Questions Scientifiques, 177:271–287, 2006.
  * [9] V. Basios and C. G. Antonopoulos.  Hyperchaos and Labyrinth chaos: Revisiting Thomas - Rössler systems.  Journal of Theoretical Biology, 460:153–159, 2019.
  * [10] D. Thieffry and M. Kaufman.  Prologue to the special issue of JTB dedicated to the memory of René Thomas (1928 - 2017): A journey through biological circuits, logical puzzles and complex dynamics.  Journal of Theoretical Biology, 474:42–47, 2019.
  * [11] G. Nicolis and C. Nicolis.  Foundations of Complex Systems - Nonlinear Dynamics, Statistical Physics, Information and Prediction.  World Scientific, 2010.
  * [12] C. G. Antonopoulos, V. Basios, J. Demongeot, P. Nardone, and R. Thomas.  Linear and nonlinear arabesques: A study of closed chains of negative 2-element circuits.  Int. J. Bifurcation and Chaos, 23:1330033, 2013.
  * [13] J. C. Sprott and K. E. Chlouverakis.  Labyrinth chaos.  Int. J. Bifurcation and Chaos, 17(6):2097–2108, 2007.
  * [14] O. Rössler.  An equation for hyperchaos.  Physics Letters A, 71, 155-157, 1979.
  * [15] L. Chunguang and C. Guanrong.  Chaos and hyperchaos in the fractional-order Rössler equations.  Physica A: Statistical Mechanics and its Applications, 341:55–61, 2004.
  * [16] A. V. Andreev, A. G. Balanov, T. M. Fromhold, M. T. Greenaway, A. E. Hramov, W. Li, V. V. Makarov, and A. M. Zagoskin.  Chaos and hyperchaos in driven interacting quantum systems, 2019.
  * [17] O. Rössler and C. Letellier.  Chaos: The World of Nonperiodic Oscillations.  Understanding complex systems series, Springer, 2020.
  * [18] C. Letellier and O. Rössler.  Hyperchaos.  Scholarpedia, 2(8):1936, 2007.
  * [19] Y. Kuramoto and D. Battogtokh.  Coexistence of coherence and incoherence in nonlocally coupled phase oscillators.  Nonlinear Phenomena Complex Systems, 5(4):380 – 385, 2002.
  * [20] D. M. Abrams and S. H. Strogatz.  Chimera states for coupled oscillators.  Phys. Rev. Lett., 93(17):174102, 2004.
  * [21] M. J. Panaggio and D. M. Abrams.  Chimera states: Coexistence of coherence and incoherence in networks of coupled oscillators.  Nonlinearity, 28(3):R67–R87, 2015.
  * [22] J. Hizanidis, N. Kouvaris, G. Zamora-López, A. DÍan Guilera, and C. G. Antonopoulos.  Chimera-like states in modular neural networks.  Sci. Rep., 6, 2016.
  * [23] F. Parastesh, S. Jafari, H. Azarnoush, B. Hatef, and A. Bountis.  Imperfect chimeras in a ring of four-dimensional simplified Lorenz systems.  Chaos, Solitons and Fractals, 110, 2018.
  * [24] A. Schmidt, T. Kasimatis, J. Hizanidis, A. Provata, and P. Hövel.  Chimera patterns in two-dimensional networks of coupled neurons.  Phys. Rev. E, 95:032224, 2017.
  * [25] A Schmidt, T Kasimatis, J Hizanidis, A Provata, and P Hövel.  Erratum –chimera patterns in two-dimensional networks of coupled neurons.  Phys. Rev. E, 100:069901, 2019.
  * [26] E. Sigalas, A. Bezerianos J. Li, and C. G. Antonopoulos.  Emergence of chimera-like states in prefrontal-cortex macaque intracranial recordings.  In International Workshop on Pattern Recognition in Neuroimaging (PRNI), Singapore, 2018.
  * [27] K. Kaneko.  From globally coupled maps to complex-systems biology.  Chaos, 25(097608), 2015.
  * [28] C. Soulé.  Graphic requirements for multistationarity.  ComPlexUs, 1:123–133, 2003.
  * [29] M. Kaufman, C. Soulé, and R. Thomas.  A new necessary condition on interaction graphs for multistationarity.  Journal of Theoretical Biology, 248:675–685, 2007.
  * [30] R. Thomas and P. Nardone.  A further understanding of phase space partition diagrams.  International Journal of Bifurcation and Chaos, 19(03):785–804, 2009.
  * [31] G. Nicolis and V. Basios.  Chaos, Information Processing and Paradoxical Games.  World Scientific, 2015.
  * [32] M. Kaufman and R. Thomas.  Emergence of complex behavior from simple circuit structures.  C. R. Acad. Sci. Paris Biologies, 326:205–214, 2003.
  * [33] E. M. S. Hitzer.  Imaginary eigenvalues and complex eigenvectors explained by real geometry. In: Dorst L., Doran C., Lasenby J. (eds) Applications of Geometric Algebra in Computer Science and Engineering.  Birkhäuser, Boston, MA, 2002.
  * [34] K. E. Chlouverakis and J. C. Sprott.  Hyperlabyrinth chaos: From chaotic walks to spatiotemporal chaos.  Chaos, 17:023110, 2007.
  * [35] L. Zarfaty, A. Peletskyi, E. Barkai, and S. Denisov.  Infinite horizon billiards: Transport at the border between Gauss and Lévy universality classes.  Phys. Rev. E, 100:042140, 2019.
  * [36] G. Teschl.  Ordinary Differential Equations and Dynamical Systems.  American Mathematical Society, 2012.
  * [37] R. Abraham and J. E. Marsden.  Foundations of Mechanics.  Addison-Wesley, 2nd ed. edition, 1987.
  * [38] D. D. Nolte.  Introduction to Modern Dynamics: Chaos, Networks, Space and Time.  Oxford University Press, 2015.
  * [39] H. Goldstein.  Classical Mechanics.  (chap.8), Addison-Wesley, Boston, 1980.
  * [40] V. I. Arnold, V. Kozlov, A. I. Neishtadt, and E. Khukhro.  Mathematical Aspects of Classical and Celestial Mechanics.  Encyclopaedia of Mathematical Sciences. Springer, 3rd ed. edition, 2010\.
  * [41] T. Manos, Ch. Skokos, and C. G. Antonopoulos.  Probing the local dynamics of periodic orbits by the Generalized Alignment Index (GALI) method.  Int. J. Bifurcation and Chaos, 22(09):1250218, 2012.
  * [42] Ch. Skokos, C. Antonopoulos, T. C. Bountis, and M. N. Vrahatis.  Detecting order and chaos in Hamiltonian systems by the SALI method.  J. Phys. A: Math. Gen., 36:6269–6284, 2004.
  * [43] Ch. Skokos, T. C. Bountis, and C. Antonopoulos.  Geometrical properties of local dynamics in Hamiltonian systems: The Generalized Alignment Index (GALI) method.  Physica D, 231:30–54, 2007.
  * [44] J. Hofbauer and K. Sigmund.  Evolutionary Games and Population Dynamics.  Cambridge University Press, 1998.
  * [45] T. Dombre, U. Frisch, J. M. Greene, M. Hénon, A. Mehr, and A. M. Soward.  Chaotic streamlines in the ABC flows.  J. Fluid Mech., 167:353–391, 1986.
  * [46] S. Friedlander, A. Gilbert, and M. Vishik.  Hydrodynamic instability for certain ABC flows.  Geophys. Astrophys. Fluid Dynamics, 73:97–107, 1993.
  * [47] K. Ikeda and K. Matsumoto.  High-dimensional chaotic behavior in systems with time-delayed feedback.  Physica D, 29(1-2):223–235, 1987.
  * [48] Glenn J. Martyna, Michael L. Klein, and Mark Tuckerman.  Nosé–hoover chains: The canonical ensemble via continuous dynamics.  The Journal of Chemical Physics, 97(4):2635–2643, 1992.
  * [49] M. E. Tuckerman, C. J. Mundy, and G. J. Martyna.  On the classical statistical mechanics of non-Hamiltonian systems.  Europhys. Lett., 45:149–155, 1999.
  * [50] M. E. Tuckerman, Y. Liu, G. Ciccotti, and G. J. Martyna.  Non-Hamiltonian molecular dynamics: Generalizing Hamiltonian phase space principles to non-Hamiltonian systems.  J. Chem. Phys., 115:1678–1702, 2001.
  * [51] B. J. Hiley.  From the Heisenberg picture to Bohm: A new perspective on active information and its relation to Shannon information.  In A. Khrennikov, editor, Quantum Theory: Reconsidering of Foundations, pages 141–162. Växjö University Press, 2002.
  * [52] B. J. Hiley.  Active Information and Teleportation.  In Epistemological and Experimental Perspectives on Quantum Physics, pages 113–126, Kluwer, Netherlands, 1999. Edited by D. Greenberger, et al.
  * [53] A. F. Stevenson.  Note on the existence and determination of a vector potential.  Quarterly of Applied Mathematics, 12(2):194–198, 1954.