# Direct learning of home vector direction
for insect-inspired robot navigation

Michiel V.M. Firlefyn, Jesse J. Hagenaars and Guido C.H.E. de Croon This
publication is part of NWA ACT, which is funded by NWO (NWA.1292.19.298). All
authors are with the Micro Air Vehicle Laboratory, Faculty of Aerospace
Engineering, Delft University of Technology, The Netherlands. Correspondence:
<EMAIL_ADDRESS>
###### Abstract

Insects have long been recognized for their ability to navigate and return
home using visual cues from their nest’s environment. However, the precise
mechanism underlying this remarkable homing skill remains a subject of ongoing
investigation. Drawing inspiration from the learning flights of honey bees and
wasps, we propose a robot navigation method that directly learns the home
vector direction from visual percepts during a learning flight in the vicinity
of the nest. After learning, the robot will travel away from the nest, come
back by means of odometry, and eliminate the resultant drift by inferring the
home vector orientation from the currently experienced view. Using a compact
convolutional neural network, we demonstrate successful learning in both
simulated and real forest environments, as well as successful homing control
of a simulated quadrotor. The average errors of the inferred home vectors in
general stay well below the $90^{\circ}$ required for successful homing, and
below $24^{\circ}$ if all images contain sufficient texture and illumination.
Moreover, we show that the trajectory followed during the initial learning
flight has a pronounced impact on the network’s performance. A higher density
of sample points in proximity to the nest results in a more consistent return.
Code and data are available at https://mavlab.tudelft.nl/learning_to_home.

##  I INTRODUCTION

Figure 1: During a learning trajectory (1), we capture omnidirectional images
(2) and train a convolutional neural network (3) to estimate the home vector
direction in a simulated or real environment (4).

Imagine yourself in the middle of a forest. How accurately can you point to
the direction of your home without external aids? Flying insects like honey
bees and wasps have evolved into experts regarding such high-level navigation
tasks, with foraging trips covering distances of up to $13.5$ km from their
nests [1]. Throughout such trips, they retain a sense of the home direction
through path integration [2] based on the Sun’s polarization pattern and the
optical flow across the retinal image [3]. To ensure consistent returns, any
odometric drift accumulated over this journey will have to be corrected for.
To do so, honey bees emerging from the nest for the first time perform
‘learning’ flights around the hive, acquiring visual memories that can be
interpreted in later flights in order to successfully navigate back to the
hive [4, 5].

Two main models of insect-inspired navigation that could accomplish this have
been proposed. In the ‘snapshot model’ [6] a stored snapshot view is compared
with the current view, and the robot moves to minimize their difference [7,
8]. The robot can determine a direction home if the views contain sufficient
similar elements for matching, restricting homing success to a limited
_catchment area_ around the snapshot. For navigating larger distances,
snapshots can be linked together in a sequence, allowing for route following.
The ‘familiarity model’ is an alternative model for route following that
proposes to move into the most familiar direction when returning to the nest
[9]. Typically, a neural network is trained to discriminate between views that
have been perceived before and views that are novel. Various methods, from
rotating physically to using repulsion and attraction, have been proposed to
find the most familiar direction on the way back to the nest [9, 10, 11].
Route following has two main problems: (i) getting off-route, e.g., by one
failure in homing to a snapshot or determining the familiar direction, can
lead to a complete navigation loss, and (ii) the outbound route may be very
tortuous and hence suboptimal for following as an inbound route.

Instead of route following, many insects return to their nest in a straight
trajectory on a novel path [12, 13]. In [13] it has been argued that insects
navigate by means of a nest-centered coordinate frame, always maintaining an
internal representation of a vector to the nest. This view is also supported
by the fact that honey bees convey this type of information to other bees when
communicating about the location of flower fields through waggle dances [14].

Interestingly, honey bees perform ‘learning flights’ of increasing area when
they first leave the hive [15]. After a few learning flights, they are able to
navigate to far-away flower fields and return successfully. We hypothesize
that honey bees use these learning flights to learn a _direct_ mapping from
visual percepts to a vector representation of home, consisting of direction
and possibly distance. Based on this hypothesis, we introduce a novel insect-
inspired navigation algorithm that uses learning flights around the home to
learn how to directly map the currently experienced view to the home
direction.

This algorithm can be abstractly implemented on a robot, making use of the
available insect sensor suite, but not implementing their specific neural
architecture. More specifically, we propose to train a compact convolutional
neural network (CNN) that takes rectified catadioptric omnidirectional images
around the nest as input and returns a unit vector representing the relative
home direction. Fig. 1 shows an overview of this. After learning, a robot can
travel far outside the range of its learning flights and return based on
odometry. Upon return, the trained visual homing network can compensate for
any odometric drift as long as this drift falls within the region explored
during the learning flights. We demonstrate successful learning in both
simulated and real forest environments, as well as successful homing control
of a simulated quadrotor. This approach is very promising for navigation in
large environments, since the learning flight effectively creates a very large
homing catchment area. Depending on the odometry drift, robots can travel much
farther than this area’s diameter. Moreover, in contrast to previous
approaches, our method allows robots to travel back home with odometry on
straight, previously unseen paths. Finally, it is a very memory efficient
navigation method as the homing knowledge is stored in a relatively small
neural network.

##  II RELATED WORK

The first known model to consider insect-inspired navigation is the snapshot
model by Cartwright and Collett [6]. Based on their own experiments, which
involved differing arrays of landmarks, they observed that visual homing in
honey bees relies on visual memories acquired at the nest location, with the
bees primarily concerned with the apparent size and relative bearing of
landmarks. They theorized a framework where the currently perceived horizontal
view is compared to the remembered one, resulting in a home vector that is
oriented towards the hive.

Lambrinos et al. [7] conducted a robotic experiment of the snapshot hypothesis
using a wheeled robot equipped with an omnidirectional camera. Their set-up
consisted of black cylinders positioned on a flat $40\times 40$ m desert
terrain. Although their tested variants address many limitations of the
original snapshot model, they still face challenges when visual noise or
distant landmarks are introduced to the testing area [6]. Cues that have a
large impact on home vector direction, such as landmark occlusion or
mismatches, appear to be the most common obstacles.

Hafner and Möller [8] demonstrated the efficiency and simplicity of learning
the home vector by visually guiding a wheeled robot to its target location in
a simple simulated environment with black cylinder landmarks against a white
background. They employed a snapshot-inspired approach, where a multi-layer
perceptron (MLP) was trained to produce a population encoding, representing
$8$ possible homing directions, based on concatenating the current and home
views as inputs to the network. Instead of concatenating multiple views, we
propose to learn from only the currently experienced view, in more complex
environments.

Baddeley et al. [9] developed a route-following framework that instead learns
a familiarity score based on acquired images along a route. By following the
images with a relatively high familiarity score, the route can be retraced.
While scene familiarity can account for changing landscapes and distant
landmarks, it alone cannot provide an accurate indication of the nest
direction without a continual stop-scan-go routine. This could be addressed by
incorporating a repulsive score, which, when combined with the attractive
score, yields an angular error indicating the most familiar viewing direction
[10]. Van Dalen et al. [16] presented a validation of scene familiarity using
data obtained from a simulated drone, while Gattaux et al. [11] implement
familiarity-based control on a small car-like robot. Although the concept of
attractive/repulsive scene familiarity appears to resolve previous limitations
of homing algorithms, one could question the necessity of inferring two
abstract measures from an image and subsequently integrating them into a
relative angle that indicates the local home vector bearing. Alternatively,
the same outcome could be achieved by directly learning the bearing from the
input image.

Consequently, we present the first visual homing algorithm that allows a robot
to learn the home directions directly during its learning flight, while only
requiring the currently experienced omnidirectional view and odometry to
return to its starting location.

##  III METHODOLOGY

As shown in Fig. 1, we perform experiments in both simulated and real forest
environments, collecting omnidirectional images in regular-grid or spiral
patterns around a nest location. Using these images, we train a compact CNN in
a supervised manner to estimate the home vector direction relative to the
current gaze. In simulation, we investigate the generalization to unseen
locations, and use the trained network for homing of a simulated quadrotor.
The following subsections explain these parts in more detail.

###  III-A Environment

We use Flightmare [17] and Unity to create a $1000\times 1000$ m photo-
realistic forest environment with grass, flowers, different kinds of bushes
and trees, and Perlin-noise terrain elevation (Fig. 3). We spawn quadrotors in
the desired pattern at a fixed height of $1.89$ m above the terrain, and
capture six images (front, back, left, right, dorsal, ventral) with a
$90^{\circ}$-FoV camera with a resolution of $1024\times 1024$ pixels. These
are then combined into a single image using a catadioptric projection [18]
with a latus rectum and distance to the camera image plane of $0.1$ m. To get
an omnidirectional view, we rectify the image from polar to Cartesian
coordinates, setting the gaze angle for this rectification as the angle where
the circular image is cut.

The position of the nest is set at $[538,573]^{T}$. Images are captured in
either a rectangular grid or an Archimedean spiral (representing the honey
bee’s learning flights). Depending on the experiment, the $10\times 10$ m grid
is either spaced at $1$ m or $0.25$ m intervals. The spiral follows radius
$r=0.25\theta$, with $\theta\in[0,8\pi]$ in steps of $0.08\pi$.

We assume two-dimensional movement at the fixed height of $1.89$ m to simplify
rotations and avoid introducing tilt in the acquired images. Furthermore, for
labeling and control purposes, we assume a source of heading information is
available, such as a magnetic or polarization compass, as has been observed in
bees [19, 20]. Lastly, as we focus purely on navigation, we do not consider
collisions, and allow simulated quadrotors/cameras to move and spawn in
objects, possibly leading to (partial) occlusions of the view.

###  III-B Network architecture and training

Figure 2: CNN with kernel and feature map sizes and example input.
Convolutional layers are shown in teal; fully connected (FC) layers in blue.
The gaze is indicated as a black or white arrow and the home vector as a red
arrow. The last nodes have an output value in the range of $[-1,1]$, encoding
the coordinates of the relative home vector. Note that the presented gaze
angle $0^{\circ}$ corresponds to north/upwards.

We employ the compact three-layer CNN shown in Fig. 2, which consists of two
convolutional layers and one fully connected layer. The input image, grayscale
with a shape of $1\times 201\times 1800$
(channels$\times$height$\times$width), is convolved twice with $5\times 5$
kernels and a stride of $4$, doubling the channel count each time. The
resulting $4\times 12\times 112$ feature map is flattened and passed through a
fully connected layer with two outputs, which represent the x- and
y-coordinate of the egocentric home vector relative to the subject’s gaze on
the unit circle. All layers utilize the hyperbolic tangent ($\tanh$)
activation function in their respective intermediate outputs. The resulting
network consists of $11$k parameters (42 kB).

Directly predicting the direction to the nest as a single relative angle did
not give good results, likely because similar angles might have very different
values (e.g., $-175^{\circ}$ and $175^{\circ}$). To resolve this, we decompose
the relative angle into its two-dimensional components on the unit circle, as
can be seen in Fig. 2. Note that the angle formed by the home vector
coordinates is relative to the subject’s gaze, i.e., the middle of the
currently experienced omnidirectional view.

To train the network, we rectify every omnidirectional image collected during
the learning trajectory into $360$ distinct gaze directions ($1^{\circ}$
spacing) and compute the ground truth relative home vectors
$\bm{x}^{\mathrm{rel}}_{\mathrm{home}}$ from geometric relations for all these
gazes:

$\displaystyle\bm{x}^{\mathrm{rel}}_{\mathrm{home}}$
$\displaystyle=R_{z}(-\omega_{\mathrm{gaze}})\bm{x}_{\mathrm{home}}$ (1)
$\displaystyle R_{z}(-\omega_{\mathrm{gaze}})$
$\displaystyle=\begin{bmatrix}\cos\left(\omega_{\mathrm{gaze}}\right)&\sin\left(\omega_{\mathrm{gaze}}\right)\\\
-\sin\left(\omega_{\mathrm{gaze}}\right)&\cos\left(\omega_{\mathrm{gaze}}\right)\end{bmatrix}$
(2)

where $\bm{x}_{\mathrm{home}}$ is the home vector relative to north
($[0,1]^{T}$) and $\omega_{\mathrm{gaze}}$ is the gaze angle relative to
north. The latter can be computed from the gaze vector
$\bm{x}_{\mathrm{gaze}}$ and north vector using $\operatorname{arctan2}$. For
training, we use the mean squared error on the predicted and ground-truth
relative home vector as loss function. Furthermore, we use a batch size of $1$
and a learning rate of $910-4$, and train for a single epoch. Note that these
settings were chosen as a step towards online learning. All training is done
in PyTorch on CPU.111Intel Core i7-8550U running at 1.80 GHz.

