Petar Mircheski

Chimera States

Mircheski Petar

April 28, 2025

blog image

How Chimera States Adapt to Frequency Heterogeneity: A Self-Organizing Dance of Oscillators

Chimera states are one of the most intriguing phenomena in the world of nonlinear dynamics. They are composed of identical oscilaltors connected in a ring. Under certain conditions, they mysteriously split into two groups: one synchronized and coherent, the other disordered and chaotic. This coexistence of coherence and incoherence, is what defines a chimera state.

Originally reported by Kuramoto and Battogtokh in 2002 and named by Abrams and Strogatz in 2004, chimera states have fascinated scientists across disciplines. They've been observed in chemical reactions, mechanical systems, lasers, and more. Over the years, researchers have dug into how different types of connections, oscillator behaviors, and dimensional setups influence these states.

But what happens when the oscillators aren't quite identical? What if each one has a slightly different natural rhythm?

The Ring of Oscillators

We begin with a standard setup: a ring of phase oscillators, where each oscillator's behavior is described by a simple mathematical model.

θ˙(x,t)=ω(x)+∫−ππG(x−x′)sin⁡(θ(x′,t)−θ(x,t)−α)dx′.\dot{\theta}(x, t) = \omega(x) + \int_{-\pi}^{\pi} G(x-x') \sin(\theta(x', t) - \theta(x, t) - \alpha) d x'.

These oscillators are connected nonlocally, each one of them feels the influence of its neighbors, but not just the immediate ones. The coupling strength is determined by the function G(x)G(x), and the phase lag is given by α\alpha. In our case we define a ring of 256 oscillators, connected through a cosine distance kernel,

G(x−x′)=(1+Acos⁡(x−x′))12π,G(x-x') = (1 + A \cos(x - x')) \frac{1}{2\pi},

and fix the system parameters to α=π2−0.18\alpha = \frac{\pi}{2} -0.18 and A=0.995A = 0.995.

from typing import Callable

import numba
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp


def calc_dy(A: float, N: int) -> np.ndarray:
    """
    Calculate the distance kernel delta_y for a system of oscillators.
    This function generates a coupling kernel based on a cosine function,
    which is then normalized.

    :param A: A scaling parameter that affects the coupling strength between oscillators.
    :param N: The number of spatial grid points, which corresponds to the number of oscillators.
    :returns: A 1D numpy array of length N containing the computed coupling kernel values, scaled to ensure normalization.
    """
    x = np.linspace(0, 2 * np.pi, N, endpoint=False)
    delta_y = (1 + A * np.cos(x)) / (2 * np.pi)
    return np.array(delta_y * ((2 * np.pi) / N))


def calculate_dy_matrix(A: float, N: int) -> np.ndarray:
    """
    Generate a distance matrix delta_y for a system of oscillators,
    representing the pairwise coupling strengths between all oscillators.
    The matrix is optimized for use in coupling computations within the Kuramoto model.

    :param A: Scaling parameter for the cosine term in the kernel function.
    :param N: The number of spatial grid points (or oscillators).
    :returns: A 2D numpy array of shape [N x N] where each element (i, j) represents the
              coupling strength from oscillator j to oscillator i, based on the rolled distance kernel.
    """
    dy = calc_dy(A, N)
    dy_matrix = np.zeros((N, N))
    for i in range(N):
        dy_matrix[i, :] = np.roll(dy, i)
    return dy_matrix


def kuramoto_kernel(
    omegas: np.ndarray, kernel: np.ndarray, N: int, alpha: float = 0.0
) -> Callable[[np.ndarray], np.ndarray]:
    """
    Generate the Kuramoto differential equation function for a system of coupled oscillators,
    given their natural frequencies, a coupling kernel, and a phase shift parameter.
    This function uses Numba's jit compiler for efficient parallel computation.

    :param omegas: A 1D numpy array of length N containing the natural frequencies of the oscillators.
    :param kernel: A 2D numpy array of shape [N x N], representing the coupling strengths between oscillators.
    :param N: The number of oscillators in the system.
    :param alpha: A phase shift parameter. Default is 0.0.
    :returns: A function that, when provided with the state (phase) array, returns the rate of change of
              phases based on the Kuramoto model equations.
    """

    @numba.jit(nopython=True, parallel=True)
    def kuramoto(state: np.ndarray) -> np.ndarray:
        """
        The internal function generated by kuramoto_kernel. This function calculates the time derivative
        (rate of change) of the oscillator phases based on their natural frequencies, coupling kernel, and
        phase differences.

        :param state: A 1D numpy array of length N containing the current phases of the oscillators.
        :returns: A 1D numpy array of length N representing the rate of change of each oscillator's phase.
        """
        dstate = np.zeros_like(state)
        for i in numba.prange(N):
            phase_diff = state - state[i] - alpha
            dstate[i] = np.sum(
                kernel[i, :] * np.sin(phase_diff)
            )  # Sum product for row i
        return omegas + dstate

    return kuramoto  # type:ignore


