Authors: Eryn Cangi, Daniel Abrams

for Northwestern University Center for Interdisciplinary Exploration and Research in Astrophysics REU, 2015

This content was originally published at the CIERA website. Changes are no longer allowed to the original website, so it has been reproduced here by Eryn to allow for some minor content updates (equation corrections and an additional image).

The information below is also available as a poster.

Introduction and Science Goals

Our goal: to develop a nonlinear model that can be used in place of more complex models to describe the time evolution of synchronized astrophysical systems, such as moon formation, ring formation, spiral arm patterns or tidal locking of satellites.

The inspiration to approach astrophysical problems from this point of view came from Dr. Abrams’ experience with coupled nonlinear oscillators. Nonlinear dynamics provide a way to utilize “big picture” characterizations about a system alternative to classical mechanics, statistical mechanics and electrodynamics.

While theories of synchronization in two- or three-body astronomical systems are well understood, a generalization to many-bodied systems remains largely unexplored. Historically, problems of resonant capture among astronomical bodies have been treated using methods from conservative classical mechanics. We investigate the possibility of using non-conservative models together with numerical methods to understand the phenomenon of resonant capture in large-scale structures such as rings, planetary systems and galactic spiral arms. In particular, we focus on N-body (N\(\gg\) 1) dissipative systems such as circumplanetary discs and use methods drawn from the study of coupled oscillators.

One such method is inspired by the Kuramoto model, which describes mean field behavior in large ensembles of coupled nonlinear oscillators. The Kuramoto model can be modified to allow for non-mean field coupling, leading to the existence of chimera states, in which most of the oscillators synchronize. These chimera states can appear as clusters or spirals of synced oscillators, and may be suggestive of objects in astronomical contexts. We developed a mean-field model for N small particles in a dust ring around a massive planet, which can be tested through numerical experiments using code developed in MATLAB and Python. We are currently in the process of refining the code, after which the numerical experiments will begin.

The Kuramoto Model as Motivator

The Kuramoto model is a nonlinear differential equation that relates the frequency of, or rate of change in, an oscillation to the current position in the oscillation cycle. What is this “oscillation?” It can be many things. Some examples include:

Mathematics of the Kuramoto model

The Kuramoto equation is customarily written as follows: \(\frac{d\theta_i}{dt} = \omega_i + \frac{K}{N}\sum^N_{j=1} \sin(\theta_j-\theta_i)\)

The components contributing to the rate of change are stated in the right hand side of the equation.

  • \(\omega\): the object’s natural frequency–the normal speed at which it will perform its cycle, absent driving or dampening forces.
  • summation term: an adjustment to the the frequency that results from forces between the object and all the other objects in the system.

The variables in the system are further defined as below.

  • K: a constant representing the strength of the coupling (how strongly the objects affect each other).
  • N: the number of objects in the system
  • \(\theta\): the current phase of the object, that is, where it is in its cycle.

The subscript i denotes a specific object in the system. This means that every object in the system will have its behavior described by this equation. The subscript j varies between 1 and N and designates “other” objects: this is where the summation comes in. The object’s natural frequency is adjusted up or down depending on its position relative to another object in the system. The the fact that the difference occurs inside of a sine function indicates that the sign of the adjustment will change depending on whether the jth object is ahead or behind the current object i.

Visualization of behavior

128-bodies

So, what does this say about the system? The upshot is that if the coupling strength K is high enough, objects in the system–no matter if they are living creatures or inanimate objects–will adjust their cycles to match their neighbors! When plotted, this behavior creates intricately and chaotically overlapping paths before synchronization, as shown at left.

In the plot of the phase of oscillators over time, each colored line represents a different object in the system, with its own unique natural frequency. At the beginning of the simulation (or, in real life, when the objects start oscillating), all objects start at a given point in their cycle. But over time, minor adjustments to their frequencies occur due to effects from the others surrounding them. The end result is that all the objects coalesce to a common phase, which is synchronization.

If we envision the objects as particles on a circular track (for example, the unit circle), the objects will all be clustered around the phase angle of synchronization and will move around the circle together, like a group of runners in a race who are neck and neck.

We can use the fact that the objects are all at very slightly different radii around the circle to describe the degree of synchrony. By performing a vector sum of the components of the radius of each object, we obtain components for an “average” radius of all the oscillators. The maximum of this radius will be the radius of the circular track. When this maximum is equal to the “average” radius, complete synchrony is achieved. The radius approaching 0 would correspond to a state in which the objects are randomly distributed about the circle–i.e., the beginning of the simulation, when they are all at different phases. An animation illustrating this idea follows (Source: Wikipedia, credit to Stewart Heitmann).

Fireflies

The canonical example of this mathematical model is the blinking of certain firefly species, specifically Pteroptyx malaccae and Pteroptyx cribellata.

If we were to travel to the homes of these fireflies and wait for twilight to see their display, we would at first see numerous fireflies begin blinking randomly. Over time, they would eventually synchronize, blinking in perfect unison. For the fireflies, it’s a useful mating display. For mathematicians, it’s a fascinating example of non-linear dynamics.

Image credit: Lightning bugs in York, PA, Flickr user tom.arthur, 2009

 

 

N-Body Simulation and a Nonlinear Model

The initial stage of the project involved programming a custom N-body simulation (N = number of bodies) in MATLAB to solve Newton’s equation of motion: \(\frac{d\vec{r}}{dt} = -\frac{GM}{d^3}\vec{r}\)