As described, training trajectories are either a $10\times 10$ m grid or an
Archimedean spiral. With $100$ sample locations and $360$ rectified images per
location, we get a total of $35,640$ training images (excluding the nest
location), which we shuffle before training. For all results in the next
section, we evaluate the performance for a single northward gaze angle, as the
input to the network has to be rectified, and we cannot show results for all
$360$ possible gaze angles. Section IV-C is an exception: in that case, we
control the quadrotor’s gaze/heading based on the network’s output.

##  IV EXPERIMENTS

###  IV-A Learning and generalization in simulation

Training on the $10\times 10$ m grid and Archimedean spiral flight patterns,
we get the performance as shown in Fig. 3. The bearing map shows the network’s
predicted relative home vector $\bm{\hat{x}}^{\mathrm{rel}}_{\mathrm{home}}$
at every location. If this home vector is not a unit vector, this suggests
that the network has not converged to an exact solution, even if the predicted
angle is correct. We interpret this as the ‘confidence’ of the network in its
prediction. The prediction error, or angle deviation, is computed as the
absolute difference between the predicted and ground-truth relative home
direction, $\hat{\omega}$ and $\omega$, which are computed from
$\bm{\hat{x}}^{\mathrm{rel}}_{\mathrm{home}}$ and
$\bm{x}^{\mathrm{rel}}_{\mathrm{home}}$ with $\operatorname{arctan2}$.

Training the network on images captured in a rectangular grid yields a better
overall performance compared to the spiral case, with respective average
errors222We report average errors in degrees instead of the mean squared error
in m2 used during training to allow better interpretation. of $15.54^{\circ}$
and $19.72^{\circ}$. However, it is worth noting that, in the case of control,
the robot should be able to return to the vicinity of the nest in both
scenarios, since the angular errors of all locations are below $45^{\circ}$.
Furthermore, looking at the shaded green areas in the bearing maps in Fig. 3,
areas near trees that are surrounded by foliage do not show a significant
decline in performance. It seems that important landscape features remain
interpretable, even when they are partially occluded.

TABLE I: Average error $[^{\circ}]$ when evaluating on high-resolution grid. | Trained w. label noise | Trained w.o. label noise
---|---|---
Trained on grid | $16.82\pm 18.98$ | $17.66\pm 18.43$
Trained on spiral | $23.34\pm 21.77$ | $22.77\pm 21.94$

To demonstrate generalization to unseen locations, we capture images in a grid
with decreased spacing ($0.25$ m instead of $1$ m), while maintaining the
$10\times 10$ m grid size. Fig. 4 shows the performance when evaluating on
this higher-resolution grid after training on the lower-resolution grid or the
spiral trajectory. Once again, training on the grid exhibits superior
performance, as can be seen in Table I. When the home vectors are interpreted
as velocity vectors and integrated over the evaluation grid, the stream plots
shown in Fig. 4 can be generated. These plots offer insights into the
potential trajectory a controlled robot would take. As such, their endpoints
are a rendition of the algorithm’s catchment area, delineating the boundaries
for successful homing. Note that these homing streams converge toward a point
that does not precisely coincide with the nest location, as this was excluded
in training. Interestingly, in the case of the spiral (Fig. 4B), this
convergence point lies closer to the hive in all cases. This difference is a
result of the higher training location density near the nest. Moreover, a
circular convergence area takes shape instead of a single point. Such an area
creates less of an ‘artificial’ nest location, as can be seen in the singular
convergence point of the grid case.

While occlusions due to trees do not seem to affect training performance very
much (Fig. 3), they do seem to impact generalization. Looking at the bearing
map and error plots in Fig. 4, we see that whenever the network evaluates a
particular view occluded by tree foliage (home vectors located in green
shading) not encountered during training, the error will be larger. Regions
with high error in the plots of Fig. 4 can therefore be attributed to the
presence of occluding trees. The same holds for the ‘confidence’ of the
network’s prediction, which seems to correlate to some extent with the error
and the tree locations.

To examine the robustness of the algorithm, we train another set of networks
with labels perturbed by Gaussian noise with $\sigma=10^{\circ}$. The results
in Table I show that this has a negligible impact on the evaluation
performance.

Figure 3: Simulation environment and training performance for (A) a $10\times
10$ m grid and (B) an Archimedean spiral. The bearing maps show the predicted
home vectors for all sampled locations, with the nest location in red and
trees shaded in green. Figure 4: Generalization of networks trained on (A) a
$10\times 10$ m grid and (B) an Archimedean spiral to unseen locations,
showing predicted home vectors and resulting stream plots. Locations seen
during training are in magenta, the nest is in red and trees are shaded in
green.

###  IV-B Characterizing network activation

Existing insect-inspired visual homing approaches explicitly define which
landmarks or features to use. On the other hand, we employ a learning approach
without the explicit construction of such features. To get insight into what
information the network uses for prediction, we construct a simplified
environment consisting of flat terrain and three trees of identical appearance
(the ‘landmark array’), and evaluate the activations of networks trained in
this environment.

Fig. 5A/B shows the largest-magnitude gradients of the network’s output with
respect to the activations of the second convolutional layer [21]. Areas in
bright yellow and red show that the network mainly looks at a combination of
horizon and foliage, as well as some tree shadows. There seems to be little
difference in cues for forward/backward (x-coordinate) and left/right
(y-coordinate) motion. While using trees as landmarks is not unlike honey
bees, it is interesting that the networks trained here seem to consider the
contrast of trees against the sky as more useful than the tree stems. This
could be a result of the relatively low-resolution grayscale input, lacking
diversity in shades of green compared to real-world environments, or the fact
that the appearance of the sky is constant in simulation.

Furthermore, we investigate how predictions of a trained network change when
the landmark array is rotated. In Fig. 5, plots A.1/B.1 and A.2/B.2 show
predictions before and after rotation, respectively, with an overview of the
rotation itself shown in the middle. For both A.2 and B.2, the network
predicts a nest location that is shifted downwards by multiple meters. This
shift seems to occur due to the network confusing the original center tree and
rotated right tree (both marked $*$), with the latter closer to the nest.
These results are not unlike those obtained from experiments performed by
Collett and Cartwright [6], where honey bees would maintain relative positions
when landmark arrays were shifted.

Figure 5: Analysis (following [21]) of the highest gradients of the network
output with respect to the activations of the second convolutional layer when
trained on (A) a $10\times 10$ m grid and (B) an Archimedean spiral. Predicted
home vectors after training (A.1/B.1) change when evaluated in an environment
with a landmark array rotated $90^{\circ}$ CCW (A.2/B.2). Original landmarks
are in green, rotated ones in teal, and dominant landmarks are marked with an
asterisk.

###  IV-C Control of a simulated quadrotor

Figure 6: Homing performance when controlling a simulated quadrotor based on a
network trained on (A) a $10\times 10$ m grid and (B) an Archimedean spiral.
Individual trajectories are shown transparent, with the mean trajectory per
homing starting point in solid. Learning, outbound, inbound and visual homing
phases are shown in the inset. Trees are shown in shaded green, and the nest
is marked with a red cross at location $[538,573]^{T}$.

To further demonstrate generalization and usability, we use the output of the
network to control a simulated quadrotor. First, we collect training data in
either a grid or spiral learning flight and train the network as before.
Second, the quadrotor embarks on an outbound flight. Third, we transition to
an inbound phase, using path integration to fly in the nest direction. To
simulate odometric drift, we stop the inbound flight at the edge of the
learning trajectory. Fourth, we transition to visual homing, making use of the
output of the network. The insets in Fig. 6 show these steps for a single run
per homing starting location, whereas the full figure shows the visual homing
part of all runs (transparent trajectories, multiple starting locations). Each
run is made with an individually trained network that takes in an
omnidirectional view at that location, sets the quadrotor’s heading equal to
the predicted home direction, and takes a fixed step of $0.25$ m in that
direction, after which a new omnidirectional view is taken for the next step.
If the quadrotor arrives within $0.25$ m of the nest location, the run is
terminated and counted as a success.

As discussed in Section IV-A, a too low sample density in the vicinity of the
nest leads to an offset convergence point, and we see that in Fig. 6A as well.
This means that none of the runs successfully reach the nest. On the other
hand, when training on a spiral trajectory as in Fig. 6B, $10$ out of $12$
runs are successful, with on average $33.5$ steps taken. So, to ensure
consistent returns, the learning trajectory in the vicinity of the nest should
be of sufficiently high resolution.

###  IV-D Real-world validation

Figure 7: Training performance in two real forest environments, one (A) open
and bright and another (B) covered and dark. The images show grayscale
training samples and an overview of the experimental area, along with
experimental equipment and sample location markers (white arrows). Nest
locations are shown in red.

We show the robustness of our approach by validating it on data collected from
$10\times 10$ m grids in two real forest environments, shown in Fig. 7.
Location A is relatively open and bright, while location B is denser and
darker. Omnidirectional images are taken with a Oneplus One smartphone and a
Kogeto Dot panoramic lens at a height of $\sim 1.44$ m above the ground and a
spacing of $1$ m.

As can be seen, the real-world images contain artifacts not present in their
simulated counterparts; lens flare as well as no/inconsistent illumination
make the real pictures unclear. Moreover, judging from the wave-like
appearance of the images, the employed lens likely does not conform to the
assumption of polar-to-Cartesian rectification. Despite all this, it seems
that our approach transfers well to real-world data, with the bearing map and
training error plot shown in Fig. 7A comparable to the ones in simulation
(e.g., Fig. 3A). In the case of location B, performance decreased due to the
presence of more vegetation and shadows falling from the dense tree canopy on
the forest floor. The respective average errors for locations A and B after
training are $15.72^{\circ}$ and $30.52^{\circ}$, with the former again
comparable to the $15.54^{\circ}$ obtained in simulation.

##  V DISCUSSION & CONCLUSION

We proposed a novel insect-inspired visual homing approach that relies
exclusively on the currently experienced omnidirectional view. Our main
hypothesis posits that bees learn directional cues to their nest during their
learning flights, and use these cues to eliminate odometric drift when
returning from a long foraging journey. We have subsequently demonstrated
successful training of a compact convolutional neural network to predict the
home direction in both simulated and real-world forest environments.
Furthermore, we have investigated the generalization capability of such
networks to locations not seen in training, as well as how they use
environmental cues in their predictions. Successful control of a simulated
quadrotor shows the promise of using our approach for visual homing on real
robots.

This work can be extended in multiple ways, with first and most obvious the
implementation of the network on a real robot with simplified 2D velocity
control. Second, we observed that not shuffling the obtained images prior to
training would lead to catastrophic forgetting. To allow true online during
the learning trajectory (and not have to retain all images in memory), this
will have to be mitigated [22]. Third, we envision to not only learn direction
but also distance during the learning flights. This distance information can
for example be used to adapt the flight velocity, allowing to speed up homing
without losing homing precision. Lastly, to demonstrate the efficiency of our
approach, we should compare it to existing map-building and relocalization
techniques for robots like (visual) SLAM [23] and PoseNet [24].

