Amelia Carolina Sparavigna

Exercises in Archaeoastronomy -2 -The passage mound of Newgrange

Here it is given the second article of a series proposing exercises in archaeoastronomy. The reader can find exercises about the apparent motion of the sun, and its azimuth and altitude given in the horizontal coordinate system. As case study, we investigate the orientation of the prehistoric passage mound of Newgrange.

In the previous article [START_REF] Sparavigna | Exercises in Archaeoastronomy -I -Introduction[END_REF], we have introduced a series of exercises in archaeoastronomy. Let us start the series by studying the sun and its horizontal coordinates. After exercises on azimuth and altitude of the sun, we consider some prehistoric monuments. In particular we propose exercises on the passage mound of Newgrange, as a case study.

First, some definitions. Let us consider the horizon. In astronomy, and therefore in archaeoastronomy, we have to distinguish two horizons. One is the astronomical horizon, which is the horizon that would be seen if the earth's surface were perfectly smooth. It is given by the intersection with the celestial sphere of the local horizontal plane, passing through the observer. The second horizon is the natural or sensible horizon, that is, the line at which the sky and Earth appear to meet. If we use a software simulating azimuth and altitude of the sun or of the moon, such as Sollumis.com, SunCalc.net, SunCalc.org, MoonCalc.org, we have to remember that the lines representing azimuths are given on an astronomical horizon. Therefore, it is possible that, locally, the sunrise is not visible with that specific azimuth, due to the presence of hills and mountains. The same for the visibility of the stars; it is depending on the local natural horizon. As a consequence, when we are considering the astronomical orientations of architectural complexes by means of the above mentioned software, we have to remember that we are doing a preliminary analysis, which is using the local astronomical horizon. This approach allows to evidence the possible existence of alignments along the rise or setting of celestial bodies. After this study, we can make a further analysis by means of the local natural horizon. For this next analysis, we can use Google Earth, and its tool giving the elevation profile (or its tool giving the sun in the local natural horizon [START_REF] Sparavigna | Exercises in Archaeoastronomy -I -Introduction[END_REF]). Of course, if possible, a direct survey of the site will be the final goal of our work.

Before continuing, let us observe the following. We could have an ideal orientation or a natural orientation of the considered site. In the first case, the orientation could be just symbolic, not referring to the actual point of the natural horizon where the sun, the moon or some stars are rising, because the orientation is imagined in an ideal horizon. In the second case, it is prevailing the natural point of the horizon where the celestial object is rising or setting. And probably this point is related to a landmark, such as a peak or a hill, of the local landscape. Therefore, we have always to consider if the horizon had been symbolic or natural for the people that built the monument.

Being the horizon an object which is a little bit more complex of what we have previously told, I suggest to read the site en.wikipedia.org/ wiki/Horizon, which is also explaining the role of the Earth's curvature and atmosphere in the appearance of the celestial objects. For sunrise and sunset, see also wiki/Sunrise and wiki/Sunset.

In one of the software that we will use for the exercises (suncal.org, remarkable software by Torsten Hoffman www.torsten-hoffmann.de) we find azimuth and altitude of the sun for any time of the day, for any day of the year.

Azimuth and altitude are the coordinates of an object in the horizontal frame of reference. This coordinate system uses the observer's local astronomical horizon as the fundamental plane. The altitude is also known as elevation or height. In the Figure 1, we can see the altitude angle and azimuth (please read also wiki/Horizontal_coordinate_system). Let us note that the sunrise azimuth changes on any day of the year. Therefore, the sun is rising on the East only at the equinoxes. And in fact, in the Cambridge Dictionary, we find that "east" is defined as "the direction where the sun rises in the morning that is opposite west" [2], and this happens at equinoxes. From now on, let us remember that we will use the terms East and West (North and South) only for the points of the compass.

First, let us start from sunrise and sunset azimuths, using software sollumis.com. The site tells that the software is giving "The time and direction of sunset, noon, and sunrise for any location". We can click on the map or search any location. "The lines on the drawing show the direction and height (altitude) of the sun throughout the day.

Thicker and shorter lines mean the sun is higher in the sky. Longer and thinner lines mean the sun is closer to the horizon." In a panel on the left (Figure 2), we can read sunrise and sunset azimuths and the noon altitude. Noon altitude is the meridian altitude, used for navigation to calculate an observer's latitude. The meridian altitude is the altitude of the sun (or of another object) when it passes the observer's meridian. Read more at wiki/Meridian_altitude. Sollumis.com is one of the first software that used Google images, showing on them the direction of the sun. Moreover, Sollumis.com was the first software I used for investigations concerning archaeoastronomy [START_REF] Sparavigna | The solar towers of Chankillo[END_REF][START_REF] Sparavigna | The Chinese Pyramids and the Sun[END_REF][5][START_REF] Sparavigna | On the solar orientation of Ales Stenar site[END_REF]. As told in [START_REF] Sparavigna | Exercises in Archaeoastronomy -I -Introduction[END_REF], another software, which is giving the direction of sunrise and sunset on satellite images is suncalc.net. This software was first used for archaeoastronomy in [START_REF] Pankovic | Was Lepenski Vir an ancient Sun or Pleiades observatory?[END_REF], for the analysis of the site of Lepenski Vir.

Exercise 1: Using sollumis.com, find your home and determine for some days of the year, the sunrise and sunset azimuths. For the day when you are doing the exercise, compare the results given by software to the direct observation of the natural horizon. Solution: you can read sunrise and sunset azimuth in the left panel of the software (see Figure 2). Exercise 2: Using the same software, determine the altitude of the sun when it is passing your meridian. Solution:

Again, you can read information in the left panel of software. For instance, on February 12, Torino, the meridian altitude is 31 degrees. Let us note the azimuth in sollumis.com is rounded to the nearest degree. Exercise 5: In en.wikipedia.org/wiki/Newgrange, we read that the site of Newgrange "consists of a large circular mound with an inner stone passageway and chambers. … Its entrance is aligned with the rising sun on the winter solstice, when sunlight shines through a 'roofbox' and floods the inner chamber". Find the sunrise azimuth at Newgrange on December 21. Solution: First, let us consider Newgrange and use sollumis.com. We read a sunrise azimuth of 131 degrees. Then, let us use suncalc.org. To read the sunset azimuth, we have to move the cursor slowly. We have: azimuth 131.29 degrees, altitude 0.06 degrees. See also Figure 5. Read again en.wikipedia.org/wiki/Sunrise. Exercise 6: Consider the sites of Maeshowe in Orkney, Scotland, and Bryn Celli Ddu in Wales, how are their entrances aligned? Solution: use suncalc.net and visit the sites. You can see also [START_REF] Sparavigna | The Summer Solstice at Bryn Celli Ddu[END_REF] for instance. For pictures and maps, search Bryn Celli Ddu and the other site by www.megalithic.co.uk/search.php.

When the above mentioned suncalc.org/.net are used, be careful if you move towards the Arctic Circle. Consider the simulations only after a careful check of the results. The Arctic Circle marks the northernmost point at which the noon sun is just visible on the December solstice. The Circle marks also the southernmost point at which the midnight sun is just visible on the June solstice. Above the circle, we have the Arctic: as seen from this zone of the Earth, the Sun is above the horizon for 24 continuous hours at least once per year and below the horizon for 24 continuous hours at least once per year (and therefore not visible at noon). This is also true in the Antarctic region, south of the Antarctic Circle. Let us note that the position of the Arctic Circle is not fixed, because its latitude depends on the Earth's axial tilt, which fluctuates within a margin of 2 degrees over a 40,000-year period [START_REF] Berger | Obliquity and Precession for the Last 5000000 Years[END_REF]. "By definition, -wiki/Circle_of_latitude explains -the positions of the Tropic of Cancer, Tropic of Capricorn, Arctic Circle and Antarctic Circle all depend on the tilt of the Earth's axis relative to the plane of its orbit around the sun (the "obliquity of the ecliptic"). If the Earth were "upright" (its axis at right angles to the orbital plane) there would be no Arctic, Antarctic, or Tropical circles: at the poles the sun would always circle along the horizon, and at the equator the sun would always rise due east, pass directly overhead, and set due west". Exercise 7. Consider the Figure 6, and imagine the motion of the sun in the four proposed cases.

Figure 6: The motion of the sun as seen by an observer at the centre of the sphere, for the given latitudes.

So we have seen in the satellite images the passage mound of Newgrange (wiki/Newgrange). Newgrange, or Brú na Bóinne, was built during the Neolithic period, around 3200 BC, "making it older than Stonehenge and the Egyptian pyramids" as told by Wikipedia. The site is a large circular mound (diameter of 76 metres), having an inner stone passageway and chambers. It is not only a passage tomb in a mound, it is a monument where the mound possesses a retaining wall with a monumental entrance made by decorated stones.

Wikipedia tells that "There is no agreement about what the site was used for, but it is believed that it had religious significance". Of the link to the winter solstice we have already told in Exercise 5, and of the 'roofbox', through which the light of the sun is flooding the inner chamber during the sunrise on that day (more details at [START_REF] Ray | The winter solstice phenomenon at Newgrange, Ireland: accident or design?[END_REF]9]).

Wikipedia notes that "Several other passage tombs in Ireland are aligned with solstices and equinoxes, and Cairn G at Carrowkeel has a similar 'roofbox'. Newgrange also shares many similarities with other Neolithic constructions in Western Europe, such as Maeshowe in Orkney, Scotland, and Bryn Celli Ddu in Wales".

Wikipedia continues "Many archaeologists believed that the monument [Newgrange] had religious significance of some sort", linked to a "cult of the dead" or to an "astronomically-based faith. The archaeologist Michael J. O'Kelly, who led the 1962-1975 excavations at the site, … believed that Newgrange, alongside the hundreds of other passage tombs built in Ireland during the Neolithic, showed evidence for a religion which venerated the dead as one of its core principles. He believed that this "cult of the dead" was just one particular form of European Neolithic religion, and that other megalithic monuments displayed evidence for different religious beliefs which were solar, rather than death-oriented" (Wikipedia and [START_REF] Grogan | Prehistoric and Early Historic Cultural Change at Brugh na Bóinne[END_REF]). In fact, other scholars have offered an alternative interpretations, in which Newgrange was a monument designed for the capture of the sun on the Winter Solstice. "Once a year, the rising sun of the solstice shines along the long passage of the mound, illuminating the inner chamber. … Michael J. O'Kelly was the first person in modern times to observe this event on 21 December 1967" (Wikipedia). O'Kelly was the person that led the most extensive excavation of the site and also reconstructed the frontage of the mound. Actually, as Wikipedia explains, "after its initial use, Newgrange was sealed for several millennia".

However, it persisted in the Irish mythology and folklore, as a dwelling place of the gods Dagda and his son Aengus. The first antiquarian studies began during the 17th century, and some excavations took place in the years that followed.