Where \(\vec{r}\) is a vector identifying the bodies’ locations in 2D space and \( \left|\vec{r}\right|\) is the separation between any two given bodies. The simulation uses a 4th order Runge-Kutta method of integration. Because the number of equations to be solved is extremely large, the simulation utilizes vectorization of the equation for legibility and code efficiency.

A note on the choice of language: The importance of vectors in this problem influenced our decision to use MATLAB. Initial attempts at the code were made in Python using the Numpy module, but the code soon became difficult to easily manage while the building process was still underway. In the future, we intend to translate the simulation into Python as well.

N-body simulation without collisions

Plot of the orbital radii of mercury, venus, earth and mars over time
With the completion of the basic N-body simulator without collision handling, a fairly simple simulation of the orbits of the 8 planets can be run as proof of concept. The simulation tracks the distance of each planet from the sun over time, with each planet’s total distance oscillating about its mean. The plot at left shows an example of the results of this simulation (with distances in astronomical units) over a time span of three years. Note that for this particular simulation, all 8 planets and the sun were simulated, but only the four inner planets are plotted as the significant increase in distance of the outer planets would relegate all the inner planets to a very narrow region at the bottom of the plot.


An alternative, and perhaps more attractive, method for visualization of the data is with animation–especially since the simulation concerns bodies moving through space. Animation for the N-body simulation will be completed after further optimization and extension of the main features of the code. An example animation for a two-body system, the sun and the Earth, is shown above.

Collision handling

In order to simulate interactions of bodies in close proximity, such as those in a circumplanetary dust disk, the code must handle collisions. This allows simulations of larger N and at closer distances to a central body (for instance, distance < 1 AU). In the original simulation without collision handling, collisions would result in division by zero errors as d → 0. The collision is assumed to be inelastic; therefore the energy of the two-particle system is reduced by some fractional amount after the collision and the remaining energy is randomly distributed between the two bodies within a certain range. The range is required as decreasing or increasing a velocity by too much would violate conservation of momentum.

Our collision handling can be outlined as such:

  1. A condition is set for what constitutes collision–1000 km is initially used to speed debugging and ensure collisions actually occur, since the orbits are over very large distances.
  2. During integration, we use ode45’s options and odeset to configure an event function that checks for collisions at every step.
  3. In the event of a collision, integration is stopped and the following occurs:
    1. The code calls another function whose job is to solve a system of equations for the final velocities of each body. The system is defined by conservation of momentum along all axes, dissipation of energy and division of the remaining total energy to each of the two bodies using a random fractional value within a certain range.
    2. Using solution, the velocities for each body are inserted into a new set of initial conditions for all the objects in the system.
    3. Integration is restarted using the new initial conditions, picking up from the time step where it previously left off.

Allowing changes in velocity as described enables the bodies to alter the radius of their orbits around the central body. Lower orbits indicate higher kinetic energy of the body, while higher orbits indicate lower kinetic energy. The effect is that the bodies will speed up or slow down as a result of the collision–similar to the way that bodies obeying the Kuramoto model will alter their frequencies according to the presence of neighbors.

Custom nonlinear model of force interactions

We are currently working to adapt the code to use a new differential equation for the changes in orbital speed of a collection of many bodies around a much larger body—for example, a dust ring around a planet. This equation will be reminiscent of the Kuramoto model: \(\frac{d\theta_i}{dt} = \omega_i + \frac{GK}{4R^3} \int^{2\pi}_0 \frac{\rho(\theta ‘)d\theta ‘}{\sin^2(\frac{1}{2}\left|\theta – \theta ‘\right|)}\)

cartoonThe integral comes from a summation over N orbiting bodies and their interactions with one another. The sine is part of the term for the separation of two orbiting bodies and arises from the geometry (see cartoon at left). Lastly, the term [late\(x]R^3\) in the denominator arises from the \(R^2\) in the gravitational force and the \(R\) in \(R \frac{d\theta}{dt}\), the velocity of the particle, and our assumption that in a dissipative system, the velocity, rather than the acceleration, is proportional to the force. In the cartoon at right, two dust grains (gray) orbit a larger central body (perhaps a star).

Like the basic Kuramoto model, this equation expresses the idea that the frequency of an oscillator is modified from its natural frequency by a term comprised of effects from all other oscillators in the system. In this case, when N will be very large, the sum can be approximated with an integral. R represents the radius to the orbited body, G is the gravitational constant, K is a scaling constant and \(\rho\) is the density of material in the disk.

Future Work

Our current goal is code optimization and implementation of our custom mathematical model. In addition to this goal, we will continue refining the code and adding more layers of complexity, including:

  • Implementation of the custom nonlinear model in the N-body simulation code
  • Conversion to polar and spherical coordinates
  • Refinements to simulation code to improve efficiency and performance
  • Expansion to large N simulations to investigate:
    • Galactic spiral arms
    • Ring formation around massive planets
    • Circumstellar discs and planetary formation

References

Gladman, Brett, D. Dane Quinn, Philip Nicholson, and Richard Rand. “Synchronous Locking of Tidally Evolving Satellites.” Icarus 122 (1995): 166-92.

Strogatz, Steven H. “From Kuramoto to Crawford: Exploring the Onset of Synchronization in Populations of Coupled Oscillators.” Physica D 143 (2000): 1-20.

Images

Contact

If you would like to contact us about our project, please use the contact information below.

Eryn Cangi
See contact information here

Dr. Daniel Abrams
dmabrams *at* northwestern *dot* edu
Associate professor, Department of Engineering Sciences and Applied Mathematics
Northwestern University

Support & Funding

CIERA logo NSF logo Northwestern logo