## References

  * [1] K. von Frisch, _The Dance Language and Orientation of Bees_. Harvard University Press, 1993.
  * [2] T. Haferlach, J. Wessnitzer, M. Mangan, and B. Webb, “Evolving a Neural Model of Insect Path Integration,” _Adaptive Behavior_ , vol. 15, pp. 273–287, 2007.
  * [3] M. V. Srinivasan, S. W. Zhang, M. Lehrer, and T. S. Collett, “Honeybee Navigation En Route to the Goal: Visual Flight Control and Odometry,” _Journal of Experimental Biology_ , vol. 199, pp. 237–244, 1996.
  * [4] E. A. Capaldi, A. D. Smith, J. L. Osborne, S. E. Fahrbach, S. M. Farris, D. R. Reynolds, A. S. Edwards, A. Martin, G. E. Robinson, G. M. Poppy, and J. R. Riley, “Ontogeny of orientation flight in the honeybee revealed by harmonic radar,” _Nature_ , vol. 403, pp. 537–540, 2000.
  * [5] J. Zeil, “Visual homing: An insect perspective,” _Current Opinion in Neurobiology_ , vol. 22, pp. 285–293, 2012.
  * [6] B. A. Cartwright and T. S. Collett, “Landmark learning in bees,” _Journal of comparative physiology_ , vol. 151, pp. 521–543, 1983.
  * [7] D. Lambrinos, R. Möller, T. Labhart, R. Pfeifer, and R. Wehner, “A mobile robot employing insect strategies for navigation,” _Robotics and Autonomous Systems_ , vol. 30, pp. 39–64, 2000.
  * [8] V. V. Hafner and R. Möller, “Learning of Visual Navigation Strategies,” in _Proceedings of the European Workshop on Learning Robots_ , vol. 9, 2001, pp. 47–56.
  * [9] B. Baddeley, P. Graham, P. Husbands, and A. Philippides, “A Model of Ant Route Navigation Driven by Scene Familiarity,” _PLOS Computational Biology_ , vol. 8, p. e1002336, 2012.
  * [10] F. Le Möel and A. Wystrach, “Opponent processes in visual memories: A model of attraction and repulsion in navigating insects’ mushroom bodies,” _PLOS Computational Biology_ , vol. 16, p. e1007631, 2020.
  * [11] G. Gattaux, R. Vimbert, A. Wystrach, J. R. Serres, and F. Ruffier, “Antcar: Simple Route Following Task with Ants-Inspired Vision and Neural Model,” 2023.
  * [12] R. Wehner and S. Wehner, “Insect navigation: Use of maps or Ariadne’s thread?” _Ethology Ecology & Evolution_, vol. 2, pp. 27–48, 1990.
  * [13] B. Webb, “The internal maps of insects,” _Journal of Experimental Biology_ , vol. 222, p. jeb188094, 2019.
  * [14] L. Chittka, _The Mind of a Bee_. Princeton University Press, 2022.
  * [15] J. Degen, A. Kirbach, L. Reiter, K. Lehmann, P. Norton, M. Storms, M. Koblofsky, S. Winter, P. B. Georgieva, H. Nguyen, H. Chamkhi, U. Greggers, and R. Menzel, “Exploratory behaviour of honeybees during orientation flights,” _Animal Behaviour_ , vol. 102, pp. 45–57, 2015.
  * [16] G. J. J. van Dalen, K. N. McGuire, and G. C. H. E. de Croon, “Visual Homing for Micro Aerial Vehicles Using Scene Familiarity,” _Unmanned Systems_ , vol. 06, pp. 119–130, 2018.
  * [17] Y. Song, S. Naji, E. Kaufmann, A. Loquercio, and D. Scaramuzza, “Flightmare: A Flexible Quadrotor Simulator,” in _Proceedings of the 2020 Conference on Robot Learning_. PMLR, 2021, pp. 1147–1157.
  * [18] B. Berenguel-Baeta, J. Bermudez-Cameo, and J. J. Guerrero, “OmniSCV: An Omnidirectional Synthetic Image Generator for Computer Vision,” _Sensors_ , vol. 20, p. 2066, 2020.
  * [19] J. L. Gould, J. L. Kirschvink, and K. S. Deffeyes, “Bees Have Magnetic Remanence,” _Science_ , vol. 201, pp. 1026–1028, 1978.
  * [20] S. Rossel and R. Wehner, “Polarization vision in bees,” _Nature_ , vol. 323, pp. 128–131, 1986.
  * [21] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra, “Grad-CAM: Visual Explanations From Deep Networks via Gradient-Based Localization,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2017, pp. 618–626.
  * [22] J. Kirkpatrick, R. Pascanu, N. Rabinowitz, J. Veness, G. Desjardins, A. A. Rusu, K. Milan, J. Quan, T. Ramalho, A. Grabska-Barwinska, D. Hassabis, C. Clopath, D. Kumaran, and R. Hadsell, “Overcoming catastrophic forgetting in neural networks,” _Proceedings of the National Academy of Sciences_ , vol. 114, pp. 3521–3526, 2017.
  * [23] E. Kruzhkov, A. Savinykh, P. Karpyshev, M. Kurenkov, E. Yudin, A. Potapov, and D. Tsetserukou, “MeSLAM: Memory Efficient SLAM based on Neural Fields,” in _2022 IEEE International Conference on Systems, Man, and Cybernetics (SMC)_ , 2022, pp. 430–435.
  * [24] A. Kendall, M. Grimes, and R. Cipolla, “PoseNet: A Convolutional Network for Real-Time 6-DOF Camera Relocalization,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2015, pp. 2938–2946.


# Stirring the quantum vacuum: Angular Casimir Momentum of a Landau Charge

B.A. van Tiggelen<EMAIL_ADDRESS>Univ. Grenoble Alpes,
CNRS, LPMMC, 38000 Grenoble, France

(August 27, 2024)

###### Abstract

We consider the angular momentum of a charge $q$ rotating in a homogeneous
magnetic field and study the role of the electromagnetic quantum vacuum. Its
orbital angular momentum is caused by the recoil of energetic vacuum photons
that grows as $n^{2}$, i.e. faster than the kinetic angular momentum
$-2n\hbar$ of a Landau level.

Casimir effect and related phenomena, quantum electrodynamics

##  I Introduction

The radiation of the electromagnetic (EM) quantum vacuum is perfectly
isotropic and possesses an energy density $\hbar\omega
d\omega/\pi^{3}c_{0}^{3}$ in the frequency interval $d\omega$. This statement
is Lorentz-invariant milonni . The mystery of its UV-divergence is still one
of the major challenges in physics, but poses in general no problems to
calculate the Casimir force between dielectric or metallic objects Casi . Due
to its perfect isotropy, energy current and momentum density of the quantum
vacuum vanish.

This is no longer true when the vacuum interacts with matter. In a pioneering
work Feigel  Feigel predicted ”Casimir” momentum in bi-anisotropic materials.
Although the Poynting vector of the quantum vacuum still vanishes BartEPJD , a
momentum density emerges that suffers from a divergence at high energies.
Relativistic photons were already known to be relevant for the Lamb shift
milonni , in contrast to Casimir-Polder forces Casi , caused by low-energy
vacuum photons. In a microscopic bi-anisotropic quantum model, the divergency
of EM momentum is removed by mass renormalization kawka , resulting in a
Casimir momentum of the order of $\alpha^{2}$ times the classical Abraham
momentum loudon ; Geert1 . It was also demonstrated for a chiral quantum
particle in a magnetic field Donaire ; Donaire2 .

EM angular momentum is of great recent interest EMAM  and in macroscopic media
even controversial AMdiel , since matter and radiation are hard to
disentangle. In Ref. astrid  Casimir-Polder torques were discussed to stir
vortices in superfluids. How does ”Casimir” angular momentum emerge on a
microscopic scale? The simplest model for which this question can be answered
is the well-known cyclotron problem of a non-relativistic charge $q$ without
spin rotating in a uniform magnetic field $\mathbf{B}_{0}$, in quantum-
mechanics better known as the Landau problem. In the Coulomb gauge and in
Gaussian units, the vector potential is
$\mathbf{A}_{0}(\mathbf{r},t)=\mathbf{B}_{0}(t)\times\mathbf{r}/2$. With the
quantum vacuum, the non-relativistic Hamiltonian reads,

$H=\frac{1}{2\mu}\left(\mathbf{p}-\frac{q}{c_{0}}\mathbf{A}_{0}(\mathbf{r},t)-\frac{q}{c_{0}}\mathbf{A}(\mathbf{r})\right)^{2}+\sum_{{k}\mathbf{\Pi}}\hbar\omega_{k}a^{\dagger}_{k\mathbf{\Pi}}a_{k\mathbf{\Pi}}$
(1)

Here, $\mathbf{p}$ the canonical momentum satisfying
$[r_{i},p_{j}]=i\hbar\delta_{ij}$. It is customary to choose the $z$-axis
along the direction of $q\mathbf{B}_{0}$ and to introduce the cyclotron
frequency $\omega_{c}=qB_{0}/\mu c_{0}$. It is convenient to express the
vector potential of the quantum vacuum in terms of the spherical vector
harmonics $\Phi_{\mathbf{\Pi}}(\mathbf{\hat{k}})$ defined on the unit sphere
in reciprocal space Cohen ,

$\displaystyle\mathbf{A}(\mathbf{r})=\sqrt{2\pi\hbar
c_{0}}\sum_{\mathbf{k}\mathbf{\Pi}}$
$\displaystyle\frac{1}{k^{3/2}}\left(a_{k\mathbf{\Pi}}\exp(-i\mathbf{k}\cdot\mathbf{r})\Phi_{\mathbf{\Pi}}(\mathbf{\hat{k}})\right.$
(2)
$\displaystyle+\left.a^{\dagger}_{k\mathbf{\Pi}}\exp(i\mathbf{k}\cdot\mathbf{r})\bar{{\Phi}}_{\mathbf{\Pi}}(\mathbf{\hat{k}})\right)$

The continuum limit has been taken, $\sum_{\mathbf{k}}\equiv\int
d^{3}\mathbf{k}/(2\pi)^{3}$ and
$[a_{k\mathbf{\Pi}},a^{\dagger}_{k^{\prime}\mathbf{\Pi^{\prime}}}]=\delta_{kk^{\prime}}\delta_{\mathbf{\Pi},\mathbf{\Pi}^{\prime}}$.
The vector index $\mathbf{\Pi}=\\{J,M,p\\}$ summarizes the 3 discrete quantum
numbers of total EM angular momentum ($J$), its $z$-component ($M$) and two
transverse polarizations $p$, the longitudinal vector harmonic being excluded
by the Coulomb gauge. Without the quantum vacuum, the kinetic momentum of the
charge is $\mathbf{p}^{K}=\mu
d\mathbf{r}/dt=\mu[H_{0},\mathbf{r}]/i\hbar=\mathbf{p}-q\mathbf{A}_{0}(\mathbf{r},t)/c_{0}$.
The canonical angular momentum $l_{z}=(\mathbf{r}\times\mathbf{p})_{z}$
commutes with $H_{0}$ and has eigenvalues $(|m|-n)\hbar$. The Hamiltonian
$H_{0}$ has the eigenstates $|n,m,k_{z}\rangle$ of a 2D harmonic oscillator
whose energy levels
$E_{n}=(n+\frac{1}{2})\hbar\omega_{c}+\hbar^{2}k_{z}^{2}/2\mu$ are independent
of $m$ due to gauge invariance Goerbig .

In general, total angular momentum $\mathcal{L}_{z}$ is conserved, equal to
the sum of the canonical angular momentum $l_{z}$ of the charge and the
angular momentum of the transverse EM field $J_{z}$. The angular kinetic
momentum operator of the charge is not conserved and can be written as,

$l_{z}^{\mathrm{K}}=\mathcal{L}_{z}-J_{z}-\frac{q}{c_{0}}(\mathbf{r}\times\mathbf{A})_{z}-l^{L}$
(3)

The ”Lenz” term $l^{L}=\mu\omega_{c}\rho^{2}/2$ is responsible for the
electromotive force in classical electrodynamics, recently discussed quantum-
mechanically Barnettl  as well as key element in the EM momentum controversy
AM . Without the quantum vacuum, the kinetic angular momentum of a Landau
state is given by $\langle
n,m|l_{z}^{K}|n,m\rangle=-(2n+1)\hbar=-2E_{n}(\mathrm{rot})/\omega_{c}$, again
independent on $m$. When the quantum vacuum is included, $\langle
l_{z}^{K}\rangle$ will decay radiatively, but slow enough for $\langle
l_{z}^{K}\rangle$ to be well-defined and modified.

We will use time-dependent perturbation theory to calculate the different
contributions of the quantum vacuum to $l_{z}^{\mathrm{K}}$ for an excited
Landau state, and all proportional to $qB_{0}(t)$ for slow adiabatic changes.
We identify the interaction $W=-(q/\mu
c_{0})\mathbf{p}^{K}\cdot\mathbf{A}(\mathbf{r})$ between rotating charge and
quantum vacuum, and imagine to switch it on slowly like $W(t)=W\exp(\epsilon
t/\hbar)$ at $t_{0}\rightarrow-\infty$, when the total wave function is
assumed to be in the pure state
$|N\rangle=|n,m,k_{z}=0\rangle\otimes|\\{0\\}\rangle$ loudonboek . In this
processus $\mathcal{L}_{z}$ is conserved and equal to its initial value
$(|m|-n)\hbar$. The index $N^{\prime}$ refers to all possible product states
$\\{n,m,k_{z}\\}\otimes\\{n_{k\mathbf{\Pi}}\\}$ of charge plus transverse
photons. Explicit reference will be made nor to the highly degenerated levels
$m$, neither to the momentum $k_{z}$ of the charge along the magnetic field.
The degeneration of the $m$-levels is protected by gauge invariance, and the
impact of photon recoil on the longitudinal displacement is negligible. Since
$W_{NN}=0$, the wave function at $t\geq t_{0}$ is perturbed as,

$\displaystyle|\Psi_{N}(t)\rangle$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{i}{\hbar}\int^{t}_{t_{0}}dt^{\prime}\left[E_{N}+\Delta
E_{N}(t^{\prime})\right]\right)|N\rangle$ $\displaystyle+$
$\displaystyle\sum^{\prime}_{N^{\prime}}|N^{\prime}\rangle\frac{W_{N^{\prime}N}(t)}{E_{N}-E_{N^{\prime}}+i\epsilon}+\cdots$

with $\Delta
E_{N}(t)=\sum^{\prime}_{N^{\prime}}W_{NN^{\prime}}(t)W_{N^{\prime}N}(t)/(E_{N}-E_{N^{\prime}}+i\epsilon)$
the second-order perturbation of the energy level $E_{N}$, the sum
$\sum^{\prime}_{N^{\prime}}$ avoiding the initial level $N$ pertur . For an
excited state, $\int_{t_{0}}^{t}dt^{\prime}\Delta
E_{N}(t^{\prime})=\int_{t_{0}}^{t}dt^{\prime}(E^{L}_{n}+\hbar
A_{n}/2i)-\frac{1}{2}i\hbar\mathcal{N}_{n}$ with $\mathcal{N}_{n}$ a time-
independent normalization of the wave function. Lamb shift and spontaneous
emission rate are