As explained in [14], over the course of the excavation, "some of the many local visitors would often tell the O'Kelly's of a tradition, that the rising sun, at some unspecified time, would light up the triple spiral stone in the end recess of the chamber at Newgrange. Unfortunately, no one could be found who had witnessed this but it continued to be mentioned, … A conversation between the O'Kellys (Michael J. and his wife Claire) on the persistence of this tradition, planted in the mind of the Professor that," an orientation of the passageway towards the winter solstice was possible and that the local tradition was based on a real fact. Therefore, O'Kelly decided to reach Newgrange on the day before the winter solstice. "Some minutes before sunrise on the 21st of December 1967", O'Kelly was alone in the darkness of the inner chamber of the mound, waiting and considering if anything would happen. "To his amazement, minute by minute, the chamber grew steadily lighter and a beam of sunlight began to enter the passage and to travel inwards, 'lighting up everything as it came until the whole chamber -side recesses, floor and roof six metres above the floor -were all clearly illuminated'. O'Kelly stood rigid for a while, transfixed by the phenomenon and convinced and fearful in his own imagination that the Dagda, the sun god, who according to the ancient tradition had built the tomb, was about to hurl the roof upon him" [14,[START_REF] O'kelly | Newgrange: archaeology, art and legend[END_REF]. The spectacular phenomenon observed by O'Kelly is slightly different from the original one. "When Newgrange was built over 5000 years ago, the winter solstice sunbeam would have made its way to the back recess of the central chamber. Due to changes in the tilt of the Earth's axis the sunbeam now stops 2 metres from the back recess" (www.knowth.com/winter-solstice.htm). And also, "At dawn, on solstice after 9 am the Sun begin to rise across the Boyne Valley from Newgrange over a hill known as Red Mountain. For the following seventeen minutes, between 19 and 23 December, the sunbeam stretches into the narrow passage of Newgrange tomb and on into the central camber. In Neolithic times it illuminated the rear stone of the central recess of the chamber" hill 93 meters, at a distance of 2700 meters. Since the mound is 12 meters high, the entrance has an elevation of 52 meters. The radius of the mound is of about 38 meters. To see half of the sun, it is necessary that the centre of the sun has an altitude of about a degree ( calculus: arctan ( (93-52) / 2662 ) = 0,88 degree). Let us note that the sun has an angular diameter of 0.5 degrees.

Here, let us consider this value of 1 degree as the value of the altitude of the sun, to have its rays reaching the place in the inner chamber where there is the stone basin, shown in [START_REF]Newgrange and the Winter Solstice[END_REF]. If we consider the figure in this reference, we can prepare a plan of the passage as we are giving in the Figure 8. Let us suppose to consider the winter solstice of 2017, and imagine the sun when it has the same azimuth as it had at the time when Newgrange was built. Today, we know that the rays of the sun stop two meters before the back recess (and therefore let us imagine they stop two meters before the basin). It means that the sun (today) is about 0.7 degrees higher, than it was in origin. This happens because of a very slow change of the tilt of the Earth's axis. Considering a difference of 0.7 degrees, we can estimate the following: when Newgrange was built the sunrise azimuth on winter solstice was larger of an angle comprised between 1.0 and 1.5 degrees (we will discuss in more detail the role of the tilt of the Earth's axis in a future article, during a discussion of the Karnak temple in Egypt). Let us do a last simulation, using CalSky (CalSky: #1 Astronomical and Space Calendar, https://www.calsky.com/, The most complete astronomical observation and information online-calculator on this globe: make your own calendar with satellites, asteroids, comets, planets, sun, and moon). Select Sun, the Rise and Set. To setup the location, use the globe on the right of the page and search Newgrange. After select date and period. For instance, 21 December 2018. The result is given in the Table I. In the web page, we find a glossary too. "Rise and set times are for a mathematical horizon, corrected for standard refraction … ". Read the glossary carefully.

Table I From the Table, we have that the sunrise azimuth is 130.7 degrees (time 8h 41.0m) and the sunset azimuth is Table II From these results we could extrapolate that the sunrise azimuth, at the time Newgrange was built, was of about 132 degrees and the noon altitude of 12.3 degrees.

Let us conclude stressing that Newgrange is a remarkable monument. Moreover, as told in [START_REF] Ray | The winter solstice phenomenon at Newgrange, Ireland: accident or design?[END_REF]9], "Newgrange predates the astronomical structures of Stonehenge by 1,000 years and as such may be the oldest astronomically orientated structure in the world". Other references about Newgrange are given in [START_REF] Patrick | Midwinter sunrise at Newgrange[END_REF][START_REF] Stout | Newgrange and the Bend of the Boyne[END_REF][START_REF] Mitchell | Notes on some non-local cobbles at the entrances to the Passage-Graves at Newgrange and Knowth, county Meath[END_REF]

Figure 1 :

 1 Figure 1: Horizontal frame of reference (adapted from an image by Joshua Cesa).

Figure 2 :

 2 Figure 2: Thanks to sollumis.com, we can have sunrise and sunset azimuths and the noon altitude.

Figure 3 :Figure 4 :

 34 Figure 3: Thanks to suncalc.org, we can have more information concerning the sun.

Figure 5 :

 5 Figure 5: Newgrange seen in SunCalc.org. Note that the pivot is positioned so that the line of the sunrise in passing through the entrance (azimuth 131.69 and altitude 0.25 degrees). Of course, we cannot see the passageway inside, and therefore we can only conclude a possible alignment. As told in [8,9], "As early as 1909, Lockyer remarked that the passage grave was approximately aligned to the rising Sun at midwinter. He did not, however, pay much attention to the site." Let us remember the book [10] by Lockyer.

[

  16]. So, let us investigate the local elevation profile. Exercise 8. Using Google Earth, determine the elevation profile along the sunrise azimuth on solstice, in front of the Newgrange monument. Solution: It is given in the Figure 7.

Figure 7 :

 7 Figure 7: Thanks to Google Earth, we can use a straight edge to draw the red line in the image. The software is

Figure 8 :

 8 Figure 8: This image is showing a sketch of the passage inside Newgrange. The figure gives two rays passing through the roofbox, one imagined today (1.7 degrees), the other at the time when Newgrange was built (1.0 degrees). The difference is due to the very slowly change of the Earth's axis.

229. 3

 3 degrees (time 16h 06.9m). Altitude is 12.9° in Sgr. Day: 7h 25.9m long. If we use a date of 3200 BC, we see that no results are given for sunrise. We see results given from 1500 BC. The TableIIshows them. The southernmost sunrise is at azimuth 131.6 degrees (time 8h 28.6m) and the sunset azimuth is 228.4 degrees (15 h 47.7 m). Altitude 12.5° in Cap. Day: 7h 19.1m long.

Information about this Article

Published on Sunday 18th February, 2018 at 12:04:42.

The full citation for this Article is:

Sparavigna, A. C. (2018). Exercises in Archaeoastronomy -2 -The passage mound of Newgrange. PHILICA Article number 1245.


Amelia Carolina Sparavigna

The Zenith Passage of the Sun in the Plan of Brasilia

In this article we show that the Monumental Axis of Brasilia is oriented along the sunrise direction on the days during which it happens the zenith passage of the sun. This fact could be a mere coincidence. However, since the plan of the town is considered as highly symbolic, this alignment could have been intended as a good omen for the town, to reach the zenith of prosperity and power.

Brasilia is the federal capital of Brazil, founded on April 21, 1960. The capital was planned by Lucio Costa, Oscar Niemeyer designed its civic buildings, and Roberto Burle Marx was the landscape architect. In 1957, an international jury selected the layout proposed by Costa, in spite of the fact that his plan was not as detailed as the plans presented by other architects. It "was chosen by five out of six jurors because it had the features required to align the growth of a capital city" [START_REF] Epstein | Brasilia: Plan and Reality: a study of planned and spontaneous urban development[END_REF]2] The crossed layout of the town is often described as a dragonfly, an airplane or a bird [START_REF] Epstein | Brasilia: Plan and Reality: a study of planned and spontaneous urban development[END_REF]2]. In the framework of an architectural perspective, the airplane-shaped plan was certainly an homage to Le Corbusier [2]. About the symbolism of Costa's design of Brasilia, we find other information in [START_REF] Barnitz | Twentieth-Century Art of Latin America: Revised and Expanded Edition[END_REF]. There, it is told that in the cross of the town, besides the above-mentioned symbols, some scholars have also found "a reference to the cross of the early Portuguese conquerors, the bow and arrow of the native populations of Brazil, or the jet shape as a symbol for future innovation" [START_REF] Sparavigna | The Zenith Passage of the Sun and the Architectures of the Tropical Zone[END_REF]. In the case of the bow and arrow, the arrow is the Monumental Axis and the bow the Residential Axis.

In both the items of Wikipedia, [2] and [START_REF] Barnitz | Twentieth-Century Art of Latin America: Revised and Expanded Edition[END_REF], it is told that the Monumental Axis "points east to west", but this orientation is not that of the cardinal points of the compass. In the satellite images, the Monumental Axis is forming with the cardinal North-South direction an angle of about 108 degrees. 107 degrees is the angle of the azimuth of the sunrise on the days of the zenith passage of the sun, that is, of those days during which the sun reaches the zenith (Figures 3 and4). These days are February 5 and November 5 (Brasilia is between the Equator and the Tropic of Capricorn, and there the sun has two zenith passages). Actually, the Axis is so large that, if we observe the sunrise from its western end, we can see the sun at the other end of the avenue (see the simulation obtained by means of Google Earth in the Figure 5). In the tropical zone, the zenith passage of the sun had influenced the ancient architectures [7]. Therefore, we can consider that also the plan of Brasilia was influenced by the sun and its zenith passage. Of course, the fact that the directions of the Axis and of sunrise are so close could be a mere coincidence; however, since the plan of Brasilia is considered as having a highly symbolic meaning, the alignment of the Monumental Axis could have been intended as a good omen for the future life of the town, that of reaching the zenith of prosperity and power.

  . The initial plan proposed by Costa was transformed over time, but most of its features have survived [2]. In [2], the item of Wikipedia is telling that the foundation of Brasilia, as a "representation of the conquest of an extensive region in Brazil" inspired the symbolism of the plan. In it, Costa used a cross-axial design to indicate "the possession and conquest of this new place" [2,3]. The cross is made by the two principal components of the plan, the Monumental Axis and the Residential Axis (Figures 1 and 2). The Monumental Axis was designated for the political and administrative activities. For the design of the residential apartment blocks, inspiration came form Le Corbusier's Ville Radieuse of 1935 and from the superblocks on the North American Radburn layout, 1929 [4,2].

Figure 1 .

 1 Figure 1. Brasilia as seen from the International Space Station (Courtesy NASA).

Figure 2 .

 2 Figure 2. Brasilia from Google Earth, oriented as in the Pilot Plan drawn by Lucio Costa.

Figure 3 :

 3 Figure 3: Using SunCalc.org, we can determine the zenith passages of the sun. Here it is shown the sunrise direction on February 5 (yellow line). The red line represents the direction of the sunset. The shaded area contains the apparent motion of the sun during the year.