np.random.seed(seed=1338)
N = 256
A = 0.995
alpha = (np.pi / 2) - 0.18
kernel = calculate_dy_matrix(A, N)
dt = 0.1

In a perfectly uniform setup, the position of the chimera state on the ring is completely determined by the system's initial conditions. That's because the system is spatially homogeneous meaning, no location is special.

Chimera States Without Heterogeneity

We begin with this uniform system ω(x)=0\omega(x) = 0 and observe the classic chimera behavior: a stable pattern with a coherent region and an incoherent region coexisting side by side. These patterns arise spontaneously and can be visualized in terms of phase snapshots, local coherence, and even instantaneous frequencies across the ring.

evolve_time = 500.0
t_evolve_eval = np.arange(0, evolve_time + dt, dt)
t_evolve_span = (0, evolve_time + dt)

omegas = np.zeros(N)
initial_conditions = (
    6.0
    * np.exp(-0.76 * (np.linspace(-np.pi, np.pi, N, endpoint=False)) ** 2)
    * (np.random.random(N) - 0.5)
)

kuramoto = kuramoto_kernel(omegas, kernel, N, alpha=alpha)
sol = solve_ivp(
    lambda t, y: kuramoto(y),
    t_evolve_span,
    initial_conditions,
    t_eval=t_evolve_eval,
    method="RK45",
    rtol=1e-7,
    atol=1e-10,
)
init_state = sol.y.T


fig, ax = plt.subplots(2, sharex=True)
heatmap = ax[0].imshow(
    init_state % (2 * np.pi), aspect="auto", interpolation="none", cmap="hot"
)

ax[0].set_ylabel(r"$t$")
ax[0].set_xlabel(r"$i$")
ax[0].set_title(r"$a)$", loc="left")


# Scatter plot
ax[1].scatter(np.arange(N), init_state[5000] % (2 * np.pi))
ax[1].set_xlabel(r"$i$")
ax[1].set_yticks([0, np.pi, 2 * np.pi])
ax[1].set_yticklabels([r"$0$", r"$\pi$", r"$2\pi$"])
ax[1].set_title(r"$b)$", loc="left")


cax = plt.axes((0.85, 0.1, 0.075, 0.8))
cbar = fig.colorbar(heatmap, cax=cax)
cbar.set_label(r"$\phi$")
cbar.set_ticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi])
cbar.set_ticklabels(
    [r"$0$", r"$\frac{\pi}{2}$", r"$\pi$", r"$\frac{3\pi}{2}$", r"$2\pi$"]
)

plt.subplots_adjust(bottom=0.1, right=0.8, top=0.9)

plt.show()
Chimera state in a ring of oscillators

But what happens when we gently nudge the system out of uniformity?

Introducing Frequency Heterogeneity

To probe the system's adaptability, we introduce a small frequency heterogeneity—a subtle difference in the natural frequencies of the oscillators. This can be thought of as introducing a "preference" in certain parts of the ring. The mathematical form of this perturbation is:

ω(x)=ω0+ϵh(x−γ)\omega(x) = \omega_0 + \epsilon h(x - \gamma)

Here, ϵ\epsilon is a small number (think of it as a tiny tweak), and h(x)h(x) is a function that defines how the frequency varies in space. The parameter γ\gamma determines where this heterogeneity is placed on the ring.

We explore two types of heterogeneity:

1. Point heterogeneity

A spike at a single location.

h(x)=δ~(x)h(x) = \tilde\delta(x)

new_omegas = np.zeros(N)
new_omegas[0] += 0.2
kuramoto = kuramoto_kernel(new_omegas, kernel, N, alpha=alpha)

t_evolve_eval = np.arange(0, evolve_time * 6 + dt, dt)
t_evolve_span = (0, evolve_time * 6 + dt)

perturbed_sol = solve_ivp(
    lambda t, y: kuramoto(y),
    t_evolve_span,
    init_state[-1],
    t_eval=t_evolve_eval,
    method="RK45",
    rtol=1e-7,
    atol=1e-10,
)
perturbed_state = perturbed_sol.y.T

fig, ax = plt.subplots(2, sharex=True)
heatmap = ax[0].imshow(
    perturbed_state % (2 * np.pi), aspect="auto", interpolation="none", cmap="hot"
)