$E^{L}=\hbar\omega_{c}\frac{2\alpha}{3\pi}\,x\log\frac{2}{x}\ \ ,\ \
A_{n}=n\omega_{c}\frac{4\alpha}{3}x$ (4)

with $x=\hbar\omega_{c}/\mu c_{0}^{2}$ and $\alpha=q^{2}/\hbar c_{0}$ the fine
structure constant; $E^{L}$ is equal for all Landau levels whereas radiative
decay is proportional to rotational energy.

In reciprocal space, the transverse EM momentum $J_{z}$ in Eq. (3) is
expressed as
$({J}_{z})(\mathbf{k})_{ij}=-i\hbar\delta_{ij}(\mathbf{k}\times\mathbf{\nabla}_{\mathbf{k}})_{z}-i\hbar\epsilon_{zij}$,
i.e. as a sum of orbital angular momentum and spin Cohen . Neither one of them
behaves as a genuine angular momentum Enk  but this separation is physically
useful. In Hilbert space $J_{z}$ reads,

$J_{z}=\hbar\sum_{\mathbf{k}}\alpha_{i}^{\dagger}(\mathbf{k}){J}_{ij}(\mathbf{k})\alpha_{j}(\mathbf{k})$
(5)

with the photon annihilation operator
$\alpha_{i}(\mathbf{k})=k^{-1}\sum_{\mathbf{\Pi}}a_{k\mathbf{\Pi}}{\Phi}_{\mathbf{\Pi},i}(\mathbf{\hat{k}})$,
and its associated creation operator $\alpha^{\dagger}_{i}(\mathbf{k})$. The
quantum expectation of $J_{z}$ is obtained by inserting the linearly perturbed
eigenfunction (I) on both sides of the matrix element
$\langle\Psi_{N}(t)|J_{z}|\Psi_{N}(t)\rangle$. This creates either a virtual
or real photon with energy $\hbar\omega_{k}$ and angular momentum
$\mathbf{\Pi}^{\prime\prime}$ out of the quantum vacuum at the cost of
canonical angular momentum of the charge. Working out the photon operators,
leaves us with

$\displaystyle\langle J_{z}\rangle$ $\displaystyle=$
$\displaystyle\frac{2\pi\hbar q^{2}}{\mu^{2}c_{0}}\mathrm{e}^{2\epsilon
t}\sum_{\mathbf{k}}\frac{1}{k^{2}}\sum_{\mathbf{k^{\prime}}\mathbf{\Pi}^{\prime}}\sum_{\mathbf{k}^{\prime\prime}\mathbf{\Pi}^{\prime\prime}}\frac{\delta_{kk^{\prime}}\delta_{k^{\prime\prime}k}}{(k^{\prime})^{3/2}(k^{\prime\prime})^{3/2}}$
$\displaystyle\times$ $\displaystyle\langle
n|\mathbf{p}^{K}\cdot\mathbf{\Phi}_{\mathbf{\Pi^{\prime}}}(\hat{\mathbf{k}}^{\prime})\mathrm{e}^{i\mathbf{k^{\prime}r}}\frac{1}{E_{n}-H_{0}^{\prime}-\hbar\omega_{k}-i\epsilon}$
$\displaystyle\times$
$\displaystyle\bar{\mathbf{\Phi}}_{\mathbf{\Pi^{\prime}}i}(\hat{\mathbf{k}}){J}_{ij}(\mathbf{\mathbf{k}}){\mathbf{\Phi}}_{\mathbf{\Pi^{\prime\prime}}j}(\hat{\mathbf{k}})$
$\displaystyle\times$
$\displaystyle\frac{1}{E_{n}-H_{0}^{\prime}-\hbar\omega_{k}+i\epsilon}\mathbf{p}^{K}\cdot\bar{\mathbf{\Phi}}_{\mathbf{\Pi^{\prime\prime}}}(\hat{\mathbf{k}}^{\prime\prime})\mathrm{e}^{-i\mathbf{k^{\prime\prime}r}}|n\rangle$

In principle is $\mathbf{\Pi}^{\prime}=\mathbf{\Pi}^{\prime\prime}$ since the
spherical harmonics are orthogonal eigenfunctions of ${J}_{ij}(\mathbf{k})$.
The mathematics is easier by using their completeness,
$\sum_{\mathbf{\Pi}}\bar{\mathbf{\Phi}}_{\mathbf{\Pi}}(\hat{\mathbf{k}}){\mathbf{\Phi}}_{\mathbf{\Pi}}(\hat{\mathbf{k}}^{\prime})=\delta_{\mathbf{\hat{k}}\hat{\mathbf{k}}^{\prime}}\Delta(\mathbf{\hat{k}})$,
with $\Delta(\mathbf{\hat{k}})$ transverse to $\mathbf{\hat{k}}$ imposed by
the Coulomb gauge. The exponential $\exp(i\mathbf{k\cdot r})$ can be moved by
using the operator identity $\exp(i\mathbf{k\cdot
r})f(\mathbf{p})=f(\mathbf{p}-\hbar\mathbf{k})\exp(i\mathbf{k\cdot r})$, and
which induces a photon recoil in $H_{0}(\mathbf{p})$. With
$\delta_{\mathbf{kk}^{\prime}}=\delta_{\mathbf{\hat{k}\hat{k}}^{\prime}}\delta_{{kk}^{\prime}}/k^{2}$,