Figure 4 :

 4 Figure 4: A detail of the orientation, from SunCalc.org.

Figure 5 :

 5 Figure 5: Simulation of the sunrise on February 5, observed from the western end of the Axis (TV tower and buildings included).


Aylin Aydoğdu

Sean T Mcquade
email: sean.mcquade@rutgers.edu

Nastassia Pouradier Duteil
email: nastassia.pouradierduteil@rutgers.edu

Opinion Dynamics on a General Compact Riemannian Manifold

This work formulates the problem of defining a model for opinion dynamics on a general compact Riemannian manifold. Two approaches to modeling opinions on a manifold are explored. The first defines the distance between two points using the projection in the ambient Euclidean space. The second approach defines the distance as the length of the geodesic between two agents. Our analysis focuses on features such as equilibria, the long term behavior, and the energy of the system, as well as the interactions between agents that lead to these features. Simulations for specific manifolds, S 1 , S 2 , and T 2 , accompany the analysis. Trajectories given by opinion dynamics may resemble n-body Choreography and are called "social choreography". Conditions leading to various types of social choreography are investigated in R 2 .

Introduction

The emergence of a group's global behavior from local interactions among individual agents is a fascinating feature of opinion dynamics. When local rules imply global patterns in a population, we are observing a phenomenon called self-organization. Traditionally, interest focuses on understanding the complex rules of interacting opinions which lead to certain global configurations, such as classic consensus, alignment, clustering, or the less studied dancing equilibrium [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF]. For instance, in bounded-confidence models such as the one proposed by Hegselmann and Krause, the radius of interaction determines the clustering of the system [START_REF] Hegselmann | Opinion dynamics and bounded confidence: models, analysis and simulation[END_REF]. Motsch and Tadmor studied the influence of the shape of the interaction potential on the convergence to consensus of the Hegselmann-Krause system [START_REF] Motsch | Heterophilious dynamics enhances consensus[END_REF]. Ha, Ha and Kim looked at the Cucker-Smale second-order alignment model and provided a condition on the interaction potential ensuring convergence of the system to alignment [START_REF] Ha | Emergent behavior of a Cucker-Smale type particle model with nonlinear velocity couplings[END_REF]. Cristiani, Frasca and Piccoli studied the effect of anisotropic interactions on the behavior of the group [START_REF] Cristiani | Effects of anisotropic interactions on the structure of animal groups[END_REF].

The dynamics of an opinion formation system depend on the state-space [START_REF] Aydoğdu | Interaction network, state space and control in social dynamics[END_REF] and interaction network [START_REF] Hegselmann | Opinion dynamics and bounded confidence: models, analysis and simulation[END_REF]. Models on the Euclidean space in one dimension (for opinion dynamics) or in two or 1 three dimensions (with applications to groups of animals or robots) have been extensively studied and are well understood. However, such models are locally linear, which may be a limitation when one strives to capture more complex phenomena and better represent reality [START_REF] Sobkowicz | Modelling opinion formation with physics tools: Call for closer link with reality[END_REF]. In this line of thought, the Kuramoto model on the sphere S 1 addresses the problem of synchronizing a large number of oscillators [START_REF] Kuramoto | Cooperative dynamics of oscillator community a study based on lattice of rings[END_REF][START_REF] Strogatz | From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators[END_REF]. There exist numerous applications to this model [START_REF] Dörfler | Synchronization in complex oscillator networks and smart grids[END_REF][START_REF] Scardovi | Synchronization and balancing on the N-torus[END_REF][START_REF] Sepulchre | Stabilization of planar collective motion: All-to-all communication[END_REF][START_REF] Sepulchre | Stabilization of planar collective motion with limited communication[END_REF]. Similarly, applications to satellite or ground vehicle coordination have motivated the development of models on special orthogonal groups [START_REF] Sarlette | Coordinated motion design on lie groups. Automatic Control[END_REF][START_REF] Sarlette | Consensus optimization on manifolds[END_REF]: satellites evolve on SO(3) while ground vehicles evolve on SE [START_REF] Bliss | The geodesic lines on the anchor ring[END_REF] or SE [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF]. A nonlinear model of opinion formation on the sphere was also developed in [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF].

The present work defines a general model of opinion dynamics on a Riemannian manifold. We investigate how the manifold on which the model is defined affects the global configurations resulting from opinion dynamics. These are the first steps to build a robust theory of opinion dynamics on general Riemannian manifolds.

There is an inherent difficulty in defining opinion dynamics on a general Riemannian manifold. Using the Riemannian distance, an agent will move towards a point by following the manifold's geodesics, which are well defined only locally. On a larger scale, there might not exist a unique geodesic. Another challenge is the extreme complexity of computing geodesics, even on a relatively simple manifold such as the torus [START_REF] Gravesen | The cut locus of a torus of revolution[END_REF]. One way around this issue is to consider the embedding of the manifold into a Euclidean space. Each agent's velocity is defined by projection of the other agents' influence onto the tangent space at that point. This is the choice made in [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF].

Other than the mentioned practical aspect, there is an intrinsic rationale for choosing one approach over the other. When evolving along the geodesics of the manifold, one assumes that each agent has a global understanding of the manifold's geometry and is able to choose the shortest path among all possible ones. On the other hand, the approach based on the projection of the desired destination onto the tangent space implies that each agent only holds local information about the space in which it evolves. It chooses to move in the direction which locally seems to bring it closer to the target.

We explore these two specific approaches for our generalized model. The first method, Approach A, uses projections in the Euclidean space in which the manifold is embedded. The second method, Approach B, uses only geodesics defined on the manifold to define strength and direction of interaction. We exhibit properties of the interaction matrix that lead to specific kinds of equilibria.

Simulations and examples compare the two methods. Dancing equilibria for Approach B are shown (dancing equilibria were studied for Approach A in [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF]).

We use the sphere and torus as sample manifolds to evaluate these approaches. Specifically, we simulate dynamics on the following manifolds: S 1 , S 2 and T 2 . These examples allow us to directly compare the two approaches, and see if one is more appropriate for a given manifold. We show the influence of the manifold's geometry on the dynamics by examining the dynamics resulting from the same interaction matrix in S 2 , T 2 and R 2 .

Opinion dynamics trajectories can resemble n-body choreography, that is, solutions to the well known n-body problem. These dynamics drive agents along orbits which either are periodic, or have a periodic feature, and that may be shared by multiple agents. We refer to opinion dynamics trajectories along such orbits as "Social Choreography". We show that a simple example of Social Choreography in R 2 does not hold on S 2 or T 2 , see Figures 10 and11. We then focus on the case of R 2 and investigate initial conditions and properties of the interaction matrix which give rise to Social Choreography. Future investigations can consist of exploring similar properties of trajectories on general manifolds.

Choice of the model

This work will primarily discuss two approaches to define opinion dynamics on a Riemannian manifold. Let M be a Riemannian manifold. Let N ∈ N represent the number of agents with opinions evolving on M . We denote by x := (x i ) i∈{1,...,N } ∈ M N the set of opinions. For each i ∈ {1, ..., N }, ẋi ∈ T xi M . The opinions x i evolve according to the following general dynamics:

ẋi = N j=1 a ij Ψ(d(x i , x j ))ν ij (1)
where

• a ij ∈ R is the interaction coefficient of the pair of agents i and j,

• Ψ : R → R is the interaction potential,

• d(•, •) : M × M → R + represents the distance between opinions,
• ν ij ∈ T xi M is a unit vector giving the direction of the influence of j over i.

Each of these terms is further specified in the following.

Approaches

The evolution of each agent's opinion depends on the opinions of all other agents, with influences weighted by the interaction coefficients a ij . More specifically, an agent x j 's influence on x i is determined by two elements: the direction of influence ν ij ∈ T xi M and the magnitude of influence

Ψ(d(x i , x j )) ∈ R + .
We propose and study two different approaches for the choices of d and ν ij . Approach A uses the embedding of M in R n to define d(x i , x j ), whereas Approach B is intrinsic to M , with distance and direction of influence based on geodesics.

Approach A. Assume that M of dimension m is embedded in a Euclidean space R n , with n ≥ m. Agent x j acts on agent x i via a projection onto T xi M ⊂ R m . Now considering points (x i , x j ) ∈ M 2 as points of R n , the difference x j -x i is a vector of R n . Given a vector subspace Y of R n , we denote by Π Y y the projection of y ∈ R n onto Y ⊂ R n and define d P (•, •) as follows:

d P (x i , x j ) = Π Tx i M (x j -x i ) (2)
where • denotes the Euclidean norm on R n . The same projection also defines the direction of influence of x j on x i :

ν ij =    Π Tx i M (xj -xi) Π Tx i M (xj -xi) if Π Tx i M (x j -x i ) = 0 0 otherwise. (3)
With the specific choice Ψ ≡ Id, system (1)-( 2)-(3) becomes:

ẋi = N j=1 a ij Π Tx i M (x j -x i ). (4)
x j -x i x j

x i Figure 1: An example of a manifold M such that d p (x i , x j ) = d p (x j , x i ), Using system (4), an agent is subject to "local visibility", and movement of x i along T xi M (dashed line through x i ) will not bring x i closer to x j in this local sense. This is the approach used in [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF], applied to the sphere S 2 . Notice that the magnitude of influence, d P (x i , x j ), is symmetric for the sphere in the sense that d P (x i , x j ) = d P (x j , x i ), but not symmetric for a general Riemannian manifold (see Figure 1). However, it is a continuous function defined for all pairs of points (x i , x j ) ∈ M 2 . The originality of this approach is that the influence of x j on x i is not related to a notion of distance between the points. The use of the projection of x j -x i onto T xi M reflects the concept of "local visibility." For the situation of two agents evolving on a one dimensional manifold, if x j -x i ⊥ T xi M , then a local displacement of x i does not affect the distance between the points x i -x j . Indeed, a first order Taylor expansion gives:

x i (ε) = x i (0) + ε ẋi (0) + o(ε).
Supposing that x j is fixed, we have:

x i (ε) -x j 2 = x i (ε) -x j , x i (ε) -x j = x i (0) -x j , x i (0) -x j + 2ε ẋi (0), x i (0) -x j + o(ε) so if x j -x i (0) ⊥ T xi(0) M , then x i (ε) -x j 2 = x i (0) -x j 2 + o(ε).
Hence if x i only has local visibility, all directions of displacement seem equivalent (at first order), which justifies the influence of x j over x i to be zero if their difference is orthogonal to the tangent space of M at x i . This is illustrated in Figure 1.

Approach B. This second approach defines d and ν ij using the manifold M itself, and does not require any reference to the space in which M is immersed. This would make Approach B a natural way to define system dynamics, however the complete knowledge of the geodesics between any two points on the manifold may be unrealistic. Furthermore, the geometry of the manifold may introduce difficulties to the uniqueness of ν ij , particularly at the cut-locus of a point. Definition 2.1. The cut locus of a point q ∈ M is the set of points CL(q) ⊂ M for which there are multiple geodesics between q and p ∈ CL(q) (see also [START_REF] Cheeger | Comparison theorems in Riemannian geometry[END_REF]).