ax[0].set_ylabel(r"$t$")
ax[0].set_xlabel(r"$i$")
ax[0].set_title(r"$a)$", loc="left")


# Scatter plot
ax[1].scatter(np.arange(N), perturbed_state[-1] % (2 * np.pi))
ax[1].set_xlabel(r"$i$")
ax[1].set_yticks([0, np.pi, 2 * np.pi])
ax[1].set_yticklabels([r"$0$", r"$\pi$", r"$2\pi$"])
ax[1].set_title(r"$b)$", loc="left")


cax = plt.axes((0.85, 0.1, 0.075, 0.8))
cbar = fig.colorbar(heatmap, cax=cax)
cbar.set_label(r"$\phi$")
cbar.set_ticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi])
cbar.set_ticklabels(
    [r"$0$", r"$\frac{\pi}{2}$", r"$\pi$", r"$\frac{3\pi}{2}$", r"$2\pi$"]
)

plt.subplots_adjust(bottom=0.1, right=0.8, top=0.9)

plt.show()
Chimera state in a ring of oscillators

2. Cosine heterogeneity

h(x)=cos⁡(x)h(x) = \cos(x)

x = np.linspace(-np.pi, np.pi, N, endpoint=False)
new_omegas = 0.01 * np.cos(x + np.pi)
kuramoto = kuramoto_kernel(new_omegas, kernel, N, alpha=alpha)

t_evolve_eval = np.arange(0, evolve_time * 3 + dt, dt)
t_evolve_span = (0, evolve_time * 3 + dt)

perturbed_sol = solve_ivp(
    lambda t, y: kuramoto(y),
    t_evolve_span,
    init_state[-1],
    t_eval=t_evolve_eval,
    method="RK45",
    rtol=1e-7,
    atol=1e-10,
)
perturbed_state = perturbed_sol.y.T

fig, ax = plt.subplots(2, sharex=True)
heatmap = ax[0].imshow(
    perturbed_state % (2 * np.pi), aspect="auto", interpolation="none", cmap="hot"
)

ax[0].set_ylabel(r"$t$")
ax[0].set_xlabel(r"$i$")
ax[0].set_title(r"$a)$", loc="left")


# Scatter plot
ax[1].scatter(np.arange(N), perturbed_state[-1] % (2 * np.pi))
ax[1].set_xlabel(r"$i$")
ax[1].set_yticks([0, np.pi, 2 * np.pi])
ax[1].set_yticklabels([r"$0$", r"$\pi$", r"$2\pi$"])
ax[1].set_title(r"$b)$", loc="left")


cax = plt.axes((0.85, 0.1, 0.075, 0.8))
cbar = fig.colorbar(heatmap, cax=cax)
cbar.set_label(r"$\phi$")
cbar.set_ticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi])
cbar.set_ticklabels(
    [r"$0$", r"$\frac{\pi}{2}$", r"$\pi$", r"$\frac{3\pi}{2}$", r"$2\pi$"]
)

plt.subplots_adjust(bottom=0.1, right=0.8, top=0.9)

plt.show()
Chimera state in a ring of oscillators

The Self-Adaptive Behavior of Chimera States

When we introduce point heterogeneity (like a single oscillator ticking faster), the chimera state responds in a surprising but consistent way. The incoherent region starts drifting—migrating toward the faster oscillator. Meanwhile, the coherent region drifts in the opposite direction, settling in the slower area.

This movement continues until the system reaches a new equilibrium, where the incoherent region hugs the heterogeneity, and the coherent region finds its new home far away from it.

When we switch to the smoother cosine heterogeneity, the same behavior emerges only this time, more gracefully. The incoherent region moves toward the high-frequency side of the cosine wave, while the coherent region moves toward the low-frequency trough. Again, a new steady state is reached, but now perfectly aligned with the shape of the heterogeneity.

Coherence Loves Low Frequencies

The takeaway from all of this? Regardless of how the heterogeneity is shaped—spiky or smooth—the chimera state adapts in the same fundamental way:

  • Incoherent regions gravitate toward high natural frequencies.
  • Coherent regions drift toward low natural frequencies.

This self-adaptive behavior is both robust and elegant. It shows that even in complex dynamical systems, a simple principle can emerge: the system self-organizes in a way that aligns its internal structure with the external heterogeneity.

We call this phenomenon spatial phase locking to frequency heterogeneity.

Why Does This Matter?

Understanding how chimera states adapt to heterogeneity has implications in various fields. Knowing how these patterns respond to imperfections or design tweaks could help us control, harness, or stabilize such systems in the real world.