$\displaystyle\langle J_{z}\rangle$ $\displaystyle=$
$\displaystyle\frac{2\pi\hbar\mathrm{e}^{2\epsilon
t/\hbar}q^{2}}{\mu^{2}c_{0}}\sum_{\mathbf{k}\mathbf{k}^{\prime\prime}}\langle
n|{p}_{m}^{K}\frac{1}{E_{n}-H_{0}^{\prime}(\mathbf{p}-\hbar\mathbf{k})-\hbar\omega_{k}-i\epsilon}$
$\displaystyle\ \ \
\frac{1}{k}\mathrm{e}^{i\mathbf{k^{\prime}r}}\Delta_{mi}(\mathbf{\hat{k}}){{J}}_{ij}(\mathbf{\mathbf{k}})\Delta_{jl}(\mathbf{\hat{k}}^{\prime\prime})\mathrm{e}^{-i\mathbf{k^{\prime\prime}r}}\delta_{\mathbf{k}\mathbf{k}"}$
$\displaystyle\frac{1}{E_{n}-H_{0}^{\prime}(\mathbf{p}-\hbar\mathbf{k}^{\prime\prime})-\hbar\omega_{k}+i\epsilon}{p}_{l}^{K}|n\rangle$

From this expression the EM spin of the quantum vacuum can be identified as,

$\displaystyle\langle S_{z}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{3}\frac{2\pi\hbar\mathrm{e}^{2\epsilon
t/\hbar}q^{2}}{\mu^{2}c_{0}}\frac{\hbar}{i}\epsilon_{zij}\times$ (6)
$\displaystyle\sum_{\mathbf{k}}\frac{1}{k}\langle
n|{p}_{i}^{K}\frac{1}{|E_{n}-H_{0}^{\prime}-\mathcal{E}(k)+i\epsilon|^{2}}{p}_{j}^{K}|n\rangle$

We have performed the angular integral over $\mathbf{k}$ to eliminate
$\Delta(\mathbf{\hat{k}})$, neglected the photon recoil
$\mathbf{p}^{K}\cdot\hbar\mathbf{k}/\mu$ irrelevant for spin, and defined the
energy $\mathcal{E}(k)=\hbar\omega_{k}+\hbar^{2}k^{2}/2\mu$. The orbital
angular momentum is associated with a differential operator acting on
$\delta_{\mathbf{k}\mathbf{k}"}$, and an integration by parts is imposed to
perform the integral over $\mathbf{k}^{\prime\prime}$. Since this operator
acts only on angles, we find for the orbital angular momentum,

$\displaystyle\langle L_{z}\rangle$
$\displaystyle=\frac{2\pi\hbar\mathrm{e}^{2\epsilon
t/\hbar}q^{2}}{\mu^{2}c_{0}}\sum_{\mathbf{k}}\langle
n|{p}_{m}^{K}\times\frac{\hbar}{i}\epsilon_{zst}$ (7)
$\displaystyle\hat{k}_{t}\nabla_{s}\left(\Delta_{mi}(\mathbf{\hat{k}})\frac{1}{E_{n}-H_{0}^{\prime}(\mathbf{p}-\hbar\mathbf{k})-\hbar\omega_{k}-i\epsilon}\mathrm{e}^{i\mathbf{kr}}\right)$
$\displaystyle\Delta_{il}(\mathbf{\hat{k}})\mathrm{e}^{-i\mathbf{kr}}\frac{1}{E_{n}-H_{0}^{\prime}(\mathbf{p}-\hbar\mathbf{k})-\hbar\omega_{k}+i\epsilon}{p}_{l}^{K}|n\rangle$

As $k_{z}=0$ the kinetic operator $\mathbf{p}^{K}$ is located in the $xy$
plane. It is customary to write
$p^{K}_{x}+(-)ip^{K}_{y}=(2\mu\hbar\omega_{c})^{1/2}c^{({\dagger})}$ in terms
of the raising and lowering operators of the Landau levels, in terms of which
$H_{0}=\hbar\omega_{c}(c^{\dagger}c+\frac{1}{2})$. To evaluate the spin in Eq.
(6) we use
$\epsilon_{zij}p^{K}_{i}f(H_{0})p^{K}_{j}=-i\mu\hbar\omega_{c}(c^{\dagger}f(E_{n-1})c-cf(E_{n+1})c^{\dagger})$.
The first term implies the release of a real photon with energy
$\hbar\omega_{c}$. As $\epsilon\downarrow 0$, this part of $\langle
S_{z}\rangle$ is written as
$(\hbar/2)\int_{t_{0}}^{t}dt^{\prime}\exp(2\epsilon
t^{\prime}/\hbar)\delta(\hbar\omega_{c}-\hbar\omega_{k})$, so that

$\frac{d}{dt}\langle S_{z}\rangle=-\frac{1}{2}A_{n}\hbar$ (8)

with $A_{n}$ defined in Eq. (4). The second term involves virtual photons and
is finite as $\epsilon\downarrow 0$,

$\displaystyle\langle S_{z}\rangle$ $\displaystyle=$
$\displaystyle(n+1)\frac{q^{2}\hbar^{3}\omega_{c}}{3\mu
c_{0}}\int_{0}^{\infty}\frac{kdk}{(\hbar\omega_{c}+\mathcal{E}(k))^{2}}$ (9)
$\displaystyle=$
$\displaystyle\alpha\frac{\hbar}{3\pi}(n+1)\,x\log\frac{2}{x}$

In expression (7) for $\langle L_{z}\rangle$ the derivative
$\nabla_{\mathbf{k}}$ acts on 3 factors inside brackets. Its action on the
factor in the middle, caused by the photon recoil, is a factor $x$ smaller
than the rest. The action on the first factor gives
$\epsilon_{zst}\hat{k}_{t}\nabla_{s}(\Delta_{mi})\Delta_{il}=-k^{-1}\epsilon_{zlt}\hat{k}_{m}\hat{k}_{t}$
and produces an angular momentum $\langle L^{(1)}_{z}\rangle=\langle
S_{z}\rangle$. Finally, the action of $\nabla_{\mathbf{k}}$ on
$\exp(i\mathbf{k}\cdot\mathbf{r})$ leads to the expression

$\displaystyle\langle L^{(2)}_{z}\rangle$ $\displaystyle=$
$\displaystyle\frac{2\pi\hbar^{2}\mathrm{e}^{2\epsilon
t/\hbar}q^{2}}{\mu^{2}c_{0}}\times$ $\displaystyle\sum_{\mathbf{k}}\langle
n|p^{K}_{m}\Delta_{ml}(\mathbf{\hat{k}})\frac{1}{\mathcal{H}_{n}-i\epsilon}(\mathbf{r}\times\mathbf{\hat{k}})_{z}\frac{1}{\mathcal{H}_{n}+i\epsilon}p^{K}_{l}|n\rangle$

where $\mathcal{H}_{n}\equiv
E_{n}-H_{0}-\mathcal{E}(k)+{\hbar\mathbf{p}^{K}\cdot\mathbf{k}}/{\mu}$. This
time, the photon recoil cannot be ignored and we must expand either one of the
denominators which produces an integral $dk\hat{d\mathbf{k}}$ with integrand
of the type

$\epsilon_{zst}\Delta_{ml}(\mathbf{\hat{k}})k_{u}\hat{k}_{t}\times
p^{K}_{m}\frac{1}{\mathcal{H}_{n}-i\epsilon}p^{K}_{u}\frac{1}{\mathcal{H}_{n}-i\epsilon}r_{s}\frac{1}{\mathcal{H}_{n}+i\epsilon}p^{K}_{l}$

Physically, this corresponds to the creation of a virtual photonic mode with
finite orbital angular momentum. In terms of the infinitely degenerated center
$(X,Y)$ of the cyclotron orbit we associate $x=X-p_{y}^{K}/\mu\omega_{c}$ and
$y=Y+p_{x}^{K}/\mu\omega_{c}$. The operators $(X,Y)$ drop out in the vacuum
expectation value since they do not occur in pairs. Upon expressing
$(p^{K}_{x},p^{K}_{y})$ in the operators $c$ and $c^{\dagger}$, the integrand
contains four transition operators. As was the case for $\langle
S_{z}\rangle$, some contribute to spontaneous emission but are seen to be a
factor $x$ smaller than $A_{n}$. We thus focus on terms where the limit
$\epsilon\downarrow 0$ exists. For instance, the sequence

$\langle
n|c\mathcal{H}_{n}^{-1}c^{\dagger}\mathcal{H}_{n}^{-1}c\mathcal{H}_{n}^{-1}c^{\dagger}|n\rangle=\frac{(n+1)(n+2)}{\mathcal{E}(k)^{3}}$

leads to,

$\displaystyle\langle L^{(2)}_{z}\rangle$ $\displaystyle\sim$
$\displaystyle(n+1)(n+2)\frac{2\pi\hbar^{2}q^{2}}{\mu^{2}c_{0}}\frac{\hbar}{\mu}\mu\hbar^{2}\omega_{c}\sum_{\mathbf{k}}\frac{k}{\mathcal{E}(k)^{3}}$
$\displaystyle\approx$ $\displaystyle\frac{\hbar}{8\pi}\alpha(n+1)(n+2)\,x$

Upon collecting all possibilities, performing the angular integral, and adding
the complex conjugate, we obtain

$\displaystyle\langle J_{z}\rangle=\langle
S_{z}\rangle-\frac{4\hbar}{15\pi}\alpha(n+1)(n+4)\,x$ (10)

and ${d}/{dt}\langle J_{z}\rangle=-\hbar A_{n}/2$.

The last but one term in Eq. (3), the longitudinal EM angular momentum $\Delta
L_{\parallel}$, is itself linear in the vacuum field. Its leading quantum
expectation value is obtained using the linear perturbation of the
eigenfunctions and the completeness of the spherical vector harmonics,

$\displaystyle\langle\Psi_{N}|\Delta
L_{\parallel}|\Psi_{N}\rangle=-\frac{q^{2}}{\mu
c_{0}}\frac{4\pi\hbar}{3}\sum_{\mathbf{k}}\frac{1}{k}\epsilon_{zij}\langle
n|r_{i}\frac{1}{\mathcal{H}_{n}+i\epsilon}p^{K}_{j}|n\rangle$

plus its $c.c.$ Here, the $k$-integral diverges in the UV as $dk/k$. We
recognize for large $k$ the form $(\delta\mu/\mu)\langle
n|l^{K}_{z}|n\rangle$, with $\delta\mu$ the well-known Bethe-Kramers mass
renormalization milonni . Upon adding it to the kinetic mass as
$\mathbf{l}^{K}=(\mu+\delta\mu)\mathbf{r}\times d\mathbf{r}/dt$, and upon
subtracting it from the above equation, we obtain,

$\displaystyle\langle\Psi_{N}|\Delta
L_{\parallel}|\Psi_{N}\rangle=-\frac{q^{2}}{\mu
c_{0}}\frac{4\pi\hbar}{3}\sum_{\mathbf{k}}\frac{1}{k}\sum^{\prime}_{n^{\prime}}\frac{i\hbar}{\mu}\times$
$\displaystyle\epsilon_{ijz}\frac{\langle
n|p^{K}_{i}|n^{\prime}\rangle}{\Delta
E_{n^{\prime}n}+\mathcal{E}(k)-i\epsilon}\frac{\langle
n^{\prime}|p^{K}_{j}|n\rangle}{\mathcal{E}(k)}+c.c.$

where we have used the identity $\Delta E_{n^{\prime}n}\langle
n|\mathbf{r}|n^{\prime}\rangle=(i\hbar/\mu)\langle
n|\mathbf{p}^{K}|n^{\prime}\rangle$. The operators $p^{K}_{x}$ and $p^{K}_{y}$
can be expressed in terms of $c^{\dagger},\,c$, which results in,

$\langle\Delta
L_{\parallel}\rangle=\hbar\frac{4\alpha}{3\pi}\,x\log\frac{2}{x}$ (11)

which, like the Lamb shift in energy, is independent on $n$.

The last contribution of the quantum vacuum to the angular momentum is
associated with the Lenz term $l^{T}$ in Eq. (3). The quantum vacuum comes in
via $\mathcal{{N}}_{n}$ and the second term in Eq. (I). With
$l^{T}=\hbar(c^{\dagger}c+1+b^{\dagger}b-ib^{\dagger}c^{\dagger}/2+icb/2)$ in
terms of the raising ($b^{\dagger}$) and lowering ($b$) operator of the
degenerated $m$-levels Goerbig , the $m$-dependence is seen to cancel in the
sum of both terms. The second term equals $\hbar(n+1)$ times

$\displaystyle\frac{4\pi\hbar
q^{2}}{3\mu^{2}c_{0}}\frac{1}{i\hbar}\int_{t_{0}}^{t}\,dt^{\prime}\mathrm{e}^{2\epsilon
t^{\prime}/\hbar}\sum_{\mathbf{k}}\langle
n|{p}^{K}_{j}\frac{1}{\mathcal{H}_{n}+i\epsilon}{p}^{K}_{j}|n\rangle+c.c.$
$\displaystyle=-A_{n}\int_{t_{0}}^{t}\,dt^{\prime}\mathrm{e}^{2\epsilon
t^{\prime}/\hbar}-\frac{2\alpha}{3\pi}(n+1)\,x\log\frac{2}{x}$

and the first is $\hbar$ times

$\displaystyle\frac{4\pi\hbar q^{2}}{3\mu^{2}c_{0}}\mathrm{e}^{2\epsilon
t/\hbar}\sum_{\mathbf{k}}\frac{1}{k}\langle
n|{p}^{K}_{j}\frac{1}{\mathcal{H}_{n}+i\epsilon}(c^{\dagger}c+1)\frac{1}{\mathcal{H}_{n}-i\epsilon}{p}^{K}_{j}|n\rangle$
$\displaystyle=nA_{n}\int_{t_{0}}^{t}\,dt^{\prime}\mathrm{e}^{2\epsilon
t^{\prime}/\hbar}+\frac{2\alpha}{3\pi}(n+1)(n+2)x\log\frac{2}{x}$

In particular, $n^{2}$ terms also cancel and

$\langle\delta l^{T}\rangle=2\langle S_{z}\rangle\,;\frac{d}{dt}\langle\delta
l^{T}\rangle=-\hbar A_{n}$ (12)

By adding up the four contributions $\langle S_{z}\rangle$, $\langle
J_{z}\rangle$, $\langle J_{\parallel}\rangle$ and $\langle\delta l^{T}\rangle$
we find for the total angular momentum of the quantum vacuum

$\displaystyle\langle
J_{z}^{QV}\rangle=\hbar\frac{4\alpha}{3\pi}(n+2)\,x\log\frac{2}{x}-\hbar\frac{4\alpha}{15\pi}(n+1)(n+4)x$

and ${d}\langle J_{z}^{QV}\rangle/dt=-2\hbar A_{n}$. We conclude that the
quantum vacuum achieves an angular momentum that is, in units of $\hbar$,
proportional to $\alpha\times x=\alpha\times\hbar\omega_{c}/\mu
c_{0}^{2}\approx 10^{-12}/$Tesla, which is time-dependent if $B_{0}$ is. In
all momentum integrals the photon momentum $\hbar k$ takes values up to $\mu
c_{0}$ with nonetheless a significant weight of non-relativistic momenta. A
relativistic description of the rotating electron is thus relevant but should
affect only numerical coefficients in Eq. (I). Even in a relativistic picture,
Eq. (3) for the kinetic angular momentum is valid, and $\langle
l^{K}_{z}\rangle$ remains quantized to $-(2n+1)\hbar$ noteNR . In the ground
state, the existence of Casimir angular momentum makes the kinetic angular
momentum slightly more negative than $-\hbar$, in states with large $n$ it
will be slightly less negative than $-(2n+1)\hbar$, the correction growing
like $n^{2}$. Note that the gauge-invariant magnetic moment
$M_{z}=(q/2\mu)l^{K}$ of the rotating charge is subject to the same
correction. Due to the quantum vacuum, the kinetic angular $-(2n+1)\hbar$
decays to the Landau level $n-1$ with rate $A_{n}$ so that $d\langle
l^{K}\rangle/dt=+2\hbar A_{n}$ and $\langle
l_{z}^{K}+J_{z}^{QV}\rangle=\mathcal{L}_{z}$ is conserved in the decay. The
non-relativistic analysis imposes that $E_{n}\ll\mu c_{0}^{2}$, implying, for
an electron in a field of 10 Tesla, that $n\ll 10^{9}$. Pushing our theory to
this extreme synchrotron regime, the relative contribution of Casimir orbital
angular momentum would be of order $10^{-4}$.

##  II Conclusions

The main objective of this work is to establish the existence of angular
momentum of the EM quantum vacuum, induced by the presence of a rotating
charge in a magnetic field. It is instructive to look at the separate
contributions of spin, orbital angular momentum and angular momentum directly
associated with the gauge fields. All are oriented along the magnetic field
and proportional to the product of fine structure constant and the small ratio
of rotational energy to rest energy. Spin and orbital angular momentum decay
in the sam way, their coupling being large, yet the orbital angular momentum
of the quantum vacuum, induced by photon recoil, dominates angular momentum
for highly energetic Landau levels. Casimir angular momentum is governed by
virtual photons with energies up to the rest mass of the charge and, despite
the UV renormalizability of the theory, would merit a relativistic treatment.
A future challenge would be to study EM angular momentum in the fully
relativistic synchrotron problem, or to investigate it for Rydberg orbits.

## References

  * (1) P.W. Milonni, _The Quantum Vacuum_ (Academic Press, San Diego, 1994).
  * (2) K.A. Milton, _The Casimir Effect_ (World Scientific, Singapore, 2001); M. Bordag, G.L. Klimchitskaya, U. Mohideen, and V.M. Mostepanenko, _Advances in the Casimir Effect_ (Oxford University Press, New York, 2009).
  * (3) A. Feigel, Phys. Rev. Lett. 92, 020404 (2004).
  * (4) B.A. van Tiggelen, Eur. Phys. J. D 73, 196 (2019).
  * (5) B. A. van Tiggelen, S. Kawka, G. L. J. A. Rikken, Eur. Phys. J. D 66, 272 (2012).
  * (6) R. Loudon, Fortschr. Phys. 52, 1134 (2004).
  * (7) G.L.J.A. Rikken, B.A. van Tiggelen, Phys. Rev. Lett. 107, 170401 (2011).
  * (8) M. Donaire, B.A. van Tiggelen, G.L.J.A. Rikken, Phys. Rev. Lett. 111, 143602 (2013).
  * (9) M. Donaire, arXiv:1907.13518 [hep-ph] (2019).
  * (10) K.Y. Bliokh and F. Nori, Phys. Rep. 592, 1-38 (2015).
  * (11) M. Kristensen and J.P. Woerdman, Phys. Rev. Lett. 72, 2171 (1994).
  * (12) F. Impens, A.M. Contreras-Reyes, P.A. Maia Neto, D.A.R. Dalvit, R. Guérout, A. Lambrecht and S. Reynaud, EPL 92(4),40010 (2010).
  * (13) C. Cohen-Tannoudji, J. Dupont-Roc, G. Grynberg, _Photons and Atoms_ (Wiley, 1989).
  * (14) M.O. Goerbig, Quantum Hall Effects, in: Les Houches Summer School XCI, _Ultracold Gases and Quantum Information_(Oxford, University Press, 2011).
  * (15) C.R. Greenshields, R.L. Stamps, S. Franke-Arnold, S.M. Barnett, Phys. Rev. Lett. 113, 240404 (2014).
  * (16) E.A. Hinds and S.M. Barnett, Phys. Rev. Lett. 102, 050403 (2009).
  * (17) R. Loudon, T _he Quantum Theory of Light_ (Oxford University, 1991).
  * (18) We ignored the perturbation $|\delta\Psi_{N}\rangle=\sum^{\prime}_{N^{\prime}N^{\prime\prime}}|N^{\prime}\rangle{W_{N^{\prime}N^{\prime\prime}}W_{N^{\prime\prime}N}}/({(E_{NN^{\prime}}+2i\epsilon)(E_{NN^{\prime\prime}}+i\epsilon)})$. For one-photon processes, the matrix element $\langle n,\\{0\\}|l^{K}_{z}|\delta\Psi_{n\\{0\\}}\rangle$ vanishes.
  * (19) S. van Enk and G. Nienhuis, EPL 25, 497 (1994).
  * (20) The eigenfunctions are the same for the relativistic Hamiltonian without spin, and thus the Lenz moment $\langle l^{T}\rangle$ remains unaltered, that in turn determines $l^{K}_{z}$Cohen .


# Bistability and oscillatory behaviours of cyclic feedback loops

Jules Guilberteau  Sorbonne Université, CNRS, Université Paris Cité, Inria,
Laboratoire Jacques-Louis Lions (LJLL), F-75005 Paris, France.
<EMAIL_ADDRESS>
###### Abstract

In this paper, we study the stability of an Ordinary Differential Equation
(ODE) usually referred to as Cyclic Feedback Loop, which typically models a
biological network of $d$ molecules where each molecule regulates its
successor in a cycle ($A_{1}\rightarrow A_{2}\rightarrow\ldots\rightarrow
A_{d-1}\rightarrow A_{d}\rightarrow A_{1}$). Regulations, which can be either
positive or negative, are modelled by increasing or decreasing functions. We
make a complete analysis of this model for a wide range of functions
(including affine and Hill functions) by determining the parameters for which
bistability and oscillatory behaviours arise. These results encompass previous
theoretical studies of gene regulatory networks, which are particular cases of
this model.

Keywords : Gene regulatory network, Repressilator, Toggle switch, Stability
analysis, Multistability, Periodic orbit.

##  1 Introduction

We aim to characterise the stability of the ODE system

$\displaystyle\begin{cases}\dot{x}_{1}=\alpha_{1}f_{1}(x_{d})-x_{1}\\\
\dot{x}_{2}=\alpha_{2}f_{2}(x_{1})-x_{2}\\\ \vdots\\\
\dot{x}_{d}=\alpha_{d}f_{d}(x_{d-1})-x_{d}\end{cases}$ (1)

where $f_{1},...,f_{d}\in\mathcal{C}^{1}(\mathbb{R}_{+},\mathbb{R}_{+})$ (with
$\mathbb{R}_{+}=[0,+\infty)$) are non-negative functions, at least one of them
is bounded, and $\alpha_{1},...,\alpha_{d}$ are positive parameters.
Throughout this paper, we use the convention $x_{0}=x_{d}$, which allows us to
write (1) under the compacted form

$\quad\forall
i\in\\{1,...,d\\},\quad\dot{x}_{i}=\alpha_{i}f_{i}(x_{i-1})-x_{i}.$

Figure 1: Schematic representation of the cyclic feedback loop (1): The blue
circles represent the molecules of the network, and the arrows between them
regulation, which can be positive ($f_{i}$ increasing) or negative ($f_{i}$
decreasing).

This model is a generalisation of a gene regulatory network initially proposed
by Goodwin [7, 6], usually referred to as Cyclic feedback loop, which
represents interactions between genes, mRNAs, enzymes and proteins called
repressors which have the ability to inhibit the expression of some genes. In
system (1), $x_{1},...,x_{d}$ represent the concentration of each of the
molecules involved in the network (denoted $A_{1}$, …, $A_{n}$), and
$f_{1},...,f_{d}$ the regulation between them. The system is assumed to be
cyclic ($A_{i}$ regulates $A_{i+1}$ and only $A_{i+1}$, as illustrated by
Figure 1), and each regulation can be positive ($f_{i}$ increasing) or
negative ($f_{i}$ decreasing). The relevance of these cyclic models has been
established in [1] and [5] where some theoretical predictions (oscillatory
phenomena and bistability) have been observed experimentally. This highlights
the importance of understanding the dynamical behaviour of such systems i.e.
determining the number of stable equilibrium points and their basins of
attraction, as well as the possible existence of periodic solutions or chaotic
behaviours.

System (1) has been, for some specific choices of $f_{i}$, the subject of
several theoretical studies [2, 14, 16, 17, 18, 19]. In these papers,
restrictions on the functions $f_{i}$ were notably imposed by the necessity to
compute the value of the equilibrium points of the system, which is intricate
when more than two functions are not affine, and are not identical. In the
present paper, we follow a method initiated by Cherry and Adler [4] allowing
to avoid explicitly computing the equilibrium points. In the two dimensional
case, which writes

$\displaystyle\begin{cases}\dot{x}_{1}=\alpha_{1}f_{1}(x_{2})-x_{1}\\\
\dot{x}_{2}=\alpha_{2}f_{2}(x_{1})-x_{2}\end{cases},$ (2)

the following results have already been established:

  1. 1.

If $f_{1}$ and $f_{2}$ are both increasing or both decreasing, and if

$\displaystyle\underset{x>0}{\sup}\left\lvert\frac{xf_{1}^{\prime}(x)}{f_{1}(x)}\right\rvert\underset{x>0}{\sup}\left\lvert\frac{xf_{2}^{\prime}(x)}{f_{2}(x)}\right\rvert>1,$
(3)

then there exist values of
$(\alpha_{1},\alpha_{2})\in{\mathbb{R}_{+}^{*}}^{2}$ such that system (2) is
multistable i.e. there exist at least two equilibrium points which are
asymptotically stable [4].

  2. 2.

If $f_{1}$ and $f_{2}$ are both increasing or both decreasing, and if

$\frac{1}{\sqrt{\lvert
f_{1}^{\prime}\rvert}}\quad\text{and}\quad\frac{1}{\sqrt{\lvert
f_{2}^{\prime}\rvert}}\quad\text{are strictly convex },$ (4)

then system (2) is either monostable or bistable, i.e. there exist exactly one
or exactly two equilibrium points which are asymptotically stable. Moreover,
it is possible to determine, up to a set of measure zero, the set of
parameters for $(\alpha_{1},\alpha_{2})$ for which the system is monostable
and the set of parameters for which it is bistable [8].

  3. 3.

All the solutions to system (2) converge (even without assuming any
monotonicity).

This last result is a direct application of the Poincaré-Bendixson theorem and
the Dulac-Bendixson theorem [15].

A natural question at this stage is which of these properties generalise to
higher dimensions ($d\geq 3$). A major result was achieved by Mallet-Paret and
Smith in [13], showing that the Poincaré-Bendixson theorem can be adapted to
monotone feedback systems, including (1). Furthermore, a series of results of
Hirsch [9, 11, 10] generalised and summarised in [18], has shown that in the
case where an even number of functions $f_{i}$ is decreasing, the solutions of
(1) converge to an equilibrium point for almost every initial condition. It is
well-known that this property does not hold when the number of decreasing
functions is odd (in which case the system is often called ‘repressilator’),
and that there can there exist stable orbits [3, 14].

In this paper, we prove the following result, which is a generalisation of our
previous paper on the two-dimensional case [8]:

######  Theorem 1.

Let us assume that
$f_{1},...,f_{d}\in\mathcal{C}^{3}(\mathbb{R}_{+},\mathbb{R}_{+}^{*})$ are
monotonous, non-negative, and that (at least) one of them is bounded.
Moreover, let us assume that the functions $\frac{1}{\sqrt{\lvert
f^{\prime}_{1}\rvert}},...,\frac{1}{\sqrt{\lvert f^{\prime}_{d}\rvert}}$ are
defined and convex, and that (at least) one of them is strictly convex.
Lastly, let us denote by $n$ the number of these functions which are
decreasing, and let

$D:=\prod\limits_{k=1}^{d}{\underset{x>0}{\sup}\left\lvert\frac{xf_{k}^{\prime}(x)}{f_{k}(x)}\right\rvert}\in(0,+\infty]$

  1. 1.

If $n$ is even, then

    1. (i)

If $D<1$, then for any $\alpha\in\left(\mathbb{R}_{+}^{*}\right)^{d}$ system
(1) has a unique equilibrium point which is globally asymptotically stable.

    2. (ii)

If $D>1$, then there exists a non-empty set
$A_{\textrm{bis}}\subset\left(\mathbb{R}_{+}^{*}\right)^{d}$ such that

       * •

If $\alpha\in A_{\textrm{bis}}$, then system (1) has exactly two
asymptotically stable equilibria, and the union of their basins of attraction
is a dense open subset of $\mathbb{R}_{+}^{d}$, with a complement of Lebesgue
measure zero.

       * •

If $\alpha\in\overline{A_{\textrm{bis}}}^{C}$, then system (1) has a unique
equilibrium point which is globally asymptotically stable.

    3. (iii)

If $d\geq 5$ and $D>\frac{1}{\cos\left(\frac{2\pi}{d}\right)^{d}}$, then there
exists a non-empty set $A_{\textrm{per}}\subset A_{\textrm{bis}}$ such that if
$\alpha\in A_{\textrm{per}}$, then system (1) has periodic solutions.

  2. 2.

If $n$ is odd, then system (1) has a unique equilibrium point. Moreover, if
$d\geq 3$, then

    1. (i)

If $D<\frac{1}{\cos\left(\frac{\pi}{d}\right)^{d}}$, then for any
$\alpha\in\left(\mathbb{R}_{+}^{*}\right)^{d}$, this equilibrium point is
asymptotically stable, and all the solutions of (1) either converge to this
point or to a periodic orbit.

    2. (ii)

If $D>\frac{1}{\cos\left(\frac{\pi}{d}\right)^{d}}$, then there exists
$A_{\textrm{unst}}\subset\left(\mathbb{R}_{+}^{*}\right)^{d}$ a non-empty set
such that

       * •

If $\alpha\in A_{\textrm{unst}}$, then this equilibrium point is
asymptotically unstable, and there exists a finite number of periodic
solutions, among which at least one is asymptotically stable. Moreover, the
set of initial conditions for which the solution converges to a periodic
solution is a dense open subset of $\mathbb{R}_{+}^{d}$, and its complement,
which is the set of initial conditions for which the solution converges to the
equilibrium point, has Lebesgue measure zero.

       * •

If $\alpha\in\overline{A_{\textrm{unst}}}^{C}$, this equilibrium point is
asymptotically stable, and all the solutions of (1) either converge to this
point or to a periodic orbit.

In each of these cases, the sets $A_{\textrm{bis}}$, $A_{\textrm{per}}$,
$A_{\textrm{unst}}$ can be explicitly expressed, as we will show in Sections 3
and 4. Moreover, note that any Hill function (even shifted), i.e. function of
the form $x\mapsto\frac{1+\lambda x^{r}}{1+x^{r}}$, with
$\lambda\in\mathbb{R}_{+}\backslash\\{1\\}$ and $r\geq 1$, as well as linear
functions, satisfy the convexity hypothesis of this theorem, which means that
this result encompasses the other theoretical studies mentioned above [2, 14,
16, 17, 18, 19].

It is worth noting that, when $n$ is even, the set of initial conditions for
which the solution converges to a periodic orbit has Lebesgue measure zero,
and can therefore hardly be reached numerically. Nevertheless, we highlight
this result since, up to our knowledge, the question of the existence of such
periodic solution remained open, as mentioned in [18]. Moreover, the proof of
this result, which uses the stable manifold theorem, seems to us non-trivial
and worthwhile. We also note that condition 1 (iii) is sufficient, but perhaps
not necessary for the existence of periodic solution: in particular, this
question for $d\in\\{3,4\\}$ remains open. Lastly, we do not know if periodic
solutions do exist under the hypotheses of 2 (i) and in the second point of 2
(ii), but our study does not rule out this possibility.

This article, entirely dedicated to the proof of Theorem 1, is organised as
follows: after characterising the equilibrium points of system (1) in Section
2, we prove the theorem when $n$ is even (Section 3), before dealing with the
case where $n$ is odd, which is simpler, in the last section.

##  2 Characterisation of fixed points

Throughout this note, we assume that
$f_{1},...,f_{d}\in\mathcal{C}^{3}(\mathbb{R}_{+},\mathbb{R}_{+}^{*})$ are
non-negative, monotonous and that at least one of these functions is bounded.
Note that, according to the regularity of these functions, the Cauchy-
Lipschitz theorem ensures the local existence and the uniqueness of the
solution of this equation for any initial condition
$x_{0}\in\mathbb{R}^{d}_{+}$. Furthermore, the positivity of the functions
guarantees that the solutions remain in $\mathbb{R}^{d}_{+}$. Lastly, using
the fact that one of them is bounded, we easily prove the global existence of
solutions and the existence of a compact attractor set,i.e. the existence of a
compact set $K\subset\mathbb{R}^{d}_{+}$ such that for any initial condition
$x_{0}\in\mathbb{R}_{+}^{d}$, there exists $T\geq 0$ such that $x(t)\in K$ for
all $t\geq T$.

We now make the additional assumption that $f_{1},...,f_{d}$ are
$\gamma^{1/2}-$convex, and that at least one of them is strictly
$\gamma^{1/2}-$convex, i.e. they satisfy the following definition:

######  Definition 1 (($\mathbf{\gamma^{1/2}-}$convexity)).

Let $f\in\mathcal{C}^{3}\left(\mathbb{R}_{+},\mathbb{R}_{+}\right)$ a non-
negative and monotonous function. We say that $f$ is (strictly)
$\gamma^{1/2}-$ convex if $\lvert f^{\prime}\rvert>0$ and
$\frac{1}{\sqrt{\lvert f^{\prime}\rvert}}$ is (strictly) convex.

Note that this definition can be related to the definition of the Schwartzian
derivative of $f$ defined by
$S(f)=\frac{f^{\prime\prime\prime}}{f^{\prime}}-\frac{3}{2}\left(\frac{f^{\prime\prime}}{f^{\prime}}\right)^{2}$
by noting that $(\frac{1}{\sqrt{\lvert
f^{\prime}\rvert}})^{\prime\prime}=-\frac{1}{2}\frac{1}{\sqrt{\lvert
f^{\prime}\rvert}}S(f)$.

We relate here the equilibrium points of system (1) to the fixed points of an
auxiliary function $\tilde{f}$. The $\gamma^{1/2}-$convexity of the functions
$f_{i}$ ensures that the number of equilibrium points cannot exceed three, and
provides a criterion which characterises this exact number of equilibria. A
similar approach was used, (with the Schwartzian derivative) for a particular
case of this system in [14].

We start by recalling some key properties of the $\gamma^{1/2}-$convexity,
which have been established in [8].

######  Proposition 1.

Let $f,g$ be two $\gamma^{1/2}-$convex functions, $c>0$.

  1. (i)

$f\circ g$ and $cf$ are $\gamma^{1/2}-$convex. Moreover, if $f$ or $g$ is
strictly $\gamma^{1/2}-$convex, then $f\circ g$ is strictly
$\gamma^{1/2}-$convex.

Let us now assume that $f$ is strictly $\gamma^{1/2}-$convex. Then:

  1. (ii)

$cf$ is strictly $\gamma^{1/2}-$convex.

  2. (iii)

$f$ has at most three fixed points.

  3. (iv)

If all the fixed points $x$ of $f$ satisfy $f^{\prime}(x)<1$, then $f$ has a
unique fixed point.

  4. (v)

If there exists a fixed point of $f$ (denoted $x$) such that
$f^{\prime}(x)>1$, then $f$ has exactly three fixed points and the other two
fixed points (denoted $y$, $z$) satisfy $f^{\prime}(y)<1$ and
$f^{\prime}(z)<1$.

######  Example.

For any $r\geq 1$, $a,b,c,d\geq 0$ such that $ad-bc\neq 0$, the function
$x\mapsto\frac{ax^{r}+b}{cx^{r}+d}$ is $\gamma^{1/2}-$convex. Moreover, if
$r>1$, then it is strictly $\gamma^{1/2}-$convex. In particular, affine
functions and Hill functions are $\gamma^{1/2}-$convex.222In the Appendix of
[8] we show that many other usual sigmoid functions are strictly
$\gamma^{1/2}-$convex.

By definition, the point
$\bar{x}=(\bar{x}_{1},...,\bar{x}_{d})\in\mathbb{R}_{+}^{d}$ is an equilibrium
point of (1) if and only if

$\displaystyle\begin{cases}\bar{x}_{1}=\alpha_{1}f_{1}(\bar{x}_{d})\\\
\bar{x}_{2}=\alpha_{2}f_{2}(\bar{x}_{1})\\\ \vdots\\\
\bar{x}_{d}=\alpha_{d}f_{d}(\bar{x}_{d-1})\end{cases}\Longleftrightarrow\quad\begin{cases}\bar{x}_{1}=\alpha_{1}f_{1}(\bar{x}_{d})\\\
\bar{x}_{2}=\alpha_{2}f_{2}(\bar{x}_{1})\\\ \vdots\\\
\bar{x}_{d-1}=\alpha_{d-1}f_{d-1}(\bar{x}_{d-2})\\\
\bar{x}_{d}=\alpha_{d}f_{d}\circ\alpha_{d-1}f_{d-1}\circ...\circ\alpha_{1}f_{1}(\bar{x}_{d})\\\
\end{cases}.$

Thus, the number of equilibrium points of $\eqref{ODE dD}$ is equal to the
number of fixed points of
$\tilde{f}:=\alpha_{d}f_{d}\circ...\circ\alpha_{1}f_{1}$. Since
$\tilde{f}(0)>0$ and $\tilde{f}$ is bounded, it proves in particular that (1)
has at least one equilibrium point. Note that, according to the first two
properties of Proposition 1, $\tilde{f}$ is strictly $\gamma^{1/2}-$ convex,
and a direct computation shows that

$\tilde{f}^{\prime}(\bar{x}_{d})=\prod\limits_{i=1}^{d}{\alpha_{i}f_{i}^{\prime}(\bar{x}_{i-1})}.$

We can thus apply the fourth and the fifth properties of Proposition 1 to
$\tilde{f}$ to derive the following lemma:

######  Lemma 1.

Let us denote, for all $\alpha,x\in\mathbb{R}_{+}^{d}$,
$p_{x}^{\alpha}=\prod\limits_{i=1}^{d}{\alpha_{i}f_{i}^{\prime}(x_{i-1})}$.

  * •

If all the equilibrium points of system (1) satisfy $p_{\bar{x}}^{\alpha}<1$,
then this system has a unique equilibrium point.

  * •

If there exists an equilibrium point of system (1) (denoted $\bar{x}$) such
that $p_{\bar{x}}^{\alpha}>1$, then this system has exactly three equilibrium
points, and the other two points (denoted $\bar{y},\bar{z}$) satisfy
$p_{\bar{y}}^{\alpha}<1$, $p_{\bar{z}}^{\alpha}<1$.

Hence, the value of $p^{\alpha}_{\bar{x}}$ characterises the number of fixed
points the system has. In the following section, we show that it also
determines the dimension of the basin of attraction of $\bar{x}$.

##  3 Even number of decreasing functions

In the case where $n$ is even, one easily checks that system (1) is an
irreducible type K monotone system in the sense defined in [18]. As seen in
the previous section, this system has a finite number of equilibrium points
(at most three), and a compact attractor set: thus, we can apply Theorem 2.5
and Theorem 2.6 of [18] which prove that the union of the basins of attraction
of the equilibrium points is dense, and that the complement of this set has
Lebesgue measure zero.

In this section, we complete this result in two ways:

  * •

We determine, for given functions $f_{1},...,f_{d}$, a set of parameters
$A_{\textrm{bis}}$ such that system (1) is bistable (i.e. has exactly two
asymptotically stable equilibrium points) if $\alpha\in A_{\textrm{bis}}$, and
monostable (i.e. has exactly one asymptotically stable equilibrium point) if
$\alpha\in\overline{A_{\textrm{bis}}}^{C}$.

  * •

We determine a set $A_{\textrm{per}}\subset A_{\textrm{bis}}$ such that system
(1) has some periodic solutions if $\alpha\in A_{\textrm{per}}$.

Note that this last point does not mean that periodic solutions do not exist
when $\alpha\notin A_{per}$, and that, in all cases, the set of initial
conditions for which the solution converges to a periodic solution has
Lebesgue measure zero (as a corollary of [18]).

In order to prove these two points, we determine the dimension of the basin of
attraction of an equilibrium point $\bar{x}$, as a function of
$p_{\bar{x}}^{\alpha}$. For any equilibrium point $\bar{x}$, we denote its
basin of attraction $B_{\bar{x}}.$ Our reasoning is based on the stable
manifold theorem (the proof of which can be found for instance in [15]), that
we recall:

######  Theorem ((Stable manifold)).

Let $F\in\mathcal{C}^{1}(\mathbb{R}^{d},\mathbb{R}^{d})$ a vector field, and
let $\bar{x}\in\mathbb{R}^{d}$ such that $F(\bar{x})=0_{\mathbb{R}^{d}}$. If
$\bar{x}$ is a hyperbolic equilibrium point, i.e. if all the eigenvalues of
$\textrm{Jac}\,F(\bar{x})$ have a non-zero real part, then the basin of
attraction of $\bar{x}$ is a manifold of dimension $m$, where $m$ is the
number of eigenvalues of $\textrm{Jac}\,F(\bar{x})$ with a negative real part.

In order to apply this theorem, we need to compute the eigenvalues of the
Jacobian matrix associated to system (1) which writes, at a given point
$x\in\mathbb{R}^{d}$

$\displaystyle
M_{x}^{\alpha}=\begin{pmatrix}-1&0&\cdots&0&\alpha_{1}f_{1}^{\prime}(x_{d})\\\
\alpha_{2}f_{2}^{\prime}(x_{1})&-1&0&\cdots&0\\\
0&\ddots&\ddots&\ddots&\vdots\\\ \vdots&\ddots&\ddots&\ddots&0\\\
0&\cdots&0&\alpha_{d}f_{d}^{\prime}(x_{d-1})&-1\end{pmatrix}.$

Thus, the characteristic polynomial of $M_{x}^{\alpha}$ is easily computed to
be

$(-1)^{d}\left((\lambda+1)^{d}-\prod\limits_{i=1}^{d}{\alpha_{i}f_{i}^{\prime}(x_{i-1})}\right)=(-1)^{d}\left((\lambda+1)^{d}-p_{x}^{\alpha}\right)$

Since $n$ is even, $p_{x}^{\alpha}>0$ and hence the spectrum of
$M_{x}^{\alpha}$ is given by

$\mathrm{Sp}\left(M_{x}^{\alpha}\right)=\left\\{(p_{x}^{\alpha})^{1/d}e^{2k\pi
i/d}-1,k\in\\{0,...,d-1\\}\right\\}.$

We deduce that

  1. (i)

If $p_{x}^{\alpha}<1$, then all the eigenvalues of $M_{x}^{\alpha}$ have a
negative real part.

  2. (ii)

If $d\in\\{3,4\\}$ and $p_{x}^{\alpha}>1$, or if $d\geq 5$ and
$p_{x}^{\alpha}\in\left(1,\frac{1}{\cos\left(\frac{2\pi}{d}\right)^{d}}\right)$,
then $M_{x}^{\alpha}$ has exactly $d-1$ eigenvalues with a negative real part,
and one with a positive real part.

  3. (iii)

If $d\geq 5$, and if
$p_{x}^{\alpha}>\frac{1}{\cos\left(\frac{2\pi}{d}\right)^{d}}$, then
$M_{x}^{\alpha}$ has at most $d-3$ eigenvalues with a negative real part.
Moreover, if $d\leq 8$, or $d\in\llbracket 4j+1,4j+4\rrbracket$
($j\in\mathbb{N}\backslash\\{0,1\\}$), and for any $k\in\\{2,...,j\\}$,
$p_{x}^{\alpha}\neq\frac{1}{\cos\left(\frac{2\pi k}{d}\right)^{d}}$, then all
eigenvalues of $M_{x}^{\alpha}$ have a non-zero real part (where for all
$a,b\in\mathbb{R}$, $a<b$, $\llbracket a,b\rrbracket=[a,b]\cap\mathbb{N}$.)

Thus, the stable manifold theorem yields

######  Lemma 2.

Let us assume that $n$ is even, and let $\bar{x}\in\mathbb{R}^{d}$ be an
equilibrium point of (1).

  1. (i)

If $p_{\bar{x}}^{\alpha}<1$, then $\mathrm{dim}\left(B_{\bar{x}}\right)=d$,
i.e. $B_{\bar{x}}$ is an open set.

  2. (ii)

If $d\in\\{3,4\\}$ and $p_{\bar{x}}^{\alpha}>1$, or if $d\geq 5$ and
$p_{\bar{x}}^{\alpha}\in\left(1,\frac{1}{\cos\left(\frac{2\pi}{d}\right)^{d}}\right)$,
then $\mathrm{dim}\left(B_{\bar{x}}\right)=d-1$.

  3. (iii)

If $d\geq 5$,
$p_{\bar{x}}^{\alpha}>\frac{1}{\cos\left(\frac{2\pi}{d}\right)^{d}}$, and
$p_{\bar{x}}^{\alpha}\notin S_{d}$, with

$\displaystyle S_{d}:=\begin{cases}\quad\emptyset\quad&\textrm{if}\quad
d\in\\{5,6,7,8\\}\\\ \left\\{\frac{1}{\cos\left(\frac{2\pi
k}{d}\right)^{d}},k\in\\{2,...j\\}\right\\}\quad&\textrm{if}\quad
d\in\llbracket 4j+1,4j+4\rrbracket,\quad j\geq 2\end{cases},$

then $\mathrm{dim}\left(B_{\bar{x}}\right)\leq d-3$.

Moreover, we easily check that, in all cases, $M_{\bar{x}}^{\alpha}$ has an
odd number of eigenvalues with a positive real part. Thus, since $n$ is even,
we get from the main theorem of [13] that any solution converges to an
equilibrium point or to a periodic orbit.

We will now use Lemmas 2 and 3 to prove Theorem 1 in a more precise form which
specifies the sets $A_{\mathrm{bis}}$ and $A_{\mathrm{per}}$, in the case
where $n$ is even. Before stating it, we give a last lemma linking the
dimension of the basin of attraction of the unstable equilibrium point to the
existence of divergent solutions (which thus converge to periodic solutions),
in the bistable case.

######  Lemma 3.

Let ‘$\dot{x}=F(x)$’ be an ODE which has exactly three equilibrium points,
(denoted $\bar{x}$, $\bar{y}$, $\bar{z}$) assumed hyperbolic, and let us
assume that for all initial condition $x_{0}\in\mathbb{R}^{d}$, the solution
of this ODE is defined on $\mathbb{R}_{+}$ and converges. If $\bar{y}$,
$\bar{z}$ are asymptotically stable, then
$\mathrm{dim}\left(B_{\bar{x}}\right)=d-1.$

###### Proof.

First, let us note that, according to the stable manifold theorem,
$B_{\bar{y}}$ and $B_{\bar{z}}$ are two open sets, and that $B_{\bar{x}}$ is a
manifold. Since, by hypothesis, all the solutions converge,
$\mathbb{R}_{+}^{d}=B_{\bar{x}}\cup B_{\bar{y}}\cup B_{\bar{z}}$, which means
that $B_{\bar{x}}$ separates $\mathbb{R}_{+}^{d}$ in two open sets. As shown
in [12] (Corollary 1 of Theorem IV 4), and by the connectedness of
$\mathbb{R}_{+}^{d}$, this is possible only if
$\mathrm{dim}\left(B_{\bar{x}}\right)=d-1$. ∎

Before stating our theorem, let us introduce the functions $\Gamma$ and $G$,
defined for any $x\in\mathbb{R}_{+}^{d}$ by
$\Gamma(x)=\left(\frac{x_{1}}{f_{1}(x_{d})},\frac{x_{2}}{f_{2}(x_{1})},...,\frac{x_{d}}{f_{d}(x_{d-1})}\right)$
and $G(x)=\prod\limits_{i=1}^{d}{\frac{x_{i}f^{\prime}(x_{i})}{f(x_{i})}}$,
and the sets $E_{\mathrm{bis}}=\left\\{x\in\mathbb{R}_{+}^{d}:G(x)>1\right\\}$
and if $d\geq 5$,
$E_{\textrm{per}}=\left\\{x\in\mathbb{R}_{+}^{d}:G(x)>\frac{1}{\cos\left(\frac{2\pi}{d}\right)^{d}},G(x)\notin
S_{d}\right\\}$. Note that $E_{\textrm{per}}\subset E_{\textrm{bis}}$.

######  Theorem 2.

Let us assume that $n$ is even.

  1. (i)

If $\alpha\in\overline{\Gamma(E_{\mathrm{bis}})}^{C}$, then (1) has a unique
equilibrium point which is globally asymptotically stable.

  2. (ii)

If $\alpha\in\Gamma(E_{\mathrm{bis}})$, then (1) has exactly three equilibrium
points, among which two are asymptotically stable and one is asymptotically
unstable. Moreover, the union of the basins of attraction of the two stable
equilibria is a dense open subset of $\mathbb{R}_{+}^{d}$.

  3. (iii)

If $d\geq 5$, and if $\alpha\in\Gamma(E_{\mathrm{per}})$, then there exist
periodic solutions of (1).

###### Proof.

First, let us note that, according to the definitions of $\Gamma$ and $G$,
$\bar{x}\in\mathbb{R}_{+}^{d}$ is a fixed point of system (1) if and only if
$\alpha=\Gamma(\bar{x})$, and that for any equilibrium point $\bar{x}$ of (1),

$p_{\bar{x}}^{\alpha}=p_{\bar{x}}^{\Gamma(\bar{x})}=G(\bar{x}).$

  1. (i)

Let us assume that $\alpha\in\overline{\Gamma(E_{\mathrm{bis}})}^{C}$, and let
$\bar{x}$ be a fixed point of (1). Since
$\Gamma(\overline{E_{\mathrm{bis}}})\subset\overline{\Gamma(E_{\mathrm{bis}})}$,
and $\alpha=\Gamma(\bar{x})$,
$\bar{x}\in\overline{\Gamma(E_{\mathrm{bis}})}^{C}$, which means, by
definition of $E_{\mathrm{bis}}$, that $p^{\alpha}_{\bar{x}}=G(\bar{x})<1$.
Since this equality holds for any equilibrium point, we conclude by Lemma 1,
that system (1) has a unique equilibrium point, which is asymptotically
stable, by Lemma 2. Since (1) is an irreducible type K monotone system, this
unique equilibrium is in fact globally asymptotically stable [18].

  2. (ii)

Let us assume that $\alpha\in\Gamma(E_{\textrm{bis}})$. Then, there exists
$\bar{x}\in E_{\textrm{bis}}$ such that $\alpha=\Gamma(\bar{x})$. The point
$\bar{x}$ is thus an equilibrium point of (1) which satisfies
$p_{\bar{x}}^{\alpha}=G(\bar{x})>1$. Therefore, Lemmas 1, 2 and the result of
[18] mentioned at the beginning of the section yield the result.

  3. (iii)

Let us assume that $\alpha\in\Gamma(E_{\textrm{per}})$. Since
$E_{\textrm{per}}\subset E_{\textrm{bis}}$, (1) has exactly three equilibrium
points, denoted $\bar{x}$, $\bar{y}$ and $\bar{z}$, with $\bar{x}\in
E_{\textrm{per}}$, and $\bar{y}$ and $\bar{z}$ which are asymptotically
stable. By definition of $E_{\textrm{per}}$,
$\textrm{dim}\left(B_{\bar{x}}\right)\leq d-3$. By the contrapositive of Lemma
3, system (1) has some divergent solutions, which thus converge to a periodic
orbit, according to [13].

∎

######  Remark.

As mentioned above, this theorem is a more precise version of Theorem 1: we
find the statement of the latter by defining
$A_{\textrm{bis}}=\Gamma\left(E_{\textrm{bis}}\right)$ and
$A_{\textrm{per}}=\Gamma\left(E_{\textrm{per}}\right)$, and by noting that
$E_{\textrm{bis}}$ (resp. $E_{\textrm{per}}$) is empty if and only if $D\leq
1$ (resp. $D\leq\frac{1}{\cos\left(\frac{2\pi}{d}\right)^{d}}$).

##  4 Odd number of decreasing functions

We now deal with the case where $n$ is odd. This case is simpler, since the
system has a unique equilibrium point under this hypothesis. Nevertheless, we
make weaker conclusions regarding the global behaviour of the system, since it
is not an irreducible type K monotone system (see [18]). Thus, we simply study
the linearised system at the neighbourhood of the equilibrium point, and we
conclude with [13], which guarantees that the solutions either converge to
this equilibrium point, or to a periodic orbit.

Let us denote
$E_{unst}:=\left\\{x\in\mathbb{R}^{d}:G(x)<-\frac{1}{\cos\left(\frac{\pi}{d}\right)^{d}}\right\\}.$

We get the following result:

######  Theorem 3.

Let us assume that $n$ is odd, and that $d\geq 3$. Then, system (1) has a
unique equilibrium point. Moreover,

  1. (i)

If $\alpha\in\overline{\Gamma\left(E_{\textrm{unst}}\right)}^{C}$, then this
equilibrium point is asymptotically stable. Moreover, all the solutions of (1)
either converge to this point or to a periodic orbit.

  2. (ii)

If $\alpha\in\Gamma(E_{\textrm{unst}})$, then this equilibrium point is
asymptotically unstable, Moreover, the set of initial conditions for which the
solution converges to a periodic solution is a dense open subset of
$(\mathbb{R}_{+}^{*})^{d}$, and its complement, which is the set of initial
conditions for which the solution converges to the equilibrium point, has
Lebesgue measure zero.

###### Proof.

We use the same notations as in the previous section. Since $n$ is odd,
$\tilde{f}$ is decreasing with $\tilde{f}(0)>0$, it has a unique fixed point,
which implies that system (1) has a unique equilibrium point, that we denote
$\bar{x}$. This point is asymptotically stable if all the eigenvalues of
$M_{\bar{x}}^{\alpha}$ have a negative real part, and asymptotically unstable
if at least one of these eigenvalues has a positive real part. Since $n$ is
odd,

$\mathrm{Sp}\left(M_{\bar{x}}^{\alpha}\right)=\left\\{(\lvert
p_{x}^{\alpha}\rvert)^{1/d}e^{(2k+1)\pi i/d}-1,k\in\\{0,...,d-1\\}\right\\},$

which implies that $\bar{x}$ is stable if $\lvert
p_{\bar{x}}^{\alpha}\rvert<\frac{1}{\cos\left(\frac{\pi}{d}\right)^{d}}$, and
unstable if $\lvert
p_{\bar{x}}^{\alpha}\rvert>\frac{1}{\cos\left(\frac{\pi}{d}\right)^{d}}$. We
conclude by noting that $\lvert
p_{\bar{x}}^{\alpha}\rvert=-p_{\bar{x}}^{\alpha}=-p_{\bar{x}}^{\Gamma(\bar{x})}=-G(\bar{x})$,
and by applying the main theorem of [13] for the first point, and Theorem 4.3
of this same article for the second one.

∎

We recover the result of Theorem 1 by defining
$A_{\textrm{unst}}=\Gamma(E_{\textrm{unst}})$.

## Acknowledgements

The author thanks Nastassia Pouradier Duteil and Camille Pouchol for their
proofreading and their guidance throughout the writing of this paper.

## References

  * [1] David Angeli, James E Ferrell Jr, and Eduardo D Sontag.  Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems.  Proceedings of the National Academy of Sciences, 101(7):1822–1827, 2004.
  * [2] HT Banks and JM Mahaffy.  Stability of cyclic gene models for systems involving repression.  Journal of Theoretical Biology, 74(2):323–334, 1978.
  * [3] Olguta Buse, Rodrigo Pérez, and Alexey Kuznetsov.  Dynamical properties of the repressilator model.  Physical Review E, 81(6):066206, 2010.
  * [4] Joshua L Cherry and Frederick R Adler.  How to make a biological switch.  Journal of theoretical biology, 203(2):117–133, 2000.
  * [5] Timothy S Gardner, Charles R Cantor, and James J Collins.  Construction of a genetic toggle switch in escherichia coli.  Nature, 403(6767):339–342, 2000.
  * [6] Brian C Goodwin.  Oscillatory behavior in enzymatic control processes.  Advances in enzyme regulation, 3:425–437, 1965.
  * [7] Brian C Goodwin et al.  Temporal organization in cells. a dynamic theory of cellular control processes.  Temporal organization in cells. A dynamic theory of cellular control processes., 1963.
  * [8] Jules Guilberteau, Camille Pouchol, and Nastassia Pouradier Duteil.  Monostability and bistability of biological switches.  Journal of Mathematical Biology, 83(6-7):65, 2021.
  * [9] Morris W Hirsch.  Systems of differential equations which are competitive or cooperative: I. limit sets.  SIAM Journal on Mathematical Analysis, 13(2):167–179, 1982.
  * [10] Morris W Hirsch.  Differential equations and convergence almost everywhere in strongly monotone semiflows.  Contemp. Math, 17:267–285, 1983.
  * [11] Morris W Hirsch.  The dynamical systems approach to differential equations.  Bulletin of the American mathematical society, 11(1):1–64, 1984\.
  * [12] Witold Hurewicz and Henry Wallman.  Dimension Theory (PMS-4), Volume 4, volume 63.  Princeton university press, 2015.
  * [13] John Mallet-Paret and Hal Smith.  The poincaré-bendixson theorem for monotone cyclic feedback systems.  Journal of Dynamics and Differential Equations, 2(4):367–421, 1990\.
  * [14] Stefan Müller, Josef Hofbauer, Lukas Endler, Christoph Flamm, Stefanie Widder, and Peter Schuster.  A generalized model of the repressilator.  Journal of mathematical biology, 53:905–937, 2006.
  * [15] Lawrence Perko.  Differential equations and dynamical systems, volume 7.  Springer Science & Business Media, 2013.
  * [16] James F Selgrade.  Asymptotic behavior of solutions to single loop positive feedback systems.  Journal of Differential Equations, 38(1):80–103, 1980.
  * [17] Hal Smith.  Oscillations and multiple steady states in a cyclic gene model with repression.  Journal of mathematical biology, 25(2):169–190, 1987.
  * [18] Hal L Smith.  Systems of ordinary differential equations which generate an order preserving flow. a survey of results.  SIAM review, 30(1):87–113, 1988.
  * [19] Stefanie Widder, Josef Schicho, and Peter Schuster.  Dynamic patterns of gene regulation i: simple two-gene systems.  Journal of theoretical biology, 246(3):395–419, 2007.