Let γ ij : [0, 1] → M denote a geodesic connecting x i to x j , γ ij (0) = x i and γ ij (1) = x j .
We then define the distance between x j and x i as the length of a geodesic, i.e. denoting by

g y : T y M × T y M → R + the Riemannian metric at point y ∈ M , d G (x i , x j ) = 1 0 g γij (s) ( γij (s), γij (s))ds. ( 5
)
The direction of influence is determined by the same geodesic:

ν ij = 0 if x j = x i or if x j ∈ CL(x i ) γij (0) √ gx i ( γij (0), γij (0)) otherwise. (6)
Unlike in Approach A, the magnitude of influence is a symmetric function:

d G (x i , x j ) = d G (x j , x i ).
Furthermore, this approach ensures that the magnitude of influence of one agent on another is a function of the exact Riemannian distance between the agents.

Interaction networks. In finite-dimensional systems such as system (1), the set of interacting agents can be described by vertices of a graph. A directed edge exists from a vertex i to a vertex j if and only if a ij = 0. The system depends on the interaction network, and likewise, if the coefficients a ij are chosen to be functions of the state, the interaction network may change as a result of the dynamics. Two main types of interaction networks have been proposed in the literature: metric interactions and topological interactions. If interactions between agents occur only locally, only the neighbors of agent i influence agent i. Metric interactions define the set of neighbors of agent i, given a radius r > 0, as

S r i (x) = {j ∈ {1, . . . , N }, d(x i , x j ) ≤ r}, (7)
where d(•, •) can represent either the projection or the geodesic distance, as specified in each of the two approaches described above (see equations ( 2) and ( 5)). The other main type of interactions specifies that an agent is influenced by only its k closest neighbors. We call these topological interactions [START_REF] Aydoğdu | Interaction network, state space and control in social dynamics[END_REF]. We define the relative separation between two agents as

α ij = card{k : d(x i , x k ) ≤ d(x i , x j )},
The set of neighbors of agent i is then defined as the set of its k closest neighbors, i.e. for a given k ∈ N,

S k i (x) = {j ∈ {1, . . . , N }, α ij ≤ k}. (8)
Figures 2 and 3 illustrate differences between the metric and topological networks for the specific example of S 1 , with each of the approaches A and B.

Tx 1 S x 3 x 4 x 1 x 2 x 5 x 6 (a) Approach A metric, r = π 4 Tx 1 S x 3 x 4 x 1 x 2 x 5 x 6 (b) Approach A topological, k = 2 x 3 x 4 x 1 x 2 x 5 x 6 π 4 (c) Approach B metric, r = π 4 x 3 x 4 x 1 x 2 x 5 x 6 π 4 (d) Approach B topological, k = 2 Figure 2:
The set of agents that influence x 1 depends on how the interaction network is defined. In (a) and (b) the dashed lines show the projection of agents onto the tangent space of x i , (T xi S 1 ).

The agents depicted in red with larger dots influence x 1 . With the same configuration on S 1 , four combinations are possible (approach {A,B} type {Metric, Topological}). Each combination implies x 1 interacts with a different set of agents.

x 1

x 2

x 3

x 4

x 5

x 6

(a) Approach A metric, r = π 4

x 1

x 2

x 3

x 4

x 5

x 6

(b) Approach A topological, k = 2

x 1

x 2

x 3

x 4

x 5

x 6

(c) Approach B metric, r = π 4

x 1

x 2

x 3

x 4

x 5

x 6

(d) Approach B topological, k = 2
Figure 3: The agent x 1 is influenced by different agents depending on how the interaction network is defined. These networks may change as the dynamics move the agents on S 1 . Each agent x j , j ∈ {1, . . . , 6} will have a network describing which other agents influence x j . The interaction networks corresponding to systems from Figure 2.

Resolution of discontinuities. The definitions of ν ij for approaches A and B (given by equations ( 3) and ( 6)) allow discontinuities of ν ij at certain points. Thus, one must impose conditions on the interaction potential Ψ ∈ C 0 (R + , R + ), in order to ensure the continuity of the right-hand side of the system (1), and hence the existence and uniqueness of a solution. Table 1 lists the discontinuities of ν ij and gives necessary conditions on Ψ to ensure the continuity of Ψ(d(x i , x j ))ν ij . Firstly, notice that in both approaches, ν ij is discontinuous at the point x i = x j . Indeed, if x i = x j , ν ij = 0, whereas almost everywhere else, ν ij = 1. To ensure the continuity of Ψ(d(x i , x j ))ν ij at this point, we impose the following condition:

Ψ(0) = 0. (9)
In Approach A, we created a discontinuity of ν ij at the points x j ∈ N (x i ), where we denote by N (q) the set N (q) := {q ∈ M | Π TpM (q -p) = 0}. For convenience of notation, we will use interchangeably the notations N (x i ) and N i . More specifically, we have lim xj →Ni ν ij = 1 but ν ij = 0 if x j ∈ N i (see also Table 1). However, from the definition of d P (see equation ( 2)), we have lim xj →Ni d P (x i , x j ) = 0 and d(x i , x j ) = 0 for x j ∈ N i . Hence a sufficient condition for Ψ(d(x i , x j ))ν ij to be continuous is again:

Ψ(0) = 0. (10)
In Approach B, there is a discontinuity for x j ∈ CL(x i ). Denoting by B geo (p, ρ) the geodesic ball of center p and radius ρ, we require the following condition on the influence function Ψ:

Ψ(d) = 0 for all d ≥ (11)
where

:= inf{ρ > 0 | ∀p ∈ M, B geo (p, ρ) ∩ CL(p) = ∅}.
This distance , also known as injectivity radius, is known to exist and be greater than 0 for any compact Riemannian manifold (see [START_REF] Cheeger | Comparison theorems in Riemannian geometry[END_REF]).

Approach A B A and B Critical points

x j ∈ N i x j ∈ CL(x i ) x j = x i Discontinuities lim xj →Ni ν ij = 1 lim xj →CL(xi) ν ij = 1 lim xj →xi ν ij = 1 ν ij = 0 for x j ∈ N i ν ij = 0 for x j ∈ CL(x i ) ν ii = 0 Condition on Ψ Ψ(0) = 0 Ψ(d) = 0 for all d ≥ Ψ(0) = 0
Table 1: Possible discontinuities of the right-hand side of (1). The bottom row of the table show conditions for Ψ so that the system is continuous.

Notice that in the case of the geodesics approach (Approach B), the condition Ψ(d) = 0 for all d ≥ is incompatible with the use of the topological network [START_REF] Ha | Emergent behavior of a Cucker-Smale type particle model with nonlinear velocity couplings[END_REF]. Indeed, if agent j is among the k closest neighbors of agent i, the topological network would require: a ij = 0. However, the interaction between i and j would be canceled if d G (x i , x j ) > . On the other hand, the metric interaction network as defined by ( 7) is compatible with Approach A, and with Approach B if the interaction radius is smaller than the injectivity radius: r ≤ . For simplicity purposes, in the rest of this paper, we will consider that the interaction coefficients a ij are constant, thus not requiring the need to differentiate between metric and topological networks. While models with constant interaction coefficients are our focus here, these models are quite restrictive, and exclude all models with dynamic interactions.

Definitions and general results

Definition 2.2. The configuration x 1 = ... = x N is called consensus. On the sphere, S n , A configuration such that, for every j ∈ {2, . . . , N }, either x j = x 1 or x j = -x 1 , which is not a concensus is called antipodal equilibrium.

Proposition 2.1. The consensus configuration is an equilibrium for system (1).

Proof. In both approaches A and B, if x i = x j , then ν ij = 0. Hence if x 1 = ... = x N , then for all i ∈ {1, . . . , N }, ẋi = 0. Proposition 2.2. Let N > d + 1. Then for every x = (x 1 , . . . , xN ) ∈ M N , there exists a square matrix A = (a ij ) i,j∈{1,...,N } such that x is an equilibrium for system (1).

Proof. The configuration x = (x 1 , . . . , xN ) is an equilibrium if and only if

d dt xi = N j=1 a ij Ψ(d(x i , xj ))ν ij = 0. This is a system of at most N d equations in the N 2 -N unknowns a ij , i = j, notice that Ψ(d(x i , x i )ν ii = 0
, and diagonal values of A do not change the system. So if N > d + 1 there exists a nontrivial choice of the interaction coefficients for which x is an equilibrium.

Definition 2.3. The kinetic energy of System (1)-( 2)-( 3) is the quantity

E P (t) := 1 2 N i=1 ẋi (t) 2 . ( 12
)
The kinetic energy of System (1)-( 5)-( 6) is the quantity

E G (t) := 1 2 N i=1 g xi ( ẋi (t), ẋi (t)). (13)
Proposition 2.3. Let M be a general Riemannian bounded manifold. Consider the dynamics given by projection onto the tangent space (Approach A) given by (4). If the interaction matrix A = (a ij ) i,j∈{1,...,N } 2 is symmetric, then

lim t→∞ E P (t) = 0. ( 14
) Proof. Let F (t) = 1 2 N i=1 N j=1 a ij x i -x j 2 .
Using the symmetry of A, we prove that

d dt F (t) = 4E P (t). (15)
Indeed, notice that

∇ xi ( N j=1 a ij x i -x j 2 ) = 2Π Tx i M N j=1 a ij (x j -x i ) = 2 ẋi .
Then we compute

d dt F (t) = N k=1 ∇ x k 1 2 N i,j=1 a ij ( x i -x j 2 ), ẋk = N k=1 ∇ x k   1 2 N i=1 a ik ( x i -x k 2 ) + 1 2 N j=1 a kj ( x k -x j 2 )   , ẋk = N k=1 2∇ x k 1 2 N i=1 a ik ( x i -x k 2 ), ẋk = N k=1 2Π Tx k M N j=1 a kj (x j -x k ), ẋk = 2 N k=1 ẋk 2 = 4E P (t). ( 16
)
where the third equality uses the property: a ij = a ji for all i, j.

Since E P (t) ≥ 0, F (t) is a non-decreasing function. Moreover F (t) and d 2 dt F (t) are bounded, since M is a bounded manifold. Hence d dt F (t) → 0 when t → ∞, which implies that lim t→∞ E P (t) = 0. Remark 2.1. Propositions 2.

and 2.3 are generalizations of results proven for the case

M = S 2 in [3].
Remark 2.2. Proposition 2.3 assumes that Ψ ≡ Id which creates a discontinuity for Approach B (see Table 1). A result is shown for the more restricted case of M = S 2 and Ψ(•) ≡ sin(•), (Corollary 4.1). Definition 2.4. Let x solve the differential equation [START_REF] Aydoğdu | Interaction network, state space and control in social dynamics[END_REF]. A dancing equilibrium is a configuration in which for all pairs of agents (i, j), the distance

d P (x i , x j ) (in Approach A) or d G (x i , x j ) (in Approach B) is constant. Remark 2.
3. This definition is a generalization of the concept of dancing equilibrium described in [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF].

Remark 2.4. It follows immediately from definition 2.4 that the kinetic energy of a system in dancing equilibrium is constant.

3 Analysis and simulations on S 1

Models

We study both approaches A and B in the case M = S 1 , i.e. for the one-dimensional sphere embedded in

R 2 . Let (θ i ) i∈{1,...,N } ∈ [0, 2π] N such that for all i ∈ {1, ..., N }, x i = (cos θ i , sin θ i ) T .
Approach A. The projection onto an agent's tangent space can be rewritten as:

Π Tx i N j=1 a ij (x j -x i ) = N j=1 a ij cos θ j sin θ j - cos θ i sin θ i , -sin θ i cos θ i -sin θ i cos θ i = N j=1 a ij (-sin θ i cos θ j + sin θ j cos θ i ) -sin θ i cos θ i = N j=1 a ij sin(θ j -θ i ) -sin θ i cos θ i . (17)
So System (1)-( 2)-(3) becomes:

for all i ∈ {1, . . . , N }, θi -sin θ i cos θ i = N j=1 a ij Ψ(| sin(θ j -θ i )|) sgn(sin(θ j -θ i )) -sin θ i cos θ i (18)
where sgn(•) is the sign function defined by:

for all x ∈ R, sgn(x) =      1 if x > 0 -1 if x < 0 0 if x = 0. (19)
We can then specify:

for all (i, j) ∈ {1, . . . , N } 2 , d P (x i , x j ) = | sin(θ j -θ i )|, ν P ij = sgn(sin(θ j -θ i )). (20)
This gives the system of scalar equations:

for all i ∈ {1, . . . , N }, θi = N j=1 a ij Ψ(| sin(θ j -θ i )|) sgn(sin(θ j -θ i )). (21)
In particular, in the case Ψ ≡ Id, the system becomes the Kuramoto model [START_REF] Kuramoto | Cooperative dynamics of oscillator community a study based on lattice of rings[END_REF].

for all i ∈ {1, . . . , N }, θi = N j=1 a ij sin(θ j -θ i ). (22)
Approach B. For M = S 1 , the geodesics distance d G and the vector ν G ij are given by:

d G (x i , x j ) = arccos(cos(θ j -θ i )) , ν G ij = sgn(sin(θ j -θ i )). (23)
System ( 1)-( 5)-( 6) is written:

for all i ∈ {1, . . . , N }, θi = N j=1 a ij Ψ(arccos(cos(θ j -θ i ))) sgn(sin(θ j -θ i )). (24)
In order for the system to be well defined, the interaction function Ψ must satisfy the conditions given in Table 1. Notice that the injectivity radius is constant over S 1 , with = π. Possible choices involve choosing Ψ from a family of function defined as follows:

Ψ a (d) = 1 a d for d ≤ a d-π a-π for d > a ( 25
)
where a ∈ (0, π) (see Figure 8). Another possible choice is: Ψ : x → sin(x). Notice that for the specific choices Ψ = Id for Approach A and Ψ : x → sin(x) for Approach B, the two approaches A and B are equivalent.

Analysis

We first examine the different equilibria for both approaches. Theorem 3.1. Consider Approach A, System (21). Let N ∈ N be even. Suppose that for all i ∈ {1, . . . , N } for all j ∈ {1, . . . N 2 }, a ij = a i(j+ N 2 ) . Then any configuration that is centrally symmetric, i.e. for all j ∈ {1, . . . ,

N 2 }, θ j+ N 2 = θ j + π is an equilibrium.
Proof. Using the hypotheses from Theorem 3.1, we can easily compute:

θi = N j=1 a ij Ψ( sin(θ j -θ i ) ) sgn(sin(θ j -θ i )) = N/2 j=1 [a ij Ψ( sin(θ j -θ i ) ) sgn(sin(θ j -θ i )) + a i(j+ N 2 ) Ψ( sin(θ j+ N 2 -θ i ) ) sgn(sin(θ j+ N 2 -θ i ))] = N/2 j=1 [a ij Ψ( sin(θ j -θ i ) ) sgn(sin(θ j -θ i )) + a ij Ψ( sin(θ j + π -θ i ) ) sgn(sin(θ j + π -θ i ))] = 0.
Interestingly, Theorem 3.1 is not applicable to Approach B. We illustrate the different behaviors of the two systems by studying the specific example of four agents initially in a rectangular configuration. According to Theorem 3.1, this configuration is an equilibrium for Approach A, independently of the choice of interaction function Ψ. However, one can easily prove that in the geodesics-based Approach B, with N = 4 and the choice Ψ := Ψ a with a = 3π 4 , the only equilibrium for which all agents have pairwise distinct positions is obtained by a regular polygon, i.e. all agents are evenly spaced out on the circle. This is illustrated by numerical simulations shown in Figure 4.

This highlights the fundamentally different behaviors of the systems ( 1)-( 2)-( 3) and ( 1)-( 5)-( 6) in the case M = S 1 . 25)). Notice that with Approach A, initial and final positions are identical since any rectangle configuration is an equilibrium. However, with Approach B, the system reaches a square configuration, the only possible equilibrium with pairwise distinct positions.

In both approaches A and B, conditions on the interaction matrix A can be found such that the system forms a dancing equilibrium (see Definiton 2.4). Theorem 3.2. Consider the dynamics on S 1 given by:

for all i ∈ {1, . . . , N }, θi = N j=1 a ij Ψ(d(x i , x j ))ν ij (26)
where d(•, •) and ν are given either by Approach A [START_REF] Strogatz | From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators[END_REF] or Approach B (23). Let C ∈ R and suppose that for all i ∈ {1, . . . , N },

a ij = C Ψ(d(xi(0),xj (0))) ν ij if Ψ(d(x i (0), x j (0))) = 0 0 otherwise. ( 27
)
Then the system is in a dancing equilibrium.

Proof. If the interaction matrix satisfies (27), then at t = 0,

for all i ∈ {1, . . . , N }, θi (0) = N j=1 C = CN
so for all (i, j) ∈ {1, . . . , N } 2 , θi (0) -θj (0) = 0. Then d(x i , x j ) does not change in time, and ( 27) holds for all time.

Numerical simulations show the evolution of the system (26) with condition (27) for the projection or the geodesic distance, see Figures 5 and6. 4 Analysis and simulations on S 2

Models

We study both approaches A and B for M = S 2 , i.e. for a two dimensional sphere embedded in R 3 . We use spherical coordinates: let (θ i ) i∈{1,...,N } ∈ [0, 2π] N , and (φ i ) i∈{1,...,N } ∈ [0, π] N such that for all i ∈ {1, ..., N }, x i = (cos θ sin φ, sin θ sin φ, cos φ) T .

Choice of influence function

We choose an influence function Ψ(d) between two agents x i and x j so that the right-hand side of the system is continuous, the discontinuities are shown in Table 1. For Approach B, the only point in CL(x i ) for a given x i is the antipodal point (this is an end point of a diameter for which x i is the other end point.) As in the case of S 1 , for Approach B, we choose a function Ψ from a family of functions of the form Ψ a , see equation (25) (choices of Ψ are shown in Figure 8).

Approach A. On S 2 , the derivative for system (1)-( 2)-( 3) with Ψ ≡ Id reduces to the sum of all projections onto the tangent space of agent x i , weighted by the corresponding interaction term a ij . This is rewritten as:

Π Tx i N j=1 a ij (x j -x i ) = Π Tx i N j=1 a ij (x j ) = N j=1 a ij (x j -x j , x i x i ) = N j=1 a ij   cos θ j sin φ j sin θ j sin φ j cos φ j   -   cos θ j sin φ j sin θ j sin φ j cos φ j   ,   cos θ i sin φ i sin θ i sin φ i cos φ i     cos θ i sin φ i sin θ i sin φ i cos φ i   .
Approach B. The geodesic distance d G (x i , x j ) from (5) between two points x i , and x j on S 2 is given by:

d G (x i , x j ) = 2 arcsin x i -x j 2 ,
and the direction toward x j from x i is

ν ij = x j -x j , x i x i x j -x j , x i x i ,
where • is the standard norm in R 3 . Noticing that for S 2 , Approach A with Ψ ≡ Id is equivalent to Approach B with Ψ ≡ sin, we extend the results of Proposition 2.3.

Corollary 4.1. Consider the dynamics given by geodesic distance (Approach B) on S 2 , system (1)-( 5)-( 6), and let Ψ ≡ sin. If the interaction matrix A = (a ij ) i,j∈{1,...,N } 2 is symmetric, then

lim t→∞ E G (t) = 0. ( 28
)
Proof. System (1)-( 5)-( 6), with Ψ(•) ≡ sin(•) reads:

ẋi = N j=1 a ij sin(d G (x i , x j ))ν ij . ( 29
)
Considering the embedding of the system in R 3 , we notice that for all i, j ∈ {1, . . . , N },

a ij sin(d G (x i , x j ))ν ij = a ij Π Tx i M (x j -x i ), (30)
thus the system is (4), and by proposition 2.3, lim t→∞ E P (t) = 0. Finally,

lim t→∞ E P (t) = 0 =⇒ lim t→∞ ẋi = 0 for all i =⇒ lim t→∞ E G (t) = 0. (31)

Simulations

We use a fourth order Runge-Kutta scheme to approximate the trajectories. The derivative is calculated as a vector in R 3 , and then we express this vector in spherical coordinates, θt θ and φt φ where t θ and t φ are the unit vectors in the direction of the azimuth angle( θ) and the polar angle( φ) respectively for the ith agent. {t θ , t φ } form an orthonormal basis for the tangent space of x i . Using the angular derivatives avoids having to calculate the agent's trajectory in R 3 and then project onto the sphere for every iteration which would cause significant numerical errors.

For an agent x i = (θ i , φ i ) we can write t θi and t φi as

t θ =   -sin θ cos θ 0   , t φ =   cos θ cos φ sin θ cos φ -sin φ   .
We express the derivative of an agent as

ẋi = ∂x i ∂θ i θ + ∂x i ∂φ i φ. ( 32
)
By direct computation, we get:

∂x i ∂θ i =   -sin θ i sin φ i cos θ i sin φ i 0   , ∂x i ∂φ i =   cos θ i cos φ i sin θ i cos φ i -sin φ i   .
We can also express the derivative of x i as the projection of the derivative in R 3 with t θi and t φi ẋi = ẋi , t θi t θi + ẋi , t φi t φi .

It follows from ( 32) and (33) that θi = 1 sin φ i ẋi , t θi t φi and φi = ẋi , t φi t φi .

Singularities: In (34), the factor 1 sin φ causes a singularity around φ = kπ for a non-negative integer k. To avoid this practical problem, before each iteration of the RK4 scheme, we identify critical agents that have a polar angle close to 0 or π (φ = kπ), for non-negative integer k; for these critical agents, we rotate all agents π 2 around the x-axis to calculate the derivative. This is a concern for both Approach A and Approach B.

An additional concern is singularities for agents forming consensus. In equations ( 3) and ( 6), ν ij is normalized, and when agent i and agent j are very close together, dividing by d P and d G causes a singularity. We avoid this problem in our simulations by defining a minimum distance d min between agents for the sake of the normalization term in the denominator. When d(x i , x j ) < d min then

ν ij = x j -x j , x i x i d min

Examples

We ran simulations using different choices of Ψ to see how this choice can impact the system. We show two simulations, the first uses Approach A (Figure 7); the second uses Approach B (Figure 9). In both pictures the same interaction matrix is used (given in the appendix), and we see that our choice of Ψ (Figure 8) may dramatically change the system behavior. The second example shows the effect of curvature on the system (Figure 10) for comparison to T 2 and R 2 (Figure 11). Another example is the appendix shows unexpected behavior using Approach A. Five agents with a general interaction matrix A and Ψ as defined in (25) with a ∈ { π 4 , π 2 , 3π 4 }. The behavior of the system can change dramatically from our choice of Ψ (Figure 8), and dancing equilibria may arise from A with certain configurations, (Figure 9). Example 4.2. To assess the influence of the curvature of S 2 on the dynamics, observe a simple case involving 3 agents evolving according to the interaction matrix:

A =   0 1 -1 -1 0 1 1 -1 0   (35)
In Section 6.2, we prove that those dynamics in R 2 lead to periodic trajectories on a single orbit shared by all three agents, the orbit's parameters being fully determined by the initial conditions (see Theorem 6.2). However, the same dynamics on the sphere do not give rise to periodic trajectories. In sections 5.3, we also discuss the dynamics with this interactions matrix on T 2 , to assess the effect of curvature of the manifold.

Figure 10: Dynamics with Approach A on S 2 , using the interactions matrix (35). If the agents' initial positions are close enough to each other, the agents with will form trajectories that remain in a neighborhood of their initial position.

Analysis and simulations on T 2

We now study how the general dynamics given by equation ( 1) apply to the specific case of the torus T 2 ⊂ R 3 . Let (e x , e y , e z ) denote the Euclidean basis of R 3 . Let (R, r) ∈ (R + ) 2 , with R > r.

We define the manifold T 2 as the torus obtained by rotating the circle (x-R) 2 +z 2 = r 2 around the z-axis. Hence T 2 is defined by the equation (R -x 2 + y 2 ) 2 + z 2 = r 2 . The parametric equations for such a torus are:

     x = (R + r cos θ) cos φ y = (R + r cos θ) sin φ z = r sin θ for (φ, θ) ∈ [0, 2π) 2 .
The angles φ and θ are respectively referred to as the toroidal and poloidal angles. A set of points with the same toroidal angle is called a meridian.

Model

We first investigate the behavior of system (1) with Approach B (using the geodesic distance) in the case of T 2 . Unlike in the cases of S 1 and S 2 presented in the sections 3 and 4, there exists no simple expression for the geodesic distance between two points on the torus. In 1903, Bliss studied and classified the different kinds of geodesic lines on the standard torus [START_REF] Bliss | The geodesic lines on the anchor ring[END_REF], using elliptic functions. Gravesen et al. determined the structure of the cut loci of a torus of revolution [START_REF] Gravesen | The cut locus of a torus of revolution[END_REF]. Several challenges arise when defining Approach B on T 2 . Firstly, computing the Riemannian distance between two points is highly non-trivial. One could consider approximating it numerically, but in the numerical discretization of equations ( 1)-( 5)-( 6), N (N -1)/2 geodesics would have to be computed per time-step. That would require tremendous computing power.

Secondly, assuming that one is able to efficiently compute the geodesics on T 2 , one must take into account the cut-loci of each point to ensure that the dynamics (1)-( 5)-( 6) are well-defined. A method to guarantee well-defined dynamics would be to use a bounded confidence model [START_REF] Hegselmann | Opinion dynamics and bounded confidence: models, analysis and simulation[END_REF], where the neighborhood of influence for an agent x i at point p is of smaller radius than the closest element in the cut locus of p. See section 2.2 for conditions on Ψ to make the right hand side of equation ( 1) continuous.

For simplicity, we thus focus on Approach A, where the dynamics are a function of the projection of each vector x j -x i onto the tangent space at x i . We will show that some restrictions still apply to the interaction function Ψ, but they are less restrictive and more easily determined than in Approach B.

Equations ( 1)-( 2) reads:

ẋi = N j=1 a ij Ψ( Π Tx i T 2 (x j -x i ) )ν ij , i ∈ {1, . . . , N }. ( 36
)
The vector ν ij depends on the influence that x j has over x i . It is zero if Π Tx i T 2 (x j -x i ) = 0, and it is a unit vector otherwise. Let N i be the set of points that have no influence on x i (see Table 1). Then, given i, j ∈ {1, . . . , N }, ν ij has the following expression:

ν ij = Π Tx i T 2 (xj -xi) Π Tx i M (xj -xi) if x j / ∈ N i 0 if x j ∈ N i . (37)
Let x i ∈ T 2 . We start by determining the set N i . For all i, we define the vectors u φi = cos φ i e x + sin φ i e y and u θi = cos θ i u φi + sin θ i e z , so that each agent's position vector reads:

x i = Ru φi + ru θi .
With these notations, u θi is the normal to the tangent space at the point x i . A basis for the tangent space at a point x i (φ i , θ i ) is given by the two tangent vectors t φi = (-sin φ i , cos φ i , 0) and t θi = (-sin θ i cos φ i , -sin θ i sin φ i , cos θ i ). Notice that x i , t φi = 0. Hence the condition Π Tx i T 2 (x j -x i ) = 0 reads:

x j , t φi = 0 x j -x i , t θi = 0.

After computations, we get:

x j , t φi = 0 ⇐⇒ sin(φ j -φ i ) = 0 ⇐⇒ φ j = φ i + kπ, k ∈ Z.
If φ j = φ i , the second condition becomes:

x j -x i , t θi = 0 ⇐⇒ sin(θ j -θ i ) = 0 ⇐⇒ θ j = θ i + kπ, k ∈ Z.
If φ j = φ i ± π, the second condition becomes:

sin(θ i + θ j ) = - 2R r sin θ i .
Notice that this last equation only has a solution if | sin θ i | ≤ r 2R . The set of positions that have no influence on x i thus comprises up to four points on the torus, depending on the values of r, R and sin θ i . We then have:

N i = {(φ i , θ i ), (φ i , -θ i ), (-φ i , -θ i -sgn(sin θ i ) arcsin(| 2R r sin θ i |), (-φ i , π - θ i + sgn(sin θ i ) arcsin(| 2R r sin θ i |))}.
To ensure the continuity of the right-hand side of equation ( 36), one must impose the conditions of table 1.

We now go back to equation (36). We study the specific case where Ψ ≡ Id, which indeed satisfies (1). Then the system becomes:

ẋi = Π Tx i T 2   N j=1 a ij (x j -x i )   . ( 38
)
Hence the velocity reads:

ẋi = N j=1 a ij (x j -x i ) - N j=1 a ij (x j -x i ), u θi u θi = α i -α i , u θi u θi -   N j=1 a ij   x i , t θi t θi
where α i := N j=1 a ij x j is the sum of the influences of all agents on agent i. Notice that with the same notation, the system does not reduce to the simple form ẋi = α i -α i , x i x i for the same dynamics on the sphere (see [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF]). This is due to the fact that on the torus, the position vector x i does not define the normal to the tangent space at x i , unlike in the cases of S 1 and S 2 .

The velocity of each agent is given by:

ẋi =   -φi sin φ i (R + r cos θ i ) -r θi sin θ i cos φ i φi cos φ i (R + r cos θ i ) -r θi sin θ i sin φ i r θi cos θ i   = φi (R + r cos θ i )t φi + r θi t θi . (39)
From ( 38) and (39) we get the angular velocities:

φi = 1 (R+r cos θi) N j=1 a ij (x j -x i ), t φi θi = 1 r N j=1 a ij (x j -x i ), t θi . ( 40
)
Notice that unlike in the case of S 2 , see equation (34), here the derivatives φi and θi are not singular. This makes numerical simulations straightforward, not requiring the approximations described in Section 4.2.

Properties

We now analyze the dynamics (1)-( 2)-(3) on T 2 . We identify families of initial conditions that trivialize the dynamics. Proposition 5.1. Consider the dynamics (1)-( 2)-( 3) on M = T 2 . Let P z := {(x, y, z) ∈ R 3 | z = 0}. Let x i (t) be the position of the ithe agent at time t. If for all i ∈ {1, . . . , N }, x i (0) ∈ T 2 ∩ P z , then for all t ≥ 0, for all i ∈ {1, . . . , N }, x i (t) ∈ T 2 ∩ P z .

Proof. Suppose that for all i ∈ {1, . . . , N }, x i (0) ∈ T 2 ∩ P z . Then for all i ∈ {1, . . . , N }, θ i (0) = 0 or θ i (0) = π. Hence, for all i, j ∈ {1, . . . , N },

t θi (0) =   0 0 ±π   and x j (0) -x i (0) =   (R + r cos θ j ) cos φ j -(R + r cos θ i ) cos φ i (R + r cos θ j ) sin φ j -(R + r cos θ i ) sin φ i 0  
From equation (40) we get: for all i ∈ {1, . . . , N }, θi = 0. By uniqueness of solution, for all i ∈ {1, . . . , N }, θ i (t) = θ i (0). All the initial velocities belong to the plane P z . Hence all agents remain on P z at all time.

Remark 5.1. As a consequence of Proposition 5.1, if all agents are initially in T 2 ∩ P z , all agents initially on the bigger circle θ = 0 remain on the major circle at all time and all agents on the minor circle θ = π remain on the minor circle at all time. In particular, if all agents are initially all on the same circle (i.e. for all i ∈ {1, . . . , N }, θ i = 0 or for all i ∈ {1, . . . , N }, θ i = π), then the torus dynamics simplify to the dynamics on S 1 given by (21) or (22).

Proposition 5.2. Consider the dynamics (1)-( 2)-( 3) on M = T 2 . Let φ ∈ [0, 2π] and let

P φ := {(x, y, z) ∈ R 3 | y = tan( φ)x}. If for all i ∈ {1, . . . , N }, x i (0) ∈ T 2 ∩ P φ, then for all t ≥ 0, for all i ∈ {1, . . . , N }, x i (t) ∈ T 2 ∩ P φ.
Proof. Suppose without loss of generality that φ = 0. Similarly to the proof for Proposition 5.1, we can show that for all i ∈ {1, . . . , N }, φi (0) = 0. By uniqueness of solution, for all i ∈ {1, . . . , N }, φ i (t) = φ i (0). Hence all agents remain in P φ at all time.

Remark 5.2. As a consequence of Proposition 5.2, if all agents are initially in T 2 ∩ P φ, all agents initially on the circle φ = φ remain on that circle at all time and all agents on the circle φ = -φ remain on that circle at all time. In particular, if all agents are initially all on the same circle (i.e. for all i ∈ {1, . . . , N }, φ i = φ or for all i ∈ {1, . . . , N }, φ i = -φ), then the torus dynamics simplify to the dynamics on S 1 given by (21) or (22).

Simulations

To assess the influence of the curvature of the manifold on the dynamics, we compare a simple case involving 3 agents evolving according to the interaction matrix given in equation ( 35). As in the case of S 2 , the dynamics on the torus do not give rise to periodic trajectories (as opposed to the dynamics in R 2 , see Theorem 6.2). Instead, since T 2 can locally be identified with R 2 , if the initial mutual distances are small enough, the dynamics resemble those in R 2 . More specifically, the trajectories are quasi-periodic with a gradual shift of the center of mass (see Figure 11). However, if the initial distances between agents are large, the geometry and curvature of the torus changes radically the behavior of the system.

Social choreographies

As seen in Sections 3 and 5.3, when the interaction matrix A satisfies certain properties, for instance given by ( 27) on S 1 or by (35) in R 2 , then the trajectories exhibit special properties of symmetry or periodicity. In [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF], configurations on S 2 in which all mutual distances between agents remain constant were named dancing equilibrium.

In this section, we investigate systems with similar properties of periodicity symmetry. We use the term social choreography, drawing a parallel with the well-known "n-body choreographies" discovered by Moore [START_REF] Moore | Braids in classical dynamics[END_REF][START_REF] Moore | New periodic orbits for the n-body problem[END_REF] in the context of point masses subject to gravitational forces. In the n-body problem, the interaction potentials between masses are predetermined, as they depend exclusively on the masses and distances between agents. Hence the conditions for a n-body choreography to occur only depend on the initial state of the system. In the case of social choreography, there are more degrees of freedom, as we design the interaction matrix as well as to set the initial conditions.

We study sufficient conditions on the interaction matrices for the trajectories of the system to be periodic or symmetric by focusing on the Euclidean space R 2 , with the specific choice of interaction potential Ψ ≡ Id. A future direction of this paper can consist in extending these results to general Riemannian manifolds. In R 2 and with Ψ ≡ Id, both approaches A and B are equivalent and the system simply reads as:

for all i ∈ {1, . . . , N }, ẋi = N j=1 a ij (x j -x i ). (41)
We define the kinetic energy E = E G = E P as in Definition ??:

E(t) := 1 2 N i=1 ẋi (t) 2 . (42)
A simple case of social choreography is that of a system with periodic trajectories, which we define as follows: Definition 6.1. Let (x i ) i=1...N be a solution of (41). We refer to the system as having periodic trajectories if there exists τ > 0 such that for all i ∈ {1, ..., N }, for all t > 0, x i (t + τ ) = x i (t).

We will examine possible periodic behaviors of the system in sections 6.2, 6.3 and 6.4.

Rotationally invariant system

We now give sufficient conditions on the interaction matrix and on the initial conditions for the system to be invariant by rotation. for the angle θ ∈ [0, 2π). Suppose that initially, the system is invariant by rotation of angle 2kπ N , that is:

for all i ∈ {1, . . . , N }, R( 2kπ N )x i (0) = x i+k (0) if i + k ≤ N x i+k-N (0) if i + k > N .
Suppose that the interaction matrix A is invariant by change of basis, i.e. P -1 k AP k = A. Then the system remains invariant by rotation of angle 2kπ N at all time:

for all t > 0, for all i ∈ {1, . . . , N }, R( 2kπ N )x i (t) = x i+k (t) if i + k ≤ N x i+k-N (t) if i + k > N .
Proof. Let A ∈ M N (R) be the interaction matrix, i.e. A = (a ij ) i,j=1,...N , and define D = diag( j a ij ). Let x = (x 1 , . . . , x N ) denote the set of all x i 's. It is a vector of length N with entries in R 2 . Let X ∈ M N ×2 (R) denote the corresponding matrix of R N ×2 such that for all i ∈ {1, . . . , N }, for all j ∈ {1, 2}, X ij is the j-th coordinate of x i . With these notations, Ẋ = ÃX, where à = A -D. We denote by (e 1 , . . . , e N ) the canonical orthonormal basis of (R) N such that X = N i=1 e i x T i . From the definition of the matrix X, the condition

for all i ∈ {1, . . . , N }, R( 2kπ N )x i (0) = x i+k (0) if i + k ≤ N x i+k-N (0) if i + k > N can be rewritten as: P k X(0) = (R( 2kπ N )X(0) T ) T . Let Y := P k X and Z := (R( 2kπ N )X T ) T .
From the theorem's hypotheses, Y (0) = Z(0). Let us show that Y and Z have the same evolution. One can easily prove that

P -1 k ÃP k if and only if P -1 k AP k . Then notice that Ẋ = ÃX = P -1 k ÃP k X.
From that we compute:

Ẏ = P k Ẋ = P k (P -1 k ÃP k X) = ÃP k X = ÃY. Similarly, Ż = (R( 2kπ N ) ẊT ) T = (R( 2kπ N )( ÃX) T ) T = (R( 2kπ N )X T ÃT ) T = ÃZ.
Since Y and Z satisfy the same differential equation and Y (0) = Z(0), then Y (t) = Z(t) for all t ≥ 0. This implies that at all time,

for all i ∈ {1, . . . , N }, R( 2kπ N )x i (t) = x i+k (t) if i + k ≤ N x i+k-N (t) if i + k > N . -4 -3 -2 -1 0 1 2 3 4 -3 -2 -

Unique orbit

Another example of social choreography is that of a system in which all agents share one unique orbit. Such choreographies have been discovered in the context of the n-body problem, for instance the "figure 8" orbit for three equal masses [START_REF] Moore | Braids in classical dynamics[END_REF].

Definition 6.2. Let (x i ) i=1...N be a solution of (41). We say that the system has a unique orbit if the orbits of all points are identical, i.e.

for all i, j ∈ {1, ..., N }, {z ∈ M |∃t > 0, x i (t) = z} = {z ∈ M |∃t > 0, x j (t) = z}.

To illustrate Theorem 6.1, we study the evolution of N agents initially positioned at regular intervals on a circle, with an interaction matrix and initial conditions given by:

A =             0 1 0 . . . 0 -
           
and for all i ∈ {1, . . . , N },

x i (0) = cos( 2iπ N ) sin( 2iπ N )
.

Notice that à = A, and the system satisfies the conditions of Theorem 6.1 with k = 1. Hence for all i ∈ {1, . . . , N -1}, R( 2π N )x i (t) = x i+1 (t) and R( 2π N )x N (t) = x 1 (t). The 2N -dimensional system then reduces to a 2-dimensional one for the two coordinates x 11 and x 12 of x 1 , and all the other variables can be recovered by rotation of x 1 :

ẋ1 = x 2 -x N = R( 2π N )x 1 -R(- 2π N )x 1 .
This can be written as:

ẋ11 ẋ12 = 0 -2 sin( 2π N ) 2 sin( 2π N ) 0 x 11 x 12 .
Solving this linear system yields:

x 11 (t) = x 11 (0) cos(2 sin( 2π N )t) -x 12 (0) sin(2 sin( 2π N )t) = cos(2 sin( 2π N )t) x 12 (t) = x 11 (0) sin(2 sin( 2π N )t) + x 12 (0) cos(2 sin( 2π N )t) = sin(2 sin( 2π N )t)

.

This proves that all agents share one common circular orbit, and their trajectories are periodic of period 2π(2 sin( 2π N )) -1 . Figure 15 provides a numerical illustration of this behavior, with 10 agents initially positioned at regular intervals on the unit circle.

Another interesting example is that of 3 agents interacting according to the interaction matrix given previously, which, reduced to N = 3, gives:

A =   0 1 -1 -1 0 1 1 -1 0   . ( 44
)
Theorem 6.2. Let N = 3. Consider the system (41) with interaction matrix given by (44). Then there exists a unique orbit shared by all agents, and all three trajectories are periodic.

Proof. The x and y-coordinates of the systems are decoupled, so that the 6-dimensional system can be reduced to two 3-dimensional ones. Notice that à = A. Then for each coordinate j ∈ {1, 2}, the system reads: with

  x 1j x 2j x 3j   (
e tA = 1 3   1 + 2 cos( √ 3t) 1 -cos( √ 3t) + √ 3 sin( √ 3t) 1 -cos( √ 3t) - √ 3 sin( √ 3t) 1 -cos( √ 3t) + √ 3 sin( √ 3t) 1 -cos( √ 3t) - √ 3 sin( √ 3t) 1 + 2 cos( √ 3t) 1 -cos( √ 3t) - √ 3 sin( √ 3t) 1 + 2 cos( √ 3t) 1 -cos( √ 3t) + √ 3 sin( √ 3t)   .
Due to the special structure of e tA , this can be rewritten as:

  x 1j x 2j x 3j   (t) = 1 3   x 0 1j x 0 2j x 0 3j x 0 2j x 0 3j x 0 1j x 0 3j x 0 1j x 0 2j     1 + 2 cos( √ 3t) 1 -cos( √ 3t) + √ 3 sin( √ 3t) 1 -cos( √ 3t) - √ 3 sin( √ 3t)   .
This shows that all three trajectories are periodic, or period 2π √ 3 . One can compute the positions of each agent after a third of a period and notice that:

  x 1j x 2j x 3j   (t + 2π 3 √ 3 ) = 1 3   x 0 1j x 0 2j x 0 3j x 0 2j x 0 3j x 0 1j x 0 3j x 0 1j x 0 2j     1 -cos( √ 3t) - √ 3 sin( √ 3t) 1 + 2 cos( √ 3t) 1 -cos( √ 3t) + √ 3 sin( √ 3t)   =   x 2j x 3j x 1j   (t).
This shows that there is one unique shared orbit.

Coupled periodic trajectories

Other conditions on the interaction matrix A give rise to different kinds of periodic behaviors. Here we provide sufficient conditions for the system to exhibit periodic trajectories, such that each orbit is shared by two agents.

Theorem 6.3 (Coupled periodic trajectories). Let N be even. Suppose that initially, the system is invariant by rotation of angle 4π N , that is:

for all i ∈ {1, . . . , N }, R( 4π N )x i (0) = x i+2 (0) if i + 2 ≤ N x i+2-N (0) if i + 2 > N .
Let a, b > 0 and let

A =             0
. . 0 -a 0             . ( 45
)
Then the system is periodic of period τ = π √ ab sin(2π/N ) . Furthermore, if N is divisible by 4, opposite agents share orbits two by two, i.e.: for all t > 0, for all i ∈ {1, . . . , N

}, x i (t + τ ) = x i+ N 2 (t), 2
and the kinetic energy is periodic with period τ /2.

Proof. First remark that the system satisfies the hypotheses of Theorem 6.1, so for all t > 0, for all i ∈ {1, . . . , N }, R(

N )x i (t) = x i+2 (t) if i + 2 ≤ N x i+2-N (t) if i + 2 > N . 4π
Hence the system is entirely known from the positions of the first two agents, since all others can be obtained by simple rotations. We show that this 2N -dimensional problem can be rewritten as a 4-dimensional one. Indeed, using the fact that x N = R(-4π/N )x 2 and

x 3 = R(4π/N )x 1 , the system ẋ1 = a(x 2 -x 1 ) -b(x N -x 1 ) ẋ2 = b(x 3 -x 2 ) -a(x 1 -x 2 )
becomes:

           ẋ1 = ẋ11 ẋ12 = a x 21 x 22 - x 11 x 12 -b cos( 4π N ) sin( 4π N ) -sin( 4π N ) cos( 4π N ) x 21 x 22 - x 11 x 12 ẋ2 = ẋ21 ẋ22 = b cos( 4π N ) -sin( 4π N ) sin( 4π N ) cos( 4π N ) x 11 x 12 - x 21 x 22 -a x 11 x 12 - x 21 x 22 .
This can be rewritten in matrix form as:

    ẋ11 ẋ12 ẋ21 ẋ22     = A 4     x 11 x 12 x 21 x 22     (46)
where

A 4 :=     -a + b 0 a -b cos( 4π N ) -b sin( 4π N ) 0 -a + b b sin( 4π N ) a -b cos( 4π N ) -a + b cos( 4π N ) -b sin( 4π N ) a -b 0 b sin( 4π N ) -a + b cos( 4π N ) 0 a -b         x 11 x 12 x 21 x 22     . x 1 x 2 x 3 x 4 x 5
x 6

x 7

x 8

x 9

x 10

x 1

x 2

x 3

x 4

x 5

x 6

x 7

x 8 One can easily show that this reduced interaction matrix A 4 has two purely imaginary conjugate eigenvalues, iλ and -iλ, each of multiplicity 2, where λ = 2 √ ab sin( 2π N ). Hence the solution of the system (46) can be written as a weighted sum of the functions t → cos(λt) and t → sin(λt). This implies that the system is periodic, of period

τ = 2π λ = π √ ab sin( 2π N )
.

Furthermore, if N is divisible by 4, according to Theorem 6.1, x N 2 +1 = -x 1 and x N 2 +2 = -x 2 . This implies that for all t > 0, x 1 (t + τ ) = -x 1 (t) = x N 2 +1 (t) and x 2 (t + τ ) = -x 2 (t) = x N 2 +2 (t), so the agents x 1 and x N 2 +1 share an orbit, as well as all pairs of agents x i and x N 2 +i for i ∈ {1, . . . , N 2 }. As a consequence, the kinetic energy is periodic, of period τ = π/(2 √ ab sin( 2π N )). If N is divisible by 4, every half period, the system is rotated by an angle π, so the kinetic energy is periodic with period τ /2. Remark 6.1. Notice that the agents sharing orbits do not interact with one another, as shown in Figure 16.

An example of such a choreography is given in Figure 17. Remark 6.2. As a slight generalization, we provide numerical simulations illustrating a similar behavior, but with slightly different conditions: the periodic evolution of 9 agents on three distinct orbits shared three by three, see figures 18 and 19.

Helical trajectories

In sections 6.2 and 6.3, we provided conditions for the trajectories of the system to be periodic. Here, we explore further the notion of periodicity by studying systems with drift, displaying helical -0. trajectories but periodic kinetic energy. Definition 6.3. Let (x i ) i=1...N be a solution of (41). We call the corresponding trajectories helical trajectories if there exists v ∈ R 2 and τ ∈ R * such that for all i ∈ {1, ..., N }, for all t > 0, x i (t + τ ) = x i (t) + τ v.

Notice that this definition generalizes the notion of periodic trajectories recalled in Definition 6.1, which corresponds to the case v = 0. When v = 0, the system has a drift term, meaning that the relative positions between agents remain periodic but their absolute positions evolve in space. (47)

Then the system exhibits helical trajectories.

Proof. First notice that the first and second components x i1 and x i2 of the i-th agent's position are decoupled, so that the system in matrix form reads There is one eigenvector associated with λ 1 : v 1 := (1, 1, 1, 1) T . One can show that the vectors x(t) = v 1 and x(t) = v 1 t + ν are both solutions of System (48), where, denoting ∆ := bcd -abc + abd -acd, ν := 1 ∆ (ab + bc + ∆, ab -cd + ∆, ab + ad + ∆, ∆) T .

Let v 2 denote the eigenvector associated with λ 2 and let v R 2 and v I 2 denote respectively its real and imaginary components, i.e. v 2 := v R 2 + iv I 2 . Then the solution of System (48) can be written as:

x j (t) = C j 1 v 1 + C j 2 (v 1 t + ν) + C j 3 v R 2 cos(λ 2 t) -v I 2 sin(λ 2 t) + C j 4 v R 2 sin(λ 2 t) + v I 2 cos(λ 2 t)
where (C 1 , C 2 , C 3 , C 4 ) ∈ R 4 are constants depending on the initial conditions. Let τ = 2π λ2 . Then for all t > 0, for all i ∈ {1, . . . , 4}, for all j ∈ {1, 2}, x ij (t + τ ) = x ij (t) + C j 2 τ . This can be rewritten

A Appendix

An example using Approach A which shows unexpected behavior in the first example (A.1), as well as the interactions matrix and initial positions used for simulations shown in Figure 7 and9.

Example A.1. 15 agents with a general interaction matrix A. We use the same notion of a general interactions matrix as used in [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF]. The interactions matrix A is composed of integers a ij that are uniformly chosen between -5 and 5 inclusive. Ψ ≡ Id. Generally, a system with this kind of interaction matrix will exhibit simple oscillating kinetic energy, as in [START_REF] Caponigro | A nonlinear model of opinion formation on the sphere[END_REF]. We show a rare simulation using this general interactions matrix A in the appendix (Figures 22 and23). Example A.2. Five agents with a general interaction matrix A and Ψ as defined in (25). Simulations are shown in Figure 7 and 9 with a ∈ { π 4 , π 2 , 3π 4 }.

A =       -5 -1 1 0 4 -1 -5 1 1 0 -2 0 -5 -1 -2 -4 -1 -5 3 -4 -4 1 -2 2 1       X =                           1.
A =       -5 -1 1 0 4 -1 -5 1 1 0 -2 0 -5 -1 -2 -4 -1 -5 3 -4 -4 1 -2 2 1       X =       6

Figure 4 :

 4 Figure 4: Initial (empty circles) and final positions (filled circles) of 4 agents initially on the vertices of a rectangle with Approach B (left) and Approach A (right), with A = 1, Ψ ≡ Id (Approach A) and Ψ = Ψ 3π/4 (Approach B) (see equation (25)). Notice that with Approach A, initial and final positions are identical since any rectangle configuration is an equilibrium. However, with Approach B, the system reaches a square configuration, the only possible equilibrium with pairwise distinct positions.

Figure 5 :Figure 6 :

 56 Figure 5: Evolution of the system (26) with Approach A (left) Approach B (center) when the interaction matrix satisfies condition (27) for the projection distance. Right: Kinetic energy.

Figure 7 :

 7 Figure 7: A dancing equilibria for Approach A. The energy becomes constant in time after initial fluctuations.

Figure 8 :

 8 Figure 8: Choice of Ψ functions determines the energy of the system, for the case of a = 3π 4 the system forms an antipodal equilibrium.

Figure 9 :

 9 Figure 9: A comparison of the effect of the choice of influence function for Approach B. For a = 3π 4

Figure 11 :Figure 12 :

 1112 Figure 11: Trajectories of three agents interacting according to the matrix A given in (35). Left: Dynamics in R 2 , with periodic trajectories on a unique orbit. Center: Dynamics on M = T 2 with small initial mutual distances. Right: Dynamics on M = T 2 with large initial distances.

Theorem 6 . 1 .

 61 Let k ∈ N such that k divides N . Let P k = 0 I N -k I k 0 be the matrix of change of basis from (e 1 , . . . , e N ) to (e k , . . . , e N , e 1 , . . . , e k-1 ). Let R(θ) denote the rotation matrix in R 2

Figure 13 :Figure 14 :

 1314 Figure 13: Left: Evolution of 12 agents with the conditions of Theorem 6.1, with k = 3, resulting in diverging trajectories. Dark to light color scale indicates earlier to later time. Right: corresponding exploding kinetic energy. The interaction matrix A and the initial positions were generated according to a random algorithm, with the conditions of Theorem 6.1.

Figure 15 :

 15 Figure 15: Evolution of 10 agents with initial conditions and interaction matrix given in (43). The agents have periodic trajectories along one shared circular orbit.

Figure 16 :

 16 Figure 16: Left: Directed graph corresponding to the matrix A given in (43). Full arrows represent positive coefficients (a ij > 0) while dashed ones represent negative coefficients (a ij < 0). Right: Weighted directed graph corresponding to the matrix A given in (45). Thin arrows represent the weighted edges |a ij | = a while bold ones represent the weight |a ij | = b. Nodes with the same color and symbol share orbits but are not directly connected in the graph.

Theorem 6 . 4 .

 64 Sufficient conditions for helical trajectories. Let N = 4. Let (a, b, c, d) ∈ (R + ) 4 such that the interaction matrix reads

Figure 18 :Figure 19 :

 1819 Figure 18: Left: evolution of 9 agents with periodic trajectories, each orbit shared by 3 agents. Right: periodic kinetic energy.

Figure 22 :

 22 Figure22: Energy of the system using Approach A, 15 agents, and a general interaction matrix (left). A snapshot of the energy oscillations to match with trajectories in Figure23(right).

Figure 23 :

 23 Figure 23: An agent's trajectory simulated with Approach A, 15 agents, and a general interaction matrix. The trajectory in shown the top right is of a second agent. The agents oscillate with amplitudes that increase with time, eventually the trajectory approximates a great circle, after which the oscillations resume with smaller amplitudes.

  Periodic trajectories of 8 agents sharing orbits two by two, in the situation of Theorem 6.3. Matrix A from (45) was constructed with (a, b) = (1, 3). The initial positions x 1 (0) and x 2 (0) were randomly generated and the other 6 were obtained by rotation. The period is τ = 2π/ √ 6. Right: Corresponding kinetic energy, of period τ /2.

				Time = 1.7								Kinetic Energy
	0.6
								1	12
								2
								3
	0.4							4 5	10
								6
	0.2							7
								8	8
	0
									6
	-0.2
									4
	-0.4
									2
	-0.6
									0
	8	-0.6	-0.4	-0.2	0	0.2	0.4	0.6	0.8	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6
														Time
	Figure 17: Left:

as: for all i ∈ {1, . . . , 4}, for all t > 0,

Theorem 6.5. A system with helical trajectories has periodic kinetic energy.

Proof. Supose that (x i ) i=1...N has helical trajectories, i.e. there exists τ ∈ R, v ∈ R 2 such that for all i ∈ {1, . . . , N }, for all t ≥ 0,

and so E(t + τ ) = E(t).