Active matter

Active matter: self-propelled particles, flocking, polar fluids, chemotaxis, squirmers, and Active System B+.

ActiveBeadElement(label, visual=None, *, x=0.0, y=0.0, theta=0.0, v0=1.0, k=10.0, D_t=0.1, D_r=1.0)

Bases: Element

Bead of an active polymer, an active Brownian particle bonded to its neighbours along a chain by springs.

The bead at (x, y) swims at the speed v0 along its orientation theta, which diffuses with D_r, its position diffuses with D_t, and harmonic springs of stiffness k and zero rest length pull it towards its neighbours along the chain, k (r_{i+1} + r_{i-1} - 2 r_i), the Rouse chain. The beads at the ends have one neighbour.

Parameters:

Name	Type	Description	Default
label	str	The label of the bead or the chain.	required
visual	VisualSpec | None	How the bead is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial first coordinate.	0.0
y	Annotated[float, StateSpec()]	The initial second coordinate.	0.0
theta	Annotated[float, StateSpec(period=2 * pi)]	The initial orientation.	0.0
v0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The self-propulsion speed.	1.0
k	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The stiffness of the bonds.	10.0
D_t	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The translational diffusion coefficient.	0.1
D_r	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rotational diffusion coefficient.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the bead when it sets none.

Examples:

>>> from skcomplex.systems.physics.active_matter import ActiveBeadElement
>>> ActiveBeadElement('bead').k
10.0

Source code in src/skcomplex/systems/physics/active_matter/_polymers.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    y: Annotated[float, StateSpec()] = 0.0,
    theta: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
    v0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 10.0,
    D_t: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    D_r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.theta = theta
    self.v0 = v0
    self.k = k
    self.D_t = D_t
    self.D_r = D_r

ActiveParticleElement(label, visual=POINTS, *, x=0.0, theta=0.0, D_t=0.1, D_r=1.0, omega=0.0)

Bases: Element

Active particle in the plane, whose position diffuses, dx = sqrt(2 D_t) dW, and whose orientation turns.

The orientation turns at the angular velocity omega and diffuses, dtheta = omega dt + sqrt(2 D_r) dW. A SelfPropulsionInteraction adds the drift v (cos theta, sin theta) of the position. With omega = 0 it is an active Brownian particle, and otherwise a chiral active particle.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
theta	Annotated[float, StateSpec(period=2 * pi, role=direction)]	The initial orientation.	0.0
D_t	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The translational diffusion coefficient.	0.1
D_r	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rotational diffusion coefficient.	1.0
omega	Annotated[float, ParameterSpec()]	The angular velocity of a chiral particle.	0.0

Examples:

>>> from skcomplex.systems.physics.common.particles import ActiveParticleElement
>>> ActiveParticleElement('swimmer', D_r=2.0).D_r
2.0

Source code in src/skcomplex/systems/physics/common/particles/_active_particle.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    theta: Annotated[float, StateSpec(period=2 * pi, role='direction')] = 0.0,
    D_t: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    D_r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    omega: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.theta = theta
    self.D_t = D_t
    self.D_r = D_r
    self.omega = omega

ActiveScalarElement(label, visual=None, *, phi=0.0, kappa=1.0, lam=1.0, zeta=4.0, D=0.05)

Bases: Element

Site of Active System B+, a conserved scalar field phi that phase separates with active currents.

The field follows phi_t = laplacian(mu) - div(zeta laplacian(phi) grad phi) + div(sqrt(2 D) eta), with the chemical potential mu = -phi + phi^3 - kappa laplacian(phi) + lam |grad phi|^2. The active terms lam and zeta break the detailed balance of the Cahn-Hilliard equation, and a large zeta reverses Ostwald ripening.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
phi	Annotated[float, StateSpec()]	The initial field.	0.0
kappa	Annotated[float, ParameterSpec(domain=open(0, oo))]	The stiffness of the interfaces.	1.0
lam	Annotated[float, ParameterSpec()]	The active term of the chemical potential.	1.0
zeta	Annotated[float, ParameterSpec()]	The active current.	4.0
D	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the conserved noise.	0.05

Attributes:

Name	Type	Description
default_visual_		The visual of the element when it sets none.

Examples:

>>> from skcomplex.systems.physics.active_matter import ActiveScalarElement
>>> ActiveScalarElement('site').zeta
4.0

Source code in src/skcomplex/systems/physics/active_matter/_system_b.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    phi: Annotated[float, StateSpec()] = 0.0,
    kappa: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    lam: Annotated[float, ParameterSpec()] = 1.0,
    zeta: Annotated[float, ParameterSpec()] = 4.0,
    D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
) -> None:
    super().__init__(label, visual=visual)
    self.phi = phi
    self.kappa = kappa
    self.lam = lam
    self.zeta = zeta
    self.D = D

AngularAlignmentInteraction(label, *, candidate, field, angles=32, frequency=1.0, noise=0.5, state='f')

Bases: EncounterInteraction

Pairs of self-propelled particles leave a collision at their mean angle plus a wrapped normal noise.

A pair of angles meets at the rate 4 d0 v0 |sin((a - b) / 2)|, the frequency of Bertin, Droz, and Grégoire, with d0 the interaction range and v0 the speed. A pair of opposite angles leaves at either of its two mean angles with equal probability.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
candidate	Annotated[str, RoleSpec()]	The path of the field on a phase space with a periodic orientation axis.	required
field	Annotated[str, RoleSpec()]	The path of the field the candidate meets, usually the candidate itself.	required
angles	int	The number of cells of the orientation axis on [0, 2 pi).	32
frequency	float	The product d0 v0 of the interaction range and the speed.	1.0
noise	float	The standard deviation of the noise of the outcome.	0.5
state	str	The name of the state of the density.	'f'

Examples:

>>> from skcomplex.systems.physics.active_matter import AngularAlignmentInteraction
>>> align = AngularAlignmentInteraction('align', candidate='flock', field='flock', frequency=0.5, noise=0.5)
>>> align.outcome
'conservative'

Source code in src/skcomplex/systems/physics/active_matter/_kinetic_flocks.py

def __init__(
    self,
    label: str,
    *,
    candidate: Annotated[str, RoleSpec()],
    field: Annotated[str, RoleSpec()],
    angles: int = 32,
    frequency: float = 1.0,
    noise: float = 0.5,
    state: str = 'f',
) -> None:
    super().__init__(label, candidate=candidate, field=field, rate=self._rate, kernel=self._kernel, state=state)
    self.angles = angles
    self.frequency = frequency
    self.noise = noise

AngularDiffusionInteraction(label, *, candidate, angles=32, tumbling=1.0, noise=0.5, state='f')

Bases: EncounterInteraction

Self-propelled particles change their angle at a rate by a wrapped normal noise, a linear redistribution.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
candidate	Annotated[str, RoleSpec()]	The path of the field on a phase space with a periodic orientation axis.	required
angles	int	The number of cells of the orientation axis on [0, 2 pi).	32
tumbling	float	The rate of the events.	1.0
noise	float	The standard deviation of the noise of an event.	0.5
state	str	The name of the state of the density.	'f'

Examples:

>>> from skcomplex.systems.physics.active_matter import AngularDiffusionInteraction
>>> AngularDiffusionInteraction('tumble', candidate='flock', angles=8, tumbling=1.0, noise=0.5).field

Source code in src/skcomplex/systems/physics/active_matter/_kinetic_flocks.py

def __init__(
    self,
    label: str,
    *,
    candidate: Annotated[str, RoleSpec()],
    angles: int = 32,
    tumbling: float = 1.0,
    noise: float = 0.5,
    state: str = 'f',
) -> None:
    super().__init__(label, candidate=candidate, rate=self._rate, kernel=self._kernel, state=state)
    self.angles = angles
    self.tumbling = tumbling
    self.noise = noise

BinaryCollisionsInteraction(label, *, target, speed=1.0, lam=1.0, sigma=0.3, d0=0.5, substeps=10.0)

Bases: Interaction

Motion, self-diffusion, and aligning binary collisions of self-propelled particles, in steps of 1 / substeps.

In every substep of length tau each particle moves speed * tau along its heading, and with probability 1 - exp(-lam * tau) its heading jumps by a Gaussian angle of standard deviation sigma. Two particles collide when their distance falls below d0 during the substep, and both headings become the mean heading arg(exp(i theta_1) + exp(i theta_2)) plus independent Gaussian angles of standard deviation sigma.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the group of particles.	required
speed	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The self-propulsion speed.	1.0
lam	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of self-diffusion events.	1.0
sigma	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The standard deviation of the angular noise of both kinds of events.	0.3
d0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The collision distance.	0.5
substeps	Annotated[float, ParameterSpec(domain=open(0, oo))]	The number of substeps per unit of time.	10.0

Examples:

>>> from skcomplex.systems.physics.active_matter import BinaryCollisionsInteraction
>>> BinaryCollisionsInteraction('collisions', target='birds').d0
0.5

Source code in src/skcomplex/systems/physics/active_matter/_collisions.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    speed: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    lam: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
    d0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    substeps: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
) -> None:
    super().__init__(label)
    self.target = target
    self.speed = speed
    self.lam = lam
    self.sigma = sigma
    self.d0 = d0
    self.substeps = substeps

fly(system)

Move the particles for one unit of time, with self-diffusion and collisions in every substep.

Parameters:

Name	Type	Description	Default
system	SystemView	The system.	required

Source code in src/skcomplex/systems/physics/active_matter/_collisions.py

def fly(self, system: SystemView) -> None:
    """Move the particles for one unit of time, with self-diffusion and collisions in every substep.

    Args:
        system:
            The system.
    """
    agents = system.agents(self.target, FlockerView)
    size = np.asarray(system.dimensions, dtype=float)
    positions = np.array([agent.location for agent in agents], dtype=float)
    headings = np.array([agent.theta for agent in agents], dtype=float)
    speed, lam, sigma, d0 = (system.parameter(name) for name in ('speed', 'lam', 'sigma', 'd0'))
    substeps = max(1, round(system.parameter('substeps')))
    tau, rng = 1.0 / substeps, system.rng
    touching = KDTree(positions, boxsize=size).query_pairs(d0)
    for _ in range(substeps):
        positions = np.mod(positions + speed * tau * np.column_stack([np.cos(headings), np.sin(headings)]), size)
        kicked = rng.random(len(headings)) < -math.expm1(-lam * tau)
        headings[kicked] += rng.normal(0.0, sigma, int(kicked.sum()))
        now = KDTree(positions, boxsize=size).query_pairs(d0)
        BinaryCollisionsInteraction._collide(headings, sorted(now - touching), sigma, rng)
        touching = now
    for agent, position, heading in zip(agents, positions, np.mod(headings, 2 * math.pi), strict=True):
        agent.move_to(position)
        agent.theta = float(heading)

ChemotaxisElement(label, visual=None, *, rho=1.0, c=1.0, D=1.0, chi=5.0, Dc=1.0, k=1.0)

Bases: Element

Site of the Keller-Segel system, cells of density rho that climb the gradient of a chemical c they secrete.

The cells diffuse and drift up the chemical, rho_t = D laplacian(rho) - chi div(rho grad c), and the chemical diffuses, is secreted by the cells, and decays, c_t = Dc laplacian(c) + rho - k c. Above a critical mass the cells aggregate and the density blows up in finite time.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
rho	Annotated[float, StateSpec()]	The initial density of cells.	1.0
c	Annotated[float, StateSpec()]	The initial concentration of the chemical.	1.0
D	Annotated[float, ParameterSpec(domain=open(0, oo))]	The diffusion coefficient of the cells.	1.0
chi	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The chemotactic sensitivity.	5.0
Dc	Annotated[float, ParameterSpec(domain=open(0, oo))]	The diffusion coefficient of the chemical.	1.0
k	Annotated[float, ParameterSpec(domain=open(0, oo))]	The decay rate of the chemical.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the element when it sets none.

Examples:

>>> from skcomplex.systems.physics.active_matter import ChemotaxisElement
>>> ChemotaxisElement('site').chi
5.0

Source code in src/skcomplex/systems/physics/active_matter/_chemotaxis.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    rho: Annotated[float, StateSpec()] = 1.0,
    c: Annotated[float, StateSpec()] = 1.0,
    D: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    chi: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5.0,
    Dc: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.rho = rho
    self.c = c
    self.D = D
    self.chi = chi
    self.Dc = Dc
    self.k = k

CohesiveParticleElement(label, visual=None, *, x=0.0, theta=0.0, eta=0.1, alpha=1.0, beta=1.0, r_c=CLOSE, r_e=0.5, r_a=0.8)

Bases: Element

Particle of a Grégoire-Chaté flock, which aligns, repels close neighbours, and attracts distant ones.

The new heading is the direction of alpha sum_j exp(i theta_j) + beta sum_j f(r_ij) e_ij + eta n_i xi over the n_i particles within the unit distance, where e_ij points from i to j, xi is a random unit vector, and the force f is a strong repulsion below r_c, zero at r_e, and attractive up to r_a. A propulsion step moves the particle along its heading.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the particle is drawn, or None for default_visual_.	None
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
theta	Annotated[float, StateSpec(period=2 * pi, role=direction)]	The initial heading.	0.0
eta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the vectorial noise.	0.1
alpha	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the alignment.	1.0
beta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the cohesion.	1.0
r_c	Annotated[float, ParameterSpec(domain=open(0, oo))]	The distance below which the neighbours repel strongly.	CLOSE
r_e	Annotated[float, ParameterSpec(domain=open(0, oo))]	The preferred distance.	0.5
r_a	Annotated[float, ParameterSpec(domain=open(0, oo))]	The distance of the strongest attraction.	0.8

Attributes:

Name	Type	Description
default_visual_		The visual of the particle when it sets none.

Examples:

>>> from skcomplex.systems.physics.active_matter import CohesiveParticleElement
>>> CohesiveParticleElement('bird').r_e
0.5

Source code in src/skcomplex/systems/physics/active_matter/_flocks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    theta: Annotated[float, StateSpec(period=2 * pi, role='direction')] = 0.0,
    eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    r_c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = CLOSE,
    r_e: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    r_a: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.8,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.theta = theta
    self.eta = eta
    self.alpha = alpha
    self.beta = beta
    self.r_c = r_c
    self.r_e = r_e
    self.r_a = r_a

ContractileGelSiteElement(label, visual=None, *, c=1.0, u=0.0, D=0.1, gamma=1.0, eta=1.0, E=0.1, T0=4.0, c_sat=1.0)

Bases: Element

Site of a one-dimensional active viscoelastic gel, the density of contractile cells and the gel displacement.

The site holds the density c of the cells and the displacement u of a Kelvin-Voigt gel of stiffness E and viscosity eta, with friction gamma against a rigid substrate. The cells diffuse with the coefficient D and exert the active stress T0 c / (c + c_sat).

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the gel.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
c	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial density of the cells.	1.0
u	Annotated[float, StateSpec()]	The initial displacement of the gel.	0.0
D	Annotated[float, ParameterSpec(domain=open(0, oo))]	The diffusion coefficient of the cells.	0.1
gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The friction of the gel with the substrate.	1.0
eta	Annotated[float, ParameterSpec(domain=open(0, oo))]	The viscosity of the gel.	1.0
E	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The stiffness of the gel.	0.1
T0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The largest active stress of the cells.	4.0
c_sat	Annotated[float, ParameterSpec(domain=open(0, oo))]	The density at which the active stress is half its largest value.	1.0

Examples:

>>> from skcomplex.systems.physics.active_matter import ContractileGelSiteElement
>>> ContractileGelSiteElement('site').T0
4.0

Source code in src/skcomplex/systems/physics/active_matter/_gels.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    c: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    u: Annotated[float, StateSpec()] = 0.0,
    D: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    eta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    E: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    T0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
    c_sat: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.c = c
    self.u = u
    self.D = D
    self.gamma = gamma
    self.eta = eta
    self.E = E
    self.T0 = T0
    self.c_sat = c_sat

CuckerSmaleAlignmentInteraction(label, first, second, cutoff=1.0, position='x', *, K=1.0, beta=0.25, period=None, n_agents=1)

Bases: PairCouplingInteraction

Cucker-Smale alignment, x_i'' += K / N psi(r_ij) (x_j' - x_i') with psi(r) = (1 + r^2)^(-beta).

Every pair of members within the cutoff pulls their velocities, the first derivatives of their positions of order 2, towards each other with a weight that decays with their distance.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
first	str	The path of the group or particle of the first members.	required
second	str	The path of the group or particle of the second members.	required
cutoff	float	The distance beyond which pairs do not align, large for the all-to-all system.	1.0
position	str	The name of the vector state of order 2 that holds the positions of the members.	'x'
K	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the alignment.	1.0
beta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The decay exponent of the weight.	0.25
period	float | None	The side of a periodic box, for minimum image distances, or None for the period of the space.	None
n_agents	int	The number of agents N that normalizes the alignment.	1

Examples:

>>> from skcomplex.systems.physics.active_matter import CuckerSmaleAlignmentInteraction
>>> CuckerSmaleAlignmentInteraction('align', first='birds', second='birds', beta=0.25).beta
0.25

Source code in src/skcomplex/systems/physics/active_matter/_collective.py

def __init__(
    self,
    label: str,
    first: str,
    second: str,
    cutoff: float = 1.0,
    position: str = 'x',
    *,
    K: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.25,
    period: float | None = None,
    n_agents: int = 1,
) -> None:
    super().__init__(
        label,
        first,
        second,
        cutoff=cutoff,
        position=position,
        write=position,
        function=self._align,
        period=period,
    )
    self.K = K
    self.beta = beta
    self.n_agents = n_agents

DeformableParticleElement(label, visual=POINTS, *, x=0.0, v=0.0, S=0.0, gamma=1.0, a=1.0, b=1.0, kappa=0.5, noise=0.0)

Bases: Element

Deformable self-propelled particle, whose velocity and traceless deformation tensor drive each other.

The particle moves with its velocity v, which follows v' = gamma v - |v|² v - a S v, and its deformation, the traceless symmetric tensor S, follows S' = -kappa S + b (v v - |v|² I / 2), the model of Ohta and Ohkuma. The motion stretches the particle along the velocity for b > 0, and the deformation turns the velocity back. Noise of amplitude noise may kick the velocity.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the particle is drawn.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
v	Annotated[ArrayLike, StateSpec(vector=True, role=velocity)]	The initial velocity.	0.0
S	Annotated[ArrayLike, StateSpec(tensor=traceless)]	The initial deformation, given by its two independent components.	0.0
gamma	Annotated[float, ParameterSpec()]	The distance from the onset of self-propulsion.	1.0
a	Annotated[float, ParameterSpec()]	The coupling of the velocity to the deformation.	1.0
b	Annotated[float, ParameterSpec()]	The coupling of the deformation to the velocity.	1.0
kappa	Annotated[float, ParameterSpec(domain=open(0, oo))]	The relaxation rate of the deformation.	0.5
noise	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the noise of the velocity.	0.0

Examples:

>>> from skcomplex.systems.physics.active_matter import DeformableParticleElement
>>> DeformableParticleElement('cell').kappa
0.5

Source code in src/skcomplex/systems/physics/active_matter/_deformable.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    v: Annotated[npt.ArrayLike, StateSpec(vector=True, role='velocity')] = 0.0,
    S: Annotated[npt.ArrayLike, StateSpec(tensor='traceless')] = 0.0,
    gamma: Annotated[float, ParameterSpec()] = 1.0,
    a: Annotated[float, ParameterSpec()] = 1.0,
    b: Annotated[float, ParameterSpec()] = 1.0,
    kappa: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    noise: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.v = v
    self.S = S
    self.gamma = gamma
    self.a = a
    self.b = b
    self.kappa = kappa
    self.noise = noise

DelayedVicsekParticleElement(label, visual=None, *, x=0.0, theta=0.0, eta=0.1)

Bases: Element

Particle of a delayed Vicsek flock, which keeps its past headings for the delayed or predicted alignment.

The heading keeps its last 25 values, theta_lag1 to theta_lag25, which a delayed alignment reads, and which hold the initial heading before the start. A propulsion step moves the particle along its current heading.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the particle is drawn, or None for default_visual_.	None
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
theta	Annotated[float, StateSpec(period=2 * pi, role=direction, lags=MEMORY)]	The initial heading, also its value before the start.	0.0
eta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the noise, as a fraction of pi.	0.1

Attributes:

Name	Type	Description
default_visual_		The visual of the particle when it sets none.

Examples:

>>> from skcomplex.systems.physics.active_matter import DelayedVicsekParticleElement
>>> DelayedVicsekParticleElement('bird').eta
0.1

Source code in src/skcomplex/systems/physics/active_matter/_flocks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    theta: Annotated[float, StateSpec(period=2 * pi, role='direction', lags=MEMORY)] = 0.0,
    eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.theta = theta
    self.eta = eta

DepotParticleElement(label, visual=POINTS, *, x=0.0, v=0.0, e=0.0, gamma=20.0, q=10.0, c=1.0, d=10.0, D=0.01)

Bases: Element

Brownian particle pumped by an internal energy depot, as in the system of Schweitzer, Ebeling, and Tilch.

The depot takes energy at the rate q, loses it at the rate c e, and converts it into motion at the rate d e |v|^2, so de/dt = q - c e - d e |v|^2 and dv = (-gamma v + d e v) dt + sqrt(2 D) dW. Above q d = c gamma the particle moves at the speed sqrt(q / gamma - c / d).

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
v	Annotated[ArrayLike, StateSpec(vector=True, role=velocity)]	The initial velocity.	0.0
e	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial energy of the depot.	0.0
gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The friction coefficient.	20.0
q	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of energy intake.	10.0
c	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of internal energy loss.	1.0
d	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of conversion of depot energy into motion.	10.0
D	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The intensity of the velocity noise.	0.01

Examples:

>>> from skcomplex.systems.physics.active_matter import DepotParticleElement
>>> DepotParticleElement('walker').q
10.0

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    v: Annotated[npt.ArrayLike, StateSpec(vector=True, role='velocity')] = 0.0,
    e: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 20.0,
    q: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 10.0,
    c: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    d: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 10.0,
    D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.v = v
    self.e = e
    self.gamma = gamma
    self.q = q
    self.c = c
    self.d = d
    self.D = D

FlockerElement(label, visual=VisualSpec('agents', {'color_by': 'theta', 'color_range': (0.0, 2 * math.pi)}), *, theta=0.0)

Bases: Element

Self-propelled point particle with a heading theta, in a periodic square.

Parameters:

Name	Type	Description	Default
label	str	The label of the agent or group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	VisualSpec('agents', {'color_by': 'theta', 'color_range': (0.0, 2 * pi)})
theta	Annotated[float, StateSpec(role=direction)]	The initial heading.	0.0

Examples:

>>> from skcomplex.systems.physics.active_matter import FlockerElement
>>> FlockerElement('bird', theta=1.0).theta
1.0

Source code in src/skcomplex/systems/physics/active_matter/_collisions.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('agents', {'color_by': 'theta', 'color_range': (0.0, 2 * math.pi)}),
    *,
    theta: Annotated[float, StateSpec(role='direction')] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.theta = theta

FlockerView(agent, stepper)

Bases: AgentView

Flocker as rules see it, with its heading.

Parameters:

Name	Type	Description	Default
agent	Agent	The Mesa agent.	required
stepper	MesaSolver	The stepper that runs the system.	required

Source code in src/skcomplex/core/engine/solvers/_mesa.py

def __init__(self, agent: mesa.Agent, stepper: 'MesaSolver') -> None:
    self._agent = agent
    self._stepper = stepper

FlockingAgentElement(label, x=0.0, x_d1=0.0, visual=None, *, active=ACTIVE_AGENT)

Bases: PointElement

Agent of a crowd, a point with a mark that tells whether it is active or passive.

Parameters:

Name	Type	Description	Default
label	str	The label of the agent or the group.	required
x	Annotated[ArrayLike, StateSpec(vector=True, order=2, role=position)]	The initial position.	0.0
x_d1	ArrayLike	The initial velocity.	0.0
visual	VisualSpec | None	How the agent is drawn, or None for default_visual_.	None
active	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([PASSIVE_AGENT, ACTIVE_AGENT]))]	The mark, 1 for an active agent or 0 for a passive one.	ACTIVE_AGENT

Examples:

>>> from skcomplex.systems.physics.active_matter import FlockingAgentElement
>>> FlockingAgentElement('agent', active=0.0).active
0.0

Source code in src/skcomplex/systems/physics/active_matter/_active_passive.py

def __init__(
    self,
    label: str,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, order=2, role='position')] = 0.0,
    x_d1: npt.ArrayLike = 0.0,
    visual: VisualSpec | None = None,
    *,
    active: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([PASSIVE_AGENT, ACTIVE_AGENT]))] = ACTIVE_AGENT,
) -> None:
    super().__init__(label, x=x, x_d1=x_d1, visual=visual)
    self.active = active

HeadingAlignmentInteraction(label, *, target, nematic=False, neighborhood=None, latency=0, memory=5)

Bases: CouplingInteraction

Vicsek alignment, the turn of every particle to the circular mean heading of its neighbours plus noise.

Every step the heading theta of every particle becomes the circular mean of the headings within the unit distance, itself included, plus a uniform angle in [-eta pi, eta pi] with the parameter eta of the particle. With nematic, every heading is flipped to the side nearer the particle's own heading first, as for particles with an axis rather than a direction. It is the map form of coupling with the circular mean over a radius. Another neighbourhood, such as the k nearest particles or a vision cone, replaces the unit disk.

With a positive latency h, every particle reads the headings of its neighbours h steps ago. With a negative latency it reads their headings predicted |h| steps ahead, by a least-squares line through their last memory headings, unwrapped by taking every change between two steps as the smallest angle. A latent alignment needs particles that keep enough past headings, such as DelayedVicsekParticleElement.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the group of particles.	required
nematic	bool	Whether every heading is flipped to the side nearer the particle's own heading.	False
neighborhood	BaseNeighborhood | None	The neighbourhood of every particle, itself included, or None for the unit disk.	None
latency	int	The number of steps by which the headings of the neighbours are read late, or predicted ahead when negative.	0
memory	int	The number of past headings a prediction fits, at least 2.	5

Examples:

>>> from skcomplex.systems.physics.active_matter import HeadingAlignmentInteraction
>>> HeadingAlignmentInteraction('alignment', target='flock', nematic=True).reduce
'circular-mean'

Source code in src/skcomplex/systems/physics/active_matter/_flocks.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    nematic: bool = False,
    neighborhood: BaseNeighborhood | None = None,
    latency: int = 0,
    memory: int = 5,
) -> None:
    super().__init__(
        label,
        target=target,
        read='theta',
        write='theta',
        function=self._get_heading,
        reduce='circular-mean',
        weights=self._get_side,
        form='map',
        combine=self._get_noisy,
        neighborhood=RadiusNeighborhood(1.0, center=True) if neighborhood is None else neighborhood,
    )
    self.nematic = nematic
    self.latency = latency
    self.memory = memory

InertialActiveParticleElement(label, visual=POINTS, *, x=0.0, x_d1=0.0, theta=0.0, theta_d1=0.0, speed=1.0, mass=0.5, moment=0.1, D_t=0.1, D_r=1.0)

Bases: Element

Active Brownian particle with mass and moment of inertia, whose velocity and spin relax with friction.

With unit friction, the position follows m x'' = -x' + v0 n + sqrt(2 D_t) xi along the orientation n = (cos theta, sin theta), and the orientation follows J theta'' = -theta' + sqrt(2 D_r) xi_r. The velocity relaxes over the time m and the angular velocity over the time J, and both vanish to give the overdamped active Brownian particle.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the particle is drawn.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, order=2, role=position)]	The initial position.	0.0
x_d1	ArrayLike	The initial velocity.	0.0
theta	Annotated[float, StateSpec(order=2, role=direction)]	The initial orientation.	0.0
theta_d1	float	The initial angular velocity.	0.0
speed	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The self-propulsion speed v0.	1.0
mass	Annotated[float, ParameterSpec(domain=open(0, oo))]	The mass m, in units of the friction.	0.5
moment	Annotated[float, ParameterSpec(domain=open(0, oo))]	The moment of inertia J, in units of the rotational friction.	0.1
D_t	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The translational diffusion coefficient.	0.1
D_r	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rotational diffusion coefficient.	1.0

Examples:

>>> from skcomplex.systems.physics.active_matter import InertialActiveParticleElement
>>> InertialActiveParticleElement('swimmer').mass
0.5

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, order=2, role='position')] = 0.0,
    x_d1: npt.ArrayLike = 0.0,
    theta: Annotated[float, StateSpec(order=2, role='direction')] = 0.0,
    theta_d1: float = 0.0,
    speed: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    mass: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    moment: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
    D_t: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    D_r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.x_d1 = x_d1
    self.theta = theta
    self.theta_d1 = theta_d1
    self.speed = speed
    self.mass = mass
    self.moment = moment
    self.D_t = D_t
    self.D_r = D_r

LandauDeGennesEnergyInteraction(label, *, target, state='Q', a=-1.0, c=2.0, K=1.0, gamma=1.0, quadratic=None)

Bases: GradientInteraction

Gradient flow of a symmetric traceless tensor Q in the Landau-de Gennes free energy.

The free energy is F = integral a tr(Q^2) / 2 + c tr(Q^2)^2 / 4 + K |grad Q|^2 / 2. Its negative functional derivative among symmetric traceless tensors is the molecular field H = -a Q - c tr(Q^2) Q + K laplacian(Q), and the tensor relaxes as Q_t = H / gamma. A uniform tensor S (n n - I / 2) in two dimensions with a < 0 relaxes to the scalar order S = sqrt(-2 a / c).

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the field.	required
state	str	The name of the state of the field that holds the tensor.	'Q'
a	Annotated[float, ParameterSpec()]	The coefficient of tr(Q^2) / 2, negative in the ordered phase.	-1.0
c	Annotated[float, ParameterSpec(domain=open(0, oo))]	The coefficient of tr(Q^2)^2 / 4.	2.0
K	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The elastic constant.	1.0
gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The rotational viscosity, which divides the molecular field.	1.0
quadratic	Callable[[SimpleNamespace], object] | None	A function of the namespace that returns the coefficient of tr(Q^2) / 2 instead of a, such as an expression of a density, or None for a.	None

Examples:

>>> from sage.all__sagemath_symbolics import SR, matrix, var
>>> from skcomplex.systems.physics.active_matter import LandauDeGennesEnergyInteraction
>>> q, a, c = var('q, a, c')
>>> tensor = matrix(SR, [[q, 0], [0, -q]])
>>> LandauDeGennesEnergyInteraction._molecular_field(tensor, a, c, 0)[0, 0].factor()
-(2*c*q^2 + a)*q

Source code in src/skcomplex/systems/physics/active_matter/_nematics.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    state: str = 'Q',
    a: Annotated[float, ParameterSpec()] = -1.0,
    c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    quadratic: Callable[[SimpleNamespace], object] | None = None,
) -> None:
    super().__init__(label, target=target, states=state, variational=True)
    self.state = state
    self.a = a
    self.c = c
    self.K = K
    self.gamma = gamma
    self.quadratic = quadratic

MotileCellElement(label, visual=None, *, x=0.0, theta=0.0, v0=1.0, mu=1.0, tau=1.0, eta=0.6, R_eq=5.0 / 6.0, R_0=1.0, F_rep=30.0, F_adh=0.75)

Bases: Element

Motile cell of a tissue, with its position, the heading of its motility, and its force and motility parameters.

The position follows x' = v0 n + mu F, where n is the unit vector of the heading and F the sum of the pair forces. A pair force is F_rep (d - R_eq) / R_eq towards the neighbour below the distance R_eq, and F_adh (d - R_eq) / (R_0 - R_eq) towards it up to the distance R_0, the cell diameter.

Parameters:

Name	Type	Description	Default
label	str	The label of the cell or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
theta	Annotated[float, StateSpec(period=2 * pi, role=direction)]	The initial heading.	0.0
v0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The speed of the motility.	1.0
mu	Annotated[float, ParameterSpec(domain=open(0, oo))]	The mobility.	1.0
tau	Annotated[float, ParameterSpec(domain=open(0, oo))]	The relaxation time of the heading.	1.0
eta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the angular noise, whose variance per unit time is eta^2 / 12.	0.6
R_eq	Annotated[float, ParameterSpec(domain=open(0, oo))]	The equilibrium distance of a pair.	5.0 / 6.0
R_0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The distance beyond which a pair exerts no force.	1.0
F_rep	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The largest repulsive force.	30.0
F_adh	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The largest adhesive force.	0.75

Examples:

>>> from skcomplex.systems.physics.active_matter import MotileCellElement
>>> MotileCellElement('cell').F_rep
30.0

Source code in src/skcomplex/systems/physics/active_matter/_cells.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    theta: Annotated[float, StateSpec(period=2 * pi, role='direction')] = 0.0,
    v0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    mu: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.6,
    R_eq: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 5.0 / 6.0,
    R_0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    F_rep: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 30.0,
    F_adh: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.75,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.theta = theta
    self.v0 = v0
    self.mu = mu
    self.tau = tau
    self.eta = eta
    self.R_eq = R_eq
    self.R_0 = R_0
    self.F_rep = F_rep
    self.F_adh = F_adh

NematicDensitySiteElement(label, visual=None, *, rho=1.0, Q=0.0, beta=0.5)

Bases: Element

Site of a dry active nematic, whose alignment tensor grows with the traceless curvature of the density.

The rate of the tensor from the density is beta (partial_i partial_j rho - delta_ij laplacian(rho) / 2). The Landau-de Gennes energy and the current of the density are interactions of the system.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
rho	Annotated[ArrayLike, StateSpec()]	The initial density.	1.0
Q	Annotated[ArrayLike, StateSpec(tensor=traceless)]	The initial alignment tensor, its independent components Q_xx and Q_xy.	0.0
beta	Annotated[float, ParameterSpec()]	The coupling of the tensor to the curvature of the density.	0.5

Attributes:

Name	Type	Description
default_visual_		The visual of the element when it sets none.

Examples:

>>> from skcomplex.systems.physics.active_matter import NematicDensitySiteElement
>>> NematicDensitySiteElement('site').beta
0.5

Source code in src/skcomplex/systems/physics/active_matter/_nematics.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    rho: Annotated[npt.ArrayLike, StateSpec()] = 1.0,
    Q: Annotated[npt.ArrayLike, StateSpec(tensor='traceless')] = 0.0,
    beta: Annotated[float, ParameterSpec()] = 0.5,
) -> None:
    super().__init__(label, visual=visual)
    self.rho = rho
    self.Q = Q
    self.beta = beta

NematicGelSiteElement(label, visual=None, *, Q=0.0, a=-1.0, c=2.0, K=1.0, xi=0.0, gamma=1.0, eta=1.0, zeta=0.5)

Bases: Element

Site of an active nematic gel, a symmetric traceless tensor Q carried by an overdamped Stokes flow.

The tensor follows Q_t = -corotational(Q, u, xi) + H / gamma, with the molecular field H of the Landau-de Gennes energy. The velocity solves eta laplacian(u) - grad(P) + div(sigma) = 0, div(u) = 0, with the stress sigma = -xi H + Q H - H Q - zeta Q. An activity zeta > 0 is extensile and an activity zeta < 0 contractile.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
Q	Annotated[ArrayLike, StateSpec(tensor=traceless)]	The initial tensor, its independent components Q_xx and Q_xy.	0.0
a	Annotated[float, ParameterSpec()]	The coefficient of tr(Q^2) / 2 of the energy, negative in the ordered phase.	-1.0
c	Annotated[float, ParameterSpec(domain=open(0, oo))]	The coefficient of tr(Q^2)^2 / 4 of the energy.	2.0
K	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The elastic constant.	1.0
xi	Annotated[float, ParameterSpec()]	The flow-alignment parameter.	0.0
gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The rotational viscosity.	1.0
eta	Annotated[float, ParameterSpec(domain=open(0, oo))]	The viscosity of the flow.	1.0
zeta	Annotated[float, ParameterSpec()]	The activity.	0.5

Attributes:

Name	Type	Description
default_visual_		The visual of the element when it sets none.

Examples:

>>> from skcomplex.systems.physics.active_matter import NematicGelSiteElement
>>> NematicGelSiteElement('site').zeta
0.5

Source code in src/skcomplex/systems/physics/active_matter/_nematics.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    Q: Annotated[npt.ArrayLike, StateSpec(tensor='traceless')] = 0.0,
    a: Annotated[float, ParameterSpec()] = -1.0,
    c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    xi: Annotated[float, ParameterSpec()] = 0.0,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    eta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    zeta: Annotated[float, ParameterSpec()] = 0.5,
) -> None:
    super().__init__(label, visual=visual)
    self.Q = Q
    self.a = a
    self.c = c
    self.K = K
    self.xi = xi
    self.gamma = gamma
    self.eta = eta
    self.zeta = zeta

NematicParticleElement(label, visual=None, *, x=0.0, theta=0.0, eta=0.1)

Bases: Element

Polar particle with nematic alignment, which aligns with the axis of every neighbour on its nearer side.

HeadingAlignmentInteraction with nematic=True sets the new heading to the circular mean of the headings within the unit distance, itself included, each flipped to the side nearer the particle's own heading, plus a uniform angle in [-eta pi, eta pi]. A propulsion step moves the particle along its heading.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the particle is drawn, or None for default_visual_.	None
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
theta	Annotated[float, StateSpec(period=2 * pi, role=direction)]	The initial heading.	0.0
eta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the noise, as a fraction of pi.	0.1

Attributes:

Name	Type	Description
default_visual_		The visual of the particle when it sets none.

Examples:

>>> from skcomplex.systems.physics.active_matter import NematicParticleElement
>>> NematicParticleElement('rod', eta=0.2).eta
0.2

Source code in src/skcomplex/systems/physics/active_matter/_flocks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    theta: Annotated[float, StateSpec(period=2 * pi, role='direction')] = 0.0,
    eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.theta = theta
    self.eta = eta

NonreciprocalMoverElement(label, visual=SPECIES, *, x=0.0, theta=0.0, species=0.0, eta=0.3, J_aa=1.0, J_ab=1.0, J_ba=1.0, J_bb=1.0, radius=1.0)

Bases: Element

Vicsek mover of species A (0) or B (1), which aligns with the neighbours within radius by species couplings.

Every step the heading becomes the direction of sum_j J(species_i, species_j) exp(i theta_j) over the movers within the radius, itself included, plus uniform noise of width eta. A propulsion step moves the mover along it. J_aa and J_bb align a species with itself, and J_ab and J_ba align A with B and B with A, so J_ab != J_ba makes the alignment nonreciprocal.

Parameters:

Name	Type	Description	Default
label	str	The label of the mover.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	SPECIES
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The position.	0.0
theta	Annotated[float, StateSpec(period=2 * pi, role=direction)]	The heading.	0.0
species	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The species, 0 for A or 1 for B.	0.0
eta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The width of the uniform noise of the heading.	0.3
J_aa	Annotated[float, ParameterSpec()]	The alignment of A with A.	1.0
J_ab	Annotated[float, ParameterSpec()]	The alignment of A with B.	1.0
J_ba	Annotated[float, ParameterSpec()]	The alignment of B with A.	1.0
J_bb	Annotated[float, ParameterSpec()]	The alignment of B with B.	1.0
radius	float	The radius of the neighbourhood.	1.0

Examples:

>>> from skcomplex.systems.physics.active_matter import NonreciprocalMoverElement
>>> NonreciprocalMoverElement('mover', J_ba=-1.0).J_ba
-1.0

Source code in src/skcomplex/systems/physics/active_matter/_nonreciprocal.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = SPECIES,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    theta: Annotated[float, StateSpec(period=2 * pi, role='direction')] = 0.0,
    species: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
    eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
    J_aa: Annotated[float, ParameterSpec()] = 1.0,
    J_ab: Annotated[float, ParameterSpec()] = 1.0,
    J_ba: Annotated[float, ParameterSpec()] = 1.0,
    J_bb: Annotated[float, ParameterSpec()] = 1.0,
    radius: float = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.theta = theta
    self.species = species
    self.eta = eta
    self.J_aa = J_aa
    self.J_ab = J_ab
    self.J_ba = J_ba
    self.J_bb = J_bb
    self.radius = radius

OrnsteinUhlenbeckParticleElement(label, visual=POINTS, *, x=0.0, u=0.0, tau=1.0, D=1.0, D_t=0.0)

Bases: Element

Active Ornstein-Uhlenbeck particle, whose propulsion velocity u is an Ornstein-Uhlenbeck process.

The laws are dx = u dt + sqrt(2 D_t) dW and tau du = -u dt + sqrt(2 D) dW, so the propulsion keeps its direction for the persistence time tau and the particle diffuses with D + D_t at long times.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
u	Annotated[ArrayLike, StateSpec(vector=True)]	The initial propulsion velocity.	0.0
tau	Annotated[float, ParameterSpec(domain=open(0, oo))]	The persistence time.	1.0
D	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The active diffusion coefficient.	1.0
D_t	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The translational diffusion coefficient.	0.0

Examples:

>>> from skcomplex.systems.physics.active_matter import OrnsteinUhlenbeckParticleElement
>>> OrnsteinUhlenbeckParticleElement('swimmer').tau
1.0

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    u: Annotated[npt.ArrayLike, StateSpec(vector=True)] = 0.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    D_t: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.u = u
    self.tau = tau
    self.D = D
    self.D_t = D_t

PhoreticAttractionInteraction(label, first, second, cutoff=1.0, position='x', *, v0=20.0, R=1.0, period=None)

Bases: PairCouplingInteraction

Phoretic attraction between active colloids, the pair potential -v0 R^2 / r of an overdamped position.

Every member of a pair drifts towards the other at the speed v0 R^2 / r^2.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
first	str	The path of the group or particle of the first members.	required
second	str	The path of the group or particle of the second members.	required
cutoff	float	The distance beyond which the attraction is neglected.	1.0
position	str	The name of the vector state that holds the positions of the members.	'x'
v0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The drift speed at the distance R.	20.0
R	Annotated[float, ParameterSpec(domain=open(0, oo))]	The particle radius.	1.0
period	float | None	The side of a periodic box, for minimum image distances, or None for the period of the space.	None

Examples:

>>> from skcomplex.systems.physics.active_matter import PhoreticAttractionInteraction
>>> PhoreticAttractionInteraction('phoresis', first='colloids', second='colloids', v0=20.0).v0
20.0

Source code in src/skcomplex/systems/physics/active_matter/_collective.py

def __init__(
    self,
    label: str,
    first: str,
    second: str,
    cutoff: float = 1.0,
    position: str = 'x',
    *,
    v0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 20.0,
    R: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    period: float | None = None,
) -> None:
    super().__init__(
        label,
        first,
        second,
        cutoff=cutoff,
        position=position,
        write=position,
        potential=self._attraction,
        period=period,
    )
    self.v0 = v0
    self.R = R

PolarAlignmentInteraction(label, first, second, cutoff=1.0, position='x', *, gamma=1.0, period=None)

Bases: PairCouplingInteraction

Alignment torque between two particles closer than the cutoff, theta_i' += gamma sin(theta_j - theta_i).

Summed over every neighbour j within the interaction range d0, the cutoff, the torque turns each particle towards the headings of its neighbours.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
first	str	The path of the group or particle of the first members.	required
second	str	The path of the group or particle of the second members.	required
cutoff	float	The interaction range d0.	1.0
position	str	The name of the vector state that holds the positions of the members.	'x'
gamma	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the alignment.	1.0
period	float | None	The side of a periodic box, for minimum image distances, or None for the period of the space.	None

Examples:

>>> from skcomplex.systems.physics.active_matter import PolarAlignmentInteraction
>>> PolarAlignmentInteraction('align', first='swimmers', second='swimmers', gamma=0.5).gamma
0.5

Source code in src/skcomplex/systems/physics/active_matter/_alignment.py

def __init__(
    self,
    label: str,
    first: str,
    second: str,
    cutoff: float = 1.0,
    position: str = 'x',
    *,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    period: float | None = None,
) -> None:
    super().__init__(
        label,
        first,
        second,
        cutoff=cutoff,
        position=position,
        write='theta',
        function=self._torque,
        period=period,
    )
    self.gamma = gamma

PolarFluidElement(label, visual=None, *, rho=1.0, vx=0.0, vy=0.0, a=1.0, rho_c=0.5, beta=1.0, lam=1.0, sigma=0.5, D=0.5, noise=0.0)

Bases: Element

Site of the Toner-Tu polar active fluid, with a density rho and a velocity (vx, vy).

The density is conserved, drho/dt + div(rho v) = 0, and the velocity follows dv/dt + lam (v . grad) v = a (rho - rho_c) v - beta |v|^2 v - sigma grad(rho) + D laplacian(v) with additive noise of intensity noise. Above the critical density the fluid moves with the speed sqrt(a (rho - rho_c) / beta).

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
rho	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial density.	1.0
vx	Annotated[float, StateSpec()]	The initial velocity along the first axis.	0.0
vy	Annotated[float, StateSpec()]	The initial velocity along the second axis.	0.0
a	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of the ordering instability per unit density.	1.0
rho_c	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The critical density of the ordering transition.	0.5
beta	Annotated[float, ParameterSpec(domain=open(0, oo))]	The saturation of the speed.	1.0
lam	Annotated[float, ParameterSpec()]	The strength of the self-advection.	1.0
sigma	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The pressure per unit density.	0.5
D	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The viscosity of the velocity.	0.5
noise	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The intensity of the velocity noise.	0.0

Examples:

>>> from skcomplex.systems.physics.active_matter import PolarFluidElement
>>> PolarFluidElement('site').rho_c
0.5

Source code in src/skcomplex/systems/physics/active_matter/_collective.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    rho: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    vx: Annotated[float, StateSpec()] = 0.0,
    vy: Annotated[float, StateSpec()] = 0.0,
    a: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    rho_c: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    beta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    lam: Annotated[float, ParameterSpec()] = 1.0,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    noise: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.rho = rho
    self.vx = vx
    self.vy = vy
    self.a = a
    self.rho_c = rho_c
    self.beta = beta
    self.lam = lam
    self.sigma = sigma
    self.D = D
    self.noise = noise

PolarGelSiteElement(label, visual=None, *, p=0.0, a=-1.0, c=1.0, K=1.0, xi=0.0, gamma=1.0, eta=1.0, zeta=-0.5)

Bases: Element

Site of an active polar gel, a polarisation p carried by an overdamped Stokes flow.

The polarisation follows p_t = -corotational(p, u, xi) + h / gamma, with the molecular field h = -a p - c |p|^2 p + K laplacian(p). The velocity solves eta laplacian(u) - grad(P) + div(sigma) = 0, div(u) = 0, with the stress sigma = -xi (p h + h p) / 2 + (p h - h p) / 2 - zeta p p, as in Kruse et al.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
p	Annotated[ArrayLike, StateSpec(vector=True)]	The initial polarisation.	0.0
a	Annotated[float, ParameterSpec()]	The coefficient of |p|^2 / 2 of the energy, negative in the polar phase.	-1.0
c	Annotated[float, ParameterSpec(domain=open(0, oo))]	The coefficient of |p|^4 / 4 of the energy.	1.0
K	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The elastic constant.	1.0
xi	Annotated[float, ParameterSpec()]	The flow-alignment parameter.	0.0
gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The rotational viscosity.	1.0
eta	Annotated[float, ParameterSpec(domain=open(0, oo))]	The viscosity of the flow.	1.0
zeta	Annotated[float, ParameterSpec()]	The activity.	-0.5

Attributes:

Name	Type	Description
default_visual_		The visual of the element when it sets none.

Examples:

>>> from skcomplex.systems.physics.active_matter import PolarGelSiteElement
>>> PolarGelSiteElement('site').xi
0.0

Source code in src/skcomplex/systems/physics/active_matter/_nematics.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    p: Annotated[npt.ArrayLike, StateSpec(vector=True)] = 0.0,
    a: Annotated[float, ParameterSpec()] = -1.0,
    c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    xi: Annotated[float, ParameterSpec()] = 0.0,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    eta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    zeta: Annotated[float, ParameterSpec()] = -0.5,
) -> None:
    super().__init__(label, visual=visual)
    self.p = p
    self.a = a
    self.c = c
    self.K = K
    self.xi = xi
    self.gamma = gamma
    self.eta = eta
    self.zeta = zeta

PropulsionStepInteraction(label, *, target, speed=1.0, state='x', heading='theta')

Bases: MapInteraction

Step x(t + 1) = x + speed (cos theta, sin theta) of a position in the plane along a heading, in discrete time.

The step reads the heading at the start of the step, as every map of the step does.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the particle or the group whose positions move.	required
speed	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The distance moved in one step.	1.0
state	str	The name of the vector state that holds the position.	'x'
heading	str	The name of the state that holds the heading.	'theta'

Examples:

>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var, vector
>>> from skcomplex.systems.physics.common.particles import PropulsionStepInteraction
>>> step = PropulsionStepInteraction('propulsion', target='flock').prepare()
>>> target = SimpleNamespace(x=vector([var('x1'), var('x2')]), theta=var('theta'))
>>> [law.value for law in step.run(SimpleNamespace(target=target, speed=var('v'))).laws_]
[(v*cos(theta) + x1, v*sin(theta) + x2)]

Source code in src/skcomplex/systems/physics/common/particles/_propulsion_step.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    speed: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    state: str = 'x',
    heading: str = 'theta',
) -> None:
    super().__init__(label, target=target, state=state, function=self._step, mode='increment')
    self.speed = speed
    self.heading = heading

SelfPropulsionInteraction(label, *, target, speed=1.0, state='x', heading='theta')

Bases: RateInteraction

Self-propulsion x' += speed (cos theta, sin theta) of a position in the plane along a heading theta.

The heading is a state of the same target, which other laws turn and diffuse.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the particle or the group whose positions move.	required
speed	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The self-propulsion speed.	1.0
state	str	The name of the vector state that holds the position.	'x'
heading	str	The name of the state that holds the heading.	'theta'

Examples:

>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var, vector
>>> from skcomplex.systems.physics.common.particles import SelfPropulsionInteraction
>>> propulsion = SelfPropulsionInteraction('propulsion', target='swimmers').prepare()
>>> target = SimpleNamespace(x=vector([var('x1'), var('x2')]), theta=var('theta'))
>>> [law.value for law in propulsion.run(SimpleNamespace(target=target, speed=var('v'))).laws_]
[(v*cos(theta), v*sin(theta))]

Source code in src/skcomplex/systems/physics/common/particles/_self_propulsion.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    speed: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    state: str = 'x',
    heading: str = 'theta',
) -> None:
    super().__init__(label, target=target, state=state)
    self.speed = speed
    self.heading = heading

SquirmerElement(label, visual=POINTS, *, x=0.0, p=0.0)

Bases: Element

Spherical squirmer, a microswimmer at the position x that swims along the unit orientation p.

The slip of its surface sets how it swims and stirs the fluid, and SquirmerSlipInteraction turns the slip of every squirmer into the motion of all of them.

Parameters:

Name	Type	Description	Default
label	str	The label of the squirmer or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
p	Annotated[ArrayLike, StateSpec(vector=True)]	The initial orientation, a unit vector.	0.0

Examples:

>>> from skcomplex.systems.physics.active_matter import SquirmerElement
>>> SquirmerElement('swimmer', p=[0.0, 0.0, 1.0]).p
[0.0, 0.0, 1.0]

Source code in src/skcomplex/systems/physics/active_matter/_squirmers.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    p: Annotated[npt.ArrayLike, StateSpec(vector=True)] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.p = p

SquirmerSlipInteraction(label, *, target, B1=1.0, B2=0.0, epsilon=1.0, radius=1.0, viscosity=1.0)

Bases: MobilityInteraction

Slip of spherical squirmers, the B1 and B2 modes that make them swim and stir each other through the mobility.

A squirmer of radius a with the tangential surface slip B1 sin(theta) + B2 sin(theta) cos(theta) about its orientation p (Lighthill 1952, Blake 1971) swims at 2 B1 / 3 along p. Far from it, the slip stirs the fluid as the stresslet 4 pi eta a^2 B2 (p p - I / 3) and the potential dipole B1 a^3 p / 3 (Ishikawa, Simmonds, and Pedley 2006). So beta = B2 / B1 < 0 makes a pusher, which pushes fluid out along its axis, and beta > 0 a puller, which draws it in. The flows of the other squirmers move every squirmer and turn it by half their vorticity, through MobilityInteraction. A Weeks-Chandler-Andersen repulsion of strength epsilon keeps the spheres apart, within 2^(1/6) diameters of the current radius.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the group of squirmers.	required
B1	Annotated[float, ParameterSpec()]	The first slip mode, which sets the swimming speed 2 B1 / 3.	1.0
B2	Annotated[float, ParameterSpec()]	The second slip mode, negative for pushers and positive for pullers.	0.0
epsilon	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the repulsion between spheres closer than 2^(1/6) diameters.	1.0
radius	Annotated[float, ParameterSpec(domain=open(0, oo))]	The radius of the squirmers.	1.0
viscosity	Annotated[float, ParameterSpec(domain=open(0, oo))]	The viscosity of the fluid.	1.0

Examples:

>>> from skcomplex.systems.physics.active_matter import SquirmerSlipInteraction
>>> SquirmerSlipInteraction('slip', target='squirmers', B2=-5.0).B2
-5.0

Source code in src/skcomplex/systems/physics/active_matter/_squirmers.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    B1: Annotated[float, ParameterSpec()] = 1.0,
    B2: Annotated[float, ParameterSpec()] = 0.0,
    epsilon: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    radius: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    viscosity: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(
        label,
        target=target,
        orientation='p',
        velocity=self._get_swimming,
        stresslet=self._get_stresslet_of,
        dipole=self._get_dipole_of,
        potential=self._get_repulsion,
        radius=radius,
        viscosity=viscosity,
    )
    self.B1 = B1
    self.B2 = B2
    self.epsilon = epsilon

TumblerElement(label, visual=VisualSpec('agents', {'color_by': 'theta'}), *, theta=0.0, speed=0.1, tumble_rate=0.1, width=10.0)

Bases: Element

Run-and-tumble agent in a periodic square, which runs straight and tumbles to a random direction.

Every step the agent moves speed along its direction, then tumbles with probability 1 - exp(-tumble_rate) to a direction drawn uniformly, as bacteria such as E. coli do.

Parameters:

Name	Type	Description	Default
label	str	The label of the agent or group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	VisualSpec('agents', {'color_by': 'theta'})
theta	Annotated[float, StateSpec(role=direction)]	The initial direction.	0.0
speed	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The distance run per step.	0.1
tumble_rate	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of tumbles per step.	0.1
width	Annotated[float, ParameterSpec(domain=open(0, oo))]	The side of the periodic square.	10.0

Examples:

>>> from skcomplex.systems.physics.active_matter import TumblerElement
>>> TumblerElement('bacterium').speed
0.1

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('agents', {'color_by': 'theta'}),
    *,
    theta: Annotated[float, StateSpec(role='direction')] = 0.0,
    speed: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    tumble_rate: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    width: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
) -> None:
    super().__init__(label, visual=visual)
    self.theta = theta
    self.speed = speed
    self.tumble_rate = tumble_rate
    self.width = width

run_and_tumble(agent, system)

Run one step along the direction, then tumble at random.

Parameters:

Name	Type	Description	Default
agent	TumblerView	The agent.	required
system	SystemView	The system.	required

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

def run_and_tumble(self, agent: TumblerView, system: SystemView) -> None:
    """Run one step along the direction, then tumble at random.

    Args:
        agent:
            The agent.

        system:
            The system.
    """
    position = np.asarray(agent.location, dtype=float)
    heading = np.array([math.cos(agent.theta), math.sin(agent.theta)])
    agent.move_to(np.mod(position + agent.speed * heading, agent.width))
    if system.rng.random() < -math.expm1(-agent.tumble_rate):
        agent.theta = float(system.rng.uniform(0.0, 2 * math.pi))

TumblerView(agent, stepper)

Bases: AgentView

Run-and-tumble agent as rules see it, with its direction and the parameters of its group.

Parameters:

Name	Type	Description	Default
agent	Agent	The Mesa agent.	required
stepper	MesaSolver	The stepper that runs the system.	required

Source code in src/skcomplex/core/engine/solvers/_mesa.py

def __init__(self, agent: mesa.Agent, stepper: 'MesaSolver') -> None:
    self._agent = agent
    self._stepper = stepper

VicsekParticleElement(label, visual=None, *, x=0.0, theta=0.0, eta=0.1)

Bases: Element

Particle of a Vicsek flock, whose heading becomes the mean heading of its neighbours plus noise every step.

HeadingAlignmentInteraction sets the new heading to the circular mean of the headings within the unit distance, itself included, plus a uniform angle in [-eta pi, eta pi]. A propulsion step moves the particle along its heading.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the particle is drawn, or None for default_visual_.	None
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
theta	Annotated[float, StateSpec(period=2 * pi, role=direction)]	The initial heading.	0.0
eta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the noise, as a fraction of pi.	0.1

Attributes:

Name	Type	Description
default_visual_		The visual of the particle when it sets none.

Examples:

>>> from skcomplex.systems.physics.active_matter import VicsekParticleElement
>>> VicsekParticleElement('bird').eta
0.1

Source code in src/skcomplex/systems/physics/active_matter/_flocks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    theta: Annotated[float, StateSpec(period=2 * pi, role='direction')] = 0.0,
    eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.theta = theta
    self.eta = eta

active_brownian_particles(n_particles=200, *, speed=3.0, D_t=0.1, D_r=1.0, size=50.0, random_state=0)

Return active Brownian particles in the plane, which swim at a constant speed in a diffusing direction.

Their mean squared displacement is ballistic below the persistence time 1 / D_r and diffusive above it, with the effective diffusion coefficient D_t + v^2 / (2 D_r).

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	200
speed	float	The self-propulsion speed.	3.0
D_t	float	The translational diffusion coefficient.	0.1
D_r	float	The rotational diffusion coefficient.	1.0
size	float	The side of the square in which the particles start.	50.0
random_state	int | Generator | None	The seed or generator of the initial positions and orientations.	0

Returns:

Name	Type	Description
system	System	The system, with the positions swimmers.x and orientations swimmers.theta.

Examples:

>>> from skcomplex.systems.physics.active_matter import active_brownian_particles
>>> active_brownian_particles(10).elements[0].prepare().size_
10

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'self-propulsion', 'diffusion'),
    references=('Argun et al. 2021', 'Bechinger et al. 2016', 'Volpe, Gigan, and Volpe 2014'),
    defaults={'until': 20.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def active_brownian_particles(
    n_particles: int = 200,
    *,
    speed: float = 3.0,
    D_t: float = 0.1,
    D_r: float = 1.0,
    size: float = 50.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return active Brownian particles in the plane, which swim at a constant speed in a diffusing direction.

    Their mean squared displacement is ballistic below the persistence time `1 / D_r` and diffusive above it, with
    the effective diffusion coefficient `D_t + v^2 / (2 D_r)`.

    Args:
        n_particles:
            The number of particles.

        speed:
            The self-propulsion speed.

        D_t:
            The translational diffusion coefficient.

        D_r:
            The rotational diffusion coefficient.

        size:
            The side of the square in which the particles start.

        random_state:
            The seed or generator of the initial positions and orientations.

    Returns:
        system:
            The system, with the positions `swimmers.x` and orientations `swimmers.theta`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import active_brownian_particles
        >>> active_brownian_particles(10).elements[0].prepare().size_
        10
    """
    positions, angles = _swimmers(n_particles, size, random_state)
    group = GroupElement(
        'swimmers',
        element=ActiveParticleElement,
        size=n_particles,
        values={'x': positions, 'theta': angles, 'D_t': D_t, 'D_r': D_r},
    )
    propulsion = SelfPropulsionInteraction('propulsion', target='swimmers', speed=speed)
    return System(elements=[group], interactions=[propulsion], space=ContinuousSpace(n_dim=2))

active_gels(shape=None, *, order='nematic', spacing=None, a=None, c=None, K=1.0, xi=None, gamma=1.0, eta=1.0, zeta=None, film=False, perturbation=None, random_state=0)

Return an active gel, a nematic tensor or a polarisation transported in the Stokes flow of its own stresses.

On a periodic lattice with spectral operators, the order parameter follows its co-rotational derivative in the flow, Section III of Marchetti et al. The flow solves the overdamped Stokes equation of the elastic stress and the active stress, -zeta Q for the nematic and -zeta p p for the polar gel, and the elastic energy relaxes the order. An extensile nematic, zeta > 0, is unstable to bend, and a contractile polar gel, zeta < 0, is unstable to splay above a threshold.

The defaults of the nematic are the chaotic state of Giomi et al. 2012, cited by Marchetti et al., whose free energy has the same form: a box of side 10, the critical density 3 pi / 2, the density twice that, which give a and c, the activity 3 times the squared density, and the flow alignment 0.1, in units of the elastic constant, the rotational viscosity, and the viscosity. The density stays uniform here. At these defaults pairs of +1/2 and -1/2 defects move in a chaotic flow, active turbulence.

The nematic starts from the tensor along the first axis at the order sqrt(-2 a / c), and the polar gel from the polarisation along the first axis at the length sqrt(-a / c), both with noise. With film=True the polarisation along the second axis is odd about the sites 0 and shape[1] / 2 of the second axis, and the other one even, which holds the polarisation along the first axis at these two walls, with stress-free flow there.

Parameters:

Name	Type	Description	Default
shape	Sequence[int] | None	The number of sites along every axis of the periodic lattice, in two dimensions, or None for (128, 128) for the nematic and (32, 32) for the polar gel.	None
order	str	'nematic' or 'polar'.	'nematic'
spacing	float | None	The distance between neighbouring sites, or None for a nematic box of side 10 along the first axis and 1 for the polar gel.	None
a	float | None	The coefficient of the quadratic energy, negative in the ordered phase, or None for -3 pi / 4 for the nematic and -1 for the polar gel.	None
c	float | None	The coefficient of the quartic energy, or None for 3 pi for the nematic and 1 for the polar gel.	None
K	float	The elastic constant.	1.0
xi	float | None	The flow-alignment parameter, or None for 0.1 for the nematic and 0 for the polar gel.	None
gamma	float	The rotational viscosity.	1.0
eta	float	The viscosity of the flow.	1.0
zeta	float | None	The activity, or None for 3 (3 pi)^2 for the nematic and -0.5 for the polar gel.	None
film	bool	Whether the polar gel lies between walls of fixed polarisation.	False
perturbation	float | None	The standard deviation of the initial noise of the order parameter, or None for 0.1 for the nematic and 0.01 for the polar gel.	None
random_state	int | Generator | None	The seed or generator of the initial noise.	0

Returns:

Name	Type	Description
system	System	The system, with the field gel.Q or gel.p.

Raises:

Type	Description
ValidationError	The order is unknown, or a nematic is asked to lie in a film.

Examples:

>>> from skcomplex.systems.physics.active_matter import active_gels
>>> active_gels((8, 8), order='polar').elements[0].values['zeta']
-0.5
>>> active_gels(order='polar').elements[0].space.shape
(32, 32)

Source code in src/skcomplex/systems/physics/active_matter/_nematics.py

@register_system(
    field='physics/active_matter',
    tags=('field', 'active matter', 'nematic', 'polar', 'tensor', 'gel', 'stokes flow', 'active turbulence', 'pde'),
    references=(
        'Marchetti et al. 2013',
        'arXiv:1207.2929',
        'Kruse et al. 2005',
        'Voituriez, Joanny, and Prost 2005',
        'Giomi et al. 2012',
        'arXiv:1110.4338',
    ),
    defaults={'until': 0.5, 'dt': 0.025, 'solver': RK4Solver(dt=5e-4)},
    settings=_gel_settings,
    observables=lambda system: [
        ObservableSpec('u_x', _velocity_x, group='gel', reduce=_sites),
        ObservableSpec('u_y', _velocity_y, group='gel', reduce=_sites),
    ],
)
def active_gels(
    shape: Sequence[int] | None = None,
    *,
    order: str = 'nematic',
    spacing: float | None = None,
    a: float | None = None,
    c: float | None = None,
    K: float = 1.0,
    xi: float | None = None,
    gamma: float = 1.0,
    eta: float = 1.0,
    zeta: float | None = None,
    film: bool = False,
    perturbation: float | None = None,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return an active gel, a nematic tensor or a polarisation transported in the Stokes flow of its own stresses.

    On a periodic lattice with spectral operators, the order parameter follows its co-rotational derivative in the flow,
    Section III of Marchetti et al. The flow solves the overdamped Stokes equation of the elastic stress and the active
    stress, `-zeta Q` for the nematic and `-zeta p p` for the polar gel, and the elastic energy relaxes the order. An
    extensile nematic, `zeta > 0`, is unstable to bend, and a contractile polar gel, `zeta < 0`, is unstable to splay
    above a threshold.

    The defaults of the nematic are the chaotic state of Giomi et al. 2012, cited by Marchetti et al., whose free energy
    has the same form: a box of side 10, the critical density `3 pi / 2`, the density twice that, which give `a` and
    `c`, the activity `3` times the squared density, and the flow alignment 0.1, in units of the elastic constant, the
    rotational viscosity, and the viscosity. The density stays uniform here. At these defaults pairs of `+1/2` and
    `-1/2` defects move in a chaotic flow, active turbulence.

    The nematic starts from the tensor along the first axis at the order `sqrt(-2 a / c)`, and the polar gel from the
    polarisation along the first axis at the length `sqrt(-a / c)`, both with noise. With `film=True` the polarisation
    along the second axis is odd about the sites 0 and `shape[1] / 2` of the second axis, and the other one even, which
    holds the polarisation along the first axis at these two walls, with stress-free flow there.

    Args:
        shape:
            The number of sites along every axis of the periodic lattice, in two dimensions, or `None` for `(128, 128)`
            for the nematic and `(32, 32)` for the polar gel.

        order:
            `'nematic'` or `'polar'`.

        spacing:
            The distance between neighbouring sites, or `None` for a nematic box of side 10 along the first axis and 1
            for the polar gel.

        a:
            The coefficient of the quadratic energy, negative in the ordered phase, or `None` for `-3 pi / 4` for the
            nematic and -1 for the polar gel.

        c:
            The coefficient of the quartic energy, or `None` for `3 pi` for the nematic and 1 for the polar gel.

        K:
            The elastic constant.

        xi:
            The flow-alignment parameter, or `None` for 0.1 for the nematic and 0 for the polar gel.

        gamma:
            The rotational viscosity.

        eta:
            The viscosity of the flow.

        zeta:
            The activity, or `None` for `3 (3 pi)^2` for the nematic and -0.5 for the polar gel.

        film:
            Whether the polar gel lies between walls of fixed polarisation.

        perturbation:
            The standard deviation of the initial noise of the order parameter, or `None` for 0.1 for the nematic and
            0.01 for the polar gel.

        random_state:
            The seed or generator of the initial noise.

    Returns:
        system:
            The system, with the field `gel.Q` or `gel.p`.

    Raises:
        ValidationError:
            The order is unknown, or a nematic is asked to lie in a film.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import active_gels
        >>> active_gels((8, 8), order='polar').elements[0].values['zeta']
        -0.5
        >>> active_gels(order='polar').elements[0].space.shape
        (32, 32)
    """
    if order not in GEL_ORDERS:
        error_msg = f'The order of the gel is {order!r}. Use one of {", ".join(GEL_ORDERS)}.'
        raise ValidationError(error_msg)
    if film and order != 'polar':
        error_msg = 'A film of the gel holds a polarisation at its walls. Use order=polar.'
        raise ValidationError(error_msg)
    polar = order == 'polar'
    shape = ((32, 32) if polar else (128, 128)) if shape is None else tuple(shape)
    sites = int(np.prod(shape))
    rng = np.random.default_rng(random_state)
    quadratic = (-1.0 if polar else (NEMATIC_CRITICAL_DENSITY - NEMATIC_DENSITY) / 2) if a is None else a
    cubic = (1.0 if polar else NEMATIC_DENSITY) if c is None else c
    activity = (-0.5 if polar else NEMATIC_ACTIVITY * NEMATIC_DENSITY**2) if zeta is None else zeta
    alignment = (0.0 if polar else NEMATIC_XI) if xi is None else xi
    distance = (1.0 if polar else NEMATIC_BOX / shape[0]) if spacing is None else spacing
    amplitude = (0.01 if polar else 0.1) if perturbation is None else perturbation
    noise = amplitude * rng.standard_normal((sites, 2))
    size = np.sqrt(-quadratic / cubic) if polar else np.sqrt(-2 * quadratic / cubic) / 2
    start = np.column_stack([np.full(sites, size), np.zeros(sites)]) + noise if quadratic < 0 else noise
    values = {
        'p' if polar else 'Q': start,
        'a': quadratic,
        'c': cubic,
        'K': K,
        'xi': alignment,
        'gamma': gamma,
        'eta': eta,
        'zeta': activity,
    }
    symmetry = {'p': ((None, 'even'), (None, 'odd'))} if film else None
    gel = FieldElement(
        'gel',
        element=PolarGelSiteElement if polar else NematicGelSiteElement,
        space=LatticeSpace(shape, spacing=distance),
        values=values,
        discretization='spectral',
        symmetry=symmetry,
    )
    return System(elements=[gel])

active_ornstein_uhlenbeck_particles(n_particles=200, *, tau=1.0, D=1.0, D_t=0.0, size=50.0, random_state=0)

Return active Ornstein-Uhlenbeck particles in the plane, driven by exponentially correlated propulsion.

The propulsion starts in its stationary distribution, normal with variance D / tau per component.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	200
tau	float	The persistence time.	1.0
D	float	The active diffusion coefficient.	1.0
D_t	float	The translational diffusion coefficient.	0.0
size	float	The side of the square in which the particles start.	50.0
random_state	int | Generator | None	The seed or generator of the initial positions and propulsions.	0

Returns:

Name	Type	Description
system	System	The system, with the positions swimmers.x and propulsions swimmers.u.

Examples:

>>> from skcomplex.systems.physics.active_matter import active_ornstein_uhlenbeck_particles
>>> active_ornstein_uhlenbeck_particles(10, tau=0.5).elements[0].values['tau']
0.5

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'self-propulsion', 'colored noise'),
    references=('Bechinger et al. 2016', 'Fodor et al. 2016'),
    defaults={'until': 20.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def active_ornstein_uhlenbeck_particles(
    n_particles: int = 200,
    *,
    tau: float = 1.0,
    D: float = 1.0,
    D_t: float = 0.0,
    size: float = 50.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return active Ornstein-Uhlenbeck particles in the plane, driven by exponentially correlated propulsion.

    The propulsion starts in its stationary distribution, normal with variance `D / tau` per component.

    Args:
        n_particles:
            The number of particles.

        tau:
            The persistence time.

        D:
            The active diffusion coefficient.

        D_t:
            The translational diffusion coefficient.

        size:
            The side of the square in which the particles start.

        random_state:
            The seed or generator of the initial positions and propulsions.

    Returns:
        system:
            The system, with the positions `swimmers.x` and propulsions `swimmers.u`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import active_ornstein_uhlenbeck_particles
        >>> active_ornstein_uhlenbeck_particles(10, tau=0.5).elements[0].values['tau']
        0.5
    """
    rng = np.random.default_rng(random_state)
    positions = rng.uniform(0.0, size, (n_particles, 2))
    propulsion = rng.normal(0.0, math.sqrt(D / tau), (n_particles, 2))
    group = GroupElement(
        'swimmers',
        element=OrnsteinUhlenbeckParticleElement,
        size=n_particles,
        values={'x': positions, 'u': propulsion, 'tau': tau, 'D': D, 'D_t': D_t},
    )
    return System(elements=[group], space=ContinuousSpace(n_dim=2))

active_passive_flocking(n_agents=342, *, n_active=45, size=1.0, radius=0.025, epsilon=25.0, exponent=1.5, flocking_radius=0.1, alpha=1.0, v0=0.02, mu=10.0, eta=0.1, time_step=TIME_STEP, initial_speed=1.0, random_state=0)

Return the active and passive agents of Charbonneau (2017, section 10.4) in a periodic box.

All agents repel each other within 2 r_0 with the force epsilon (1 - r / (2 r_0))^exponent, and all feel the force mu (v_t - |v|) v / |v|, with the target speed v_t = v0 for an active agent and 0 for a passive one. An active agent also feels the flocking force alpha V / |V|, with V the sum of the velocities of the agents of both kinds within the flocking radius, itself excluded, and a random force. Agents start at uniform random positions with velocity components uniform up to initial_speed. The defaults are those of figure 10.7 at the compactness C = pi N r_0^2 / size^2 = 0.672. The values exponent=2.5, epsilon=100, mu=0.05, radius=1, flocking_radius=4, and v0=1 are those of Silverberg et al. (2013).

Parameters:

Name	Type	Description	Default
n_agents	int	The number of agents.	342
n_active	int	The number of active agents.	45
size	float	The side of the periodic box.	1.0
radius	float	The radius r_0 of an agent, so agents overlap below the distance 2 r_0.	0.025
epsilon	float	The strength of the repulsion.	25.0
exponent	float	The exponent of the repulsion.	1.5
flocking_radius	float	The radius within which an active agent reads the velocities of the others.	0.1
alpha	float	The size of the flocking force on an active agent, 0 for no flocking.	1.0
v0	float	The target speed of the active agents.	0.02
mu	float	The strength of the propulsion.	10.0
eta	float	The half-width of the uniform random force per step of the book. The noise is a Wiener process with the same variance per step, of amplitude eta sqrt(time_step / 3).	0.1
time_step	float	The step of the solver that eta refers to.	TIME_STEP
initial_speed	float	The largest size of a component of an initial velocity.	1.0
random_state	int | Generator | None	The seed or generator of the initial positions and velocities.	0

Returns:

Name	Type	Description
system	System	The system, with the group agents that holds the active agents first, the positions agents.x, the velocities agents.x_d1, the marks agents.active, and the observables polarization, active_distance, and passive_distance.

Raises:

Type	Description
ValidationError	The number of active agents is not between 1 and n_agents.

Examples:

>>> from skcomplex.systems.physics.active_matter import active_passive_flocking
>>> [element.label for element in active_passive_flocking(114).elements]
['agents']

Source code in src/skcomplex/systems/physics/active_matter/_active_passive.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'flocking', 'phase separation', 'pair interaction'),
    references=('Charbonneau 2017', 'Silverberg et al. 2013'),
    defaults={'until': 50.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=TIME_STEP), 'backend': NumbaBackend()},
    observables=_observables,
)
def active_passive_flocking(
    n_agents: int = 342,
    *,
    n_active: int = 45,
    size: float = 1.0,
    radius: float = 0.025,
    epsilon: float = 25.0,
    exponent: float = 1.5,
    flocking_radius: float = 0.1,
    alpha: float = 1.0,
    v0: float = 0.02,
    mu: float = 10.0,
    eta: float = 0.1,
    time_step: float = TIME_STEP,
    initial_speed: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the active and passive agents of Charbonneau (2017, section 10.4) in a periodic box.

    All agents repel each other within `2 r_0` with the force `epsilon (1 - r / (2 r_0))^exponent`, and all feel the
    force `mu (v_t - |v|) v / |v|`, with the target speed `v_t = v0` for an active agent and 0 for a passive one. An
    active agent also feels the flocking force `alpha V / |V|`, with `V` the sum of the velocities of the agents of both
    kinds within the flocking radius, itself excluded, and a random force. Agents start at uniform random positions
    with velocity components uniform up to `initial_speed`. The defaults are those of figure 10.7 at the compactness
    `C = pi N r_0^2 / size^2 = 0.672`. The values `exponent=2.5`, `epsilon=100`, `mu=0.05`, `radius=1`,
    `flocking_radius=4`, and `v0=1` are those of Silverberg et al. (2013).

    Args:
        n_agents:
            The number of agents.

        n_active:
            The number of active agents.

        size:
            The side of the periodic box.

        radius:
            The radius `r_0` of an agent, so agents overlap below the distance `2 r_0`.

        epsilon:
            The strength of the repulsion.

        exponent:
            The exponent of the repulsion.

        flocking_radius:
            The radius within which an active agent reads the velocities of the others.

        alpha:
            The size of the flocking force on an active agent, 0 for no flocking.

        v0:
            The target speed of the active agents.

        mu:
            The strength of the propulsion.

        eta:
            The half-width of the uniform random force per step of the book. The noise is a Wiener process with the
            same variance per step, of amplitude `eta sqrt(time_step / 3)`.

        time_step:
            The step of the solver that `eta` refers to.

        initial_speed:
            The largest size of a component of an initial velocity.

        random_state:
            The seed or generator of the initial positions and velocities.

    Returns:
        system:
            The system, with the group `agents` that holds the active agents first, the positions `agents.x`, the
            velocities `agents.x_d1`, the marks `agents.active`, and the observables `polarization`, `active_distance`,
            and `passive_distance`.

    Raises:
        ValidationError:
            The number of active agents is not between 1 and `n_agents`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import active_passive_flocking
        >>> [element.label for element in active_passive_flocking(114).elements]
        ['agents']
    """
    if not 1 <= n_active <= n_agents:
        error_msg = f'The system has {n_active} active agents among {n_agents}. Give at least one, and at most all.'
        raise ValidationError(error_msg)
    rng = np.random.default_rng(random_state)
    agents = GroupElement(
        'agents',
        element=FlockingAgentElement,
        size=n_agents,
        values={
            'x': rng.uniform(0.0, size, (n_agents, 2)),
            'x_d1': rng.uniform(-initial_speed, initial_speed, (n_agents, 2)),
            'active': np.where(np.arange(n_agents) < n_active, ACTIVE_AGENT, PASSIVE_AGENT),
        },
        visual=VisualSpec(
            'points',
            {'position': 'x', 'radius': radius, 'color_by': 'active', 'color_range': (0.0, 1.0)},
        ),
    )
    energy = functools.partial(_repulsion_energy, epsilon=epsilon, radius=radius, exponent=exponent)
    return System(
        elements=[agents],
        interactions=[
            PairCouplingInteraction('repulsion', 'agents', 'agents', cutoff=2 * radius, potential=energy, period=size),
            RateInteraction(
                'propulsion',
                target='agents',
                rate=functools.partial(_propulsion, mu=mu, v0=v0),
                noise=functools.partial(_random_force, amplitude=eta * math.sqrt(time_step / 3)),
            ),
            RateInteraction(
                'flocking',
                target='agents',
                rate=functools.partial(_flocking, alpha=alpha, reach=flocking_radius),
            ),
        ],
        space=ContinuousSpace(extent=(size, size), boundary='periodic'),
    )

active_polymer(n_beads=20, *, v0=1.0, k=10.0, D_t=0.1, D_r=1.0, random_state=0)

Return an active polymer, a Rouse chain of active Brownian beads that each swim along their own orientation.

Bead i follows r_i' = v0 n_i + k (r_{i+1} + r_{i-1} - 2 r_i) + sqrt(2 D_t) xi_i, with its orientation n_i = (cos theta_i, sin theta_i) diffusing with D_r, and the end beads bonded to one neighbour. The bonds cancel in the motion of the centre of mass, whose mean squared displacement is (4 D_t t + 2 v0² (D_r t - 1 + exp(-D_r t)) / D_r²) / N, that of an active Brownian particle divided by N. The activity swells the chain: a passive chain has the mean squared end-to-end distance 2 D_t (N - 1) / k in the plane, and self-propulsion of persistence time 1 / D_r adds to the bond fluctuations and stretches it further. The beads start on a straight line with unit spacing and random orientations.

Parameters:

Name	Type	Description	Default
n_beads	int	The number of beads.	20
v0	float	The self-propulsion speed of every bead.	1.0
k	float	The stiffness of the bonds.	10.0
D_t	float	The translational diffusion coefficient.	0.1
D_r	float	The rotational diffusion coefficient.	1.0
random_state	int | Generator | None	The seed or generator of the initial orientations.	0

Returns:

Name	Type	Description
system	System	The system, with the coordinates chain.x and chain.y and the orientations chain.theta of the beads.

Examples:

>>> from skcomplex.systems.physics.active_matter import active_polymer
>>> active_polymer(5).elements[0].prepare().size_
5

Source code in src/skcomplex/systems/physics/active_matter/_polymers.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'polymer', 'self-propulsion', 'rouse chain'),
    references=('Winkler and Gompper 2020', 'Osmanović and Rabin 2017', 'Eisenstecken et al. 2016'),
    defaults={'until': 50.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.005)},
)
def active_polymer(
    n_beads: int = 20,
    *,
    v0: float = 1.0,
    k: float = 10.0,
    D_t: float = 0.1,
    D_r: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return an active polymer, a Rouse chain of active Brownian beads that each swim along their own orientation.

    Bead `i` follows `r_i' = v0 n_i + k (r_{i+1} + r_{i-1} - 2 r_i) + sqrt(2 D_t) xi_i`, with its orientation
    `n_i = (cos theta_i, sin theta_i)` diffusing with `D_r`, and the end beads bonded to one neighbour. The bonds cancel
    in the motion of the centre of mass, whose mean squared displacement is
    `(4 D_t t + 2 v0² (D_r t - 1 + exp(-D_r t)) / D_r²) / N`, that of an active Brownian particle divided by `N`. The
    activity swells the chain: a passive chain has the mean squared end-to-end distance `2 D_t (N - 1) / k` in the
    plane, and self-propulsion of persistence time `1 / D_r` adds to the bond fluctuations and stretches it further. The
    beads start on a straight line with unit spacing and random orientations.

    Args:
        n_beads:
            The number of beads.

        v0:
            The self-propulsion speed of every bead.

        k:
            The stiffness of the bonds.

        D_t:
            The translational diffusion coefficient.

        D_r:
            The rotational diffusion coefficient.

        random_state:
            The seed or generator of the initial orientations.

    Returns:
        system:
            The system, with the coordinates `chain.x` and `chain.y` and the orientations `chain.theta` of the beads.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import active_polymer
        >>> active_polymer(5).elements[0].prepare().size_
        5
    """
    angles = np.random.default_rng(random_state).uniform(0.0, 2 * np.pi, n_beads)
    chain = FieldElement(
        'chain',
        element=ActiveBeadElement,
        space=LatticeSpace(n_beads, boundary='reflecting'),
        boundary={'x': 'zero-flux', 'y': 'zero-flux'},
        values={
            'x': np.arange(n_beads, dtype=float),
            'y': np.zeros(n_beads),
            'theta': angles,
            'v0': v0,
            'k': k,
            'D_t': D_t,
            'D_r': D_r,
        },
    )
    return System(elements=[chain])

active_system_b(shape=(64, 64), *, mean=0.4, kappa=1.0, lam=1.0, zeta=4.0, D=0.05, random_state=0)

Return Active System B+, the conserved field theory of the density of active particles.

The density phi follows phi_t = laplacian(mu) - div(zeta laplacian(phi) grad(phi)) + div(sqrt(2 D) eta), with the chemical potential mu = -phi + phi^3 - kappa laplacian(phi) + lam |grad(phi)|^2 and a conserved white noise eta, Equation 4 of Fodor, Jack, and Cates with unit mobility. The terms lam and zeta are the lowest-order active currents that no free energy gives. With zeta = 0 it is Active System B, where lam shifts the coexisting densities and the phases still coarsen. With lam = 1, zeta = 4, and a majority of the dense phase, the Ostwald ripening of the dilute bubbles reverses, so they stop growing at a finite size, a microphase separation. The noise conserves the total of phi, so its mean stays fixed. The lattice has unit spacing, and an explicit step must stay below about 1/32.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of sites along every axis of the periodic lattice.	(64, 64)
mean	float	The mean of phi, which the dynamics conserve.	0.4
kappa	float	The stiffness of the interfaces.	1.0
lam	float	The active term lam |grad(phi)|^2 of the chemical potential.	1.0
zeta	float	The active current zeta laplacian(phi) grad(phi).	4.0
D	float	The temperature of the conserved noise.	0.05
random_state	int | Generator | None	The seed or generator of the initial noise of amplitude 0.05 around the mean.	0

Returns:

Name	Type	Description
system	System	The system, with the field fluid.phi.

Examples:

>>> from skcomplex.systems.physics.active_matter import active_system_b
>>> active_system_b((8, 8), zeta=0.0).elements[0].values['zeta']
0.0

Source code in src/skcomplex/systems/physics/active_matter/_system_b.py

@register_system(
    field='physics/active_matter',
    tags=('field', 'active matter', 'phase separation', 'microphase separation', 'conserved noise', 'system b', 'spde'),
    references=('Fodor, Jack, and Cates 2022', 'Wittkowski et al. 2014', 'Tjhung, Nardini, and Cates 2018'),
    defaults={'until': 1500.0, 'dt': 30.0, 'solver': EulerMaruyamaSolver(dt=0.02), 'backend': NumbaBackend()},
)
def active_system_b(
    shape: Sequence[int] = (64, 64),
    *,
    mean: float = 0.4,
    kappa: float = 1.0,
    lam: float = 1.0,
    zeta: float = 4.0,
    D: float = 0.05,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return Active System B+, the conserved field theory of the density of active particles.

    The density `phi` follows `phi_t = laplacian(mu) - div(zeta laplacian(phi) grad(phi)) + div(sqrt(2 D) eta)`, with
    the chemical potential `mu = -phi + phi^3 - kappa laplacian(phi) + lam |grad(phi)|^2` and a conserved white noise
    `eta`, Equation 4 of Fodor, Jack, and Cates with unit mobility. The terms `lam` and `zeta` are the lowest-order
    active currents that no free energy gives. With `zeta = 0` it is Active System B, where `lam` shifts the coexisting
    densities and the phases still coarsen. With `lam = 1`, `zeta = 4`, and a majority of the dense phase, the
    Ostwald ripening of the dilute bubbles reverses, so they stop growing at a finite size, a microphase separation.
    The noise conserves the total of `phi`, so its mean stays fixed. The lattice has unit spacing, and an explicit step
    must stay below about 1/32.

    Args:
        shape:
            The number of sites along every axis of the periodic lattice.

        mean:
            The mean of `phi`, which the dynamics conserve.

        kappa:
            The stiffness of the interfaces.

        lam:
            The active term `lam |grad(phi)|^2` of the chemical potential.

        zeta:
            The active current `zeta laplacian(phi) grad(phi)`.

        D:
            The temperature of the conserved noise.

        random_state:
            The seed or generator of the initial noise of amplitude 0.05 around the mean.

    Returns:
        system:
            The system, with the field `fluid.phi`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import active_system_b
        >>> active_system_b((8, 8), zeta=0.0).elements[0].values['zeta']
        0.0
    """
    shape = tuple(shape)
    rng = np.random.default_rng(random_state)
    values = {
        'phi': mean + 0.05 * rng.uniform(-1.0, 1.0, shape),
        'kappa': kappa,
        'lam': lam,
        'zeta': zeta,
        'D': D,
    }
    return System(
        elements=[FieldElement('fluid', element=ActiveScalarElement, space=LatticeSpace(shape), values=values)],
    )

active_viscoelastic_gel(n_sites=128, *, length=10.0, density=1.0, perturbation=0.01, D=0.1, gamma=1.0, eta=1.0, E=0.1, T0=4.0, c_sat=1.0, random_state=0)

Return a one-dimensional active viscoelastic gel, where contractile cells concentrate by dragging the gel.

On a periodic line, the density c of the cells and the displacement u of a Kelvin-Voigt gel follow

c_t + (v c)_x = D c_xx,    u_t = v,    (gamma - eta d_xx) v = d_x (E u_x + Ta(c)),

with the active stress Ta(c) = T0 c / (c + c_sat). The force balance holds at every step, as the modified Helmholtz inversion of the stress gradient with the screening sqrt(gamma / eta). Without stiffness the uniform state loses stability at long wavelengths above the contractile threshold c Ta'(c) > gamma D, and the density concentrates into peaks. The source leaves the active stress open, and the saturating one keeps the peaks finite. The catalog entry streams the velocity as the observable v.

Parameters:

Name	Type	Description	Default
n_sites	int	The number of sites.	128
length	float	The length of the line.	10.0
density	float	The mean density of the cells.	1.0
perturbation	float	The standard deviation of the initial density around its mean.	0.01
D	float	The diffusion coefficient of the cells.	0.1
gamma	float	The friction of the gel with the substrate.	1.0
eta	float	The viscosity of the gel.	1.0
E	float	The stiffness of the gel.	0.1
T0	float	The largest active stress of the cells.	4.0
c_sat	float	The density at which the active stress is half its largest value.	1.0
random_state	int | Generator | None	The seed or generator of the initial density.	0

Returns:

Name	Type	Description
system	System	The system, with the density gel.c and the displacement gel.u.

Examples:

>>> from skcomplex.systems.physics.active_matter import active_viscoelastic_gel
>>> active_viscoelastic_gel(16).elements[0].values['eta']
1.0

Source code in src/skcomplex/systems/physics/active_matter/_gels.py

@register_system(
    field='physics/active_matter',
    tags=('active gel', 'contractile', 'viscoelastic', 'field', 'inversion', 'pattern formation', 'pde'),
    references=('Brauns and Marchetti 2024', 'arXiv:2306.08868'),
    defaults={'until': 30.0, 'dt': 0.1, 'solver': RK4Solver(dt=0.005)},
    observables=lambda system: [ObservableSpec('v', _gel_velocity, group='gel', reduce=_sites)],
)
def active_viscoelastic_gel(
    n_sites: int = 128,
    *,
    length: float = 10.0,
    density: float = 1.0,
    perturbation: float = 0.01,
    D: float = 0.1,
    gamma: float = 1.0,
    eta: float = 1.0,
    E: float = 0.1,
    T0: float = 4.0,
    c_sat: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return a one-dimensional active viscoelastic gel, where contractile cells concentrate by dragging the gel.

    On a periodic line, the density `c` of the cells and the displacement `u` of a Kelvin-Voigt gel follow

        c_t + (v c)_x = D c_xx,    u_t = v,    (gamma - eta d_xx) v = d_x (E u_x + Ta(c)),

    with the active stress `Ta(c) = T0 c / (c + c_sat)`. The force balance holds at every step, as the modified
    Helmholtz inversion of the stress gradient with the screening `sqrt(gamma / eta)`. Without stiffness the uniform
    state loses stability at long wavelengths above the contractile threshold `c Ta'(c) > gamma D`, and the density
    concentrates into peaks. The source leaves the active stress open, and the saturating one keeps the peaks finite.
    The catalog entry streams the velocity as the observable `v`.

    Args:
        n_sites:
            The number of sites.

        length:
            The length of the line.

        density:
            The mean density of the cells.

        perturbation:
            The standard deviation of the initial density around its mean.

        D:
            The diffusion coefficient of the cells.

        gamma:
            The friction of the gel with the substrate.

        eta:
            The viscosity of the gel.

        E:
            The stiffness of the gel.

        T0:
            The largest active stress of the cells.

        c_sat:
            The density at which the active stress is half its largest value.

        random_state:
            The seed or generator of the initial density.

    Returns:
        system:
            The system, with the density `gel.c` and the displacement `gel.u`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import active_viscoelastic_gel
        >>> active_viscoelastic_gel(16).elements[0].values['eta']
        1.0
    """
    rng = np.random.default_rng(random_state)
    values = {
        'c': density + perturbation * rng.standard_normal(n_sites),
        'u': 0.0,
        'D': D,
        'gamma': gamma,
        'eta': eta,
        'E': E,
        'T0': T0,
        'c_sat': c_sat,
    }
    gel = FieldElement(
        'gel',
        element=ContractileGelSiteElement,
        space=LatticeSpace(n_sites, spacing=length / n_sites),
        values=values,
    )
    transport = AdvectionInteraction('transport', target='gel', state='c', drift=_transport, limiter=None, riemann=None)
    diffusion = DiffusionInteraction('diffusion', target='gel', state='c', coefficient=_diffusion)
    displacement = RateInteraction('displacement', target='gel', state='u', rate=_transport)
    return System(elements=[gel], interactions=[transport, diffusion, displacement])

aligning_active_brownian_particles(n_particles=200, *, size=20.0, speed=1.0, gamma=1.0, D_r=0.1, d0=1.0, D_t=0.0, random_state=0)

Return active Brownian particles in a periodic box whose headings align with the neighbours within d0.

The heading of particle i follows dtheta_i = gamma sum_{|x_j - x_i| < d0} sin(theta_j - theta_i) dt + sqrt(2 D_r) dW, and the particle swims at the constant speed along it. In the mean-field approximation of the Fokker-Planck equation, the disordered state is unstable above the density rho_c = 2 D_r / (gamma pi d0^2), so at a fixed density rho the particles move together when gamma exceeds 2 D_r / (pi d0^2 rho). The polarization |<exp(i theta)>| is near 0 in the disordered state and near 1 in the ordered one.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	200
size	float	The side of the periodic box.	20.0
speed	float	The self-propulsion speed.	1.0
gamma	float	The strength of the alignment.	1.0
D_r	float	The rotational diffusion coefficient.	0.1
d0	float	The interaction range.	1.0
D_t	float	The translational diffusion coefficient.	0.0
random_state	int | Generator | None	The seed or generator of the uniform initial positions and headings.	0

Returns:

Name	Type	Description
system	System	The system, with the positions swimmers.x, the headings swimmers.theta, and the observables mean_heading and polarization.

Examples:

>>> from skcomplex.systems.physics.active_matter import aligning_active_brownian_particles
>>> aligning_active_brownian_particles(10, gamma=0.5).interactions[0].gamma
0.5

Source code in src/skcomplex/systems/physics/active_matter/_alignment.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'flocking', 'alignment', 'pair interaction'),
    references=('Bertin 2021', 'Peruani, Deutsch, and Bär 2008', 'Chaté 2020'),
    defaults={'until': 100.0, 'dt': 1.0, 'solver': EulerMaruyamaSolver(dt=0.05)},
    observables=lambda system: [
        ObservableSpec('mean_heading', _heading, group='swimmers'),
        ObservableSpec('polarization', _polarization),
    ],
)
def aligning_active_brownian_particles(
    n_particles: int = 200,
    *,
    size: float = 20.0,
    speed: float = 1.0,
    gamma: float = 1.0,
    D_r: float = 0.1,
    d0: float = 1.0,
    D_t: float = 0.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return active Brownian particles in a periodic box whose headings align with the neighbours within `d0`.

    The heading of particle `i` follows `dtheta_i = gamma sum_{|x_j - x_i| < d0} sin(theta_j - theta_i) dt +
    sqrt(2 D_r) dW`, and the particle swims at the constant `speed` along it. In the mean-field approximation of the
    Fokker-Planck equation, the disordered state is unstable above the density `rho_c = 2 D_r / (gamma pi d0^2)`, so
    at a fixed density `rho` the particles move together when `gamma` exceeds `2 D_r / (pi d0^2 rho)`. The
    polarization `|<exp(i theta)>|` is near 0 in the disordered state and near 1 in the ordered one.

    Args:
        n_particles:
            The number of particles.

        size:
            The side of the periodic box.

        speed:
            The self-propulsion speed.

        gamma:
            The strength of the alignment.

        D_r:
            The rotational diffusion coefficient.

        d0:
            The interaction range.

        D_t:
            The translational diffusion coefficient.

        random_state:
            The seed or generator of the uniform initial positions and headings.

    Returns:
        system:
            The system, with the positions `swimmers.x`, the headings `swimmers.theta`, and the observables
            `mean_heading` and `polarization`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import aligning_active_brownian_particles
        >>> aligning_active_brownian_particles(10, gamma=0.5).interactions[0].gamma
        0.5
    """
    rng = np.random.default_rng(random_state)
    swimmers = GroupElement(
        'swimmers',
        element=ActiveParticleElement,
        size=n_particles,
        values={
            'x': rng.uniform(0.0, size, (n_particles, 2)),
            'theta': rng.uniform(0.0, 2 * math.pi, n_particles),
            'D_t': D_t,
            'D_r': D_r,
        },
        visual=HEADINGS,
    )
    alignment = PolarAlignmentInteraction(
        'align',
        first='swimmers',
        second='swimmers',
        cutoff=d0,
        period=size,
        gamma=gamma,
    )
    return System(
        elements=[swimmers],
        interactions=[alignment, SelfPropulsionInteraction('propulsion', target='swimmers', speed=speed)],
        space=ContinuousSpace(extent=(size, size), boundary='periodic'),
    )

bertin_droz_gregoire(positions=4, angles=16, *, density=1.0, speed=1.0, interaction_range=0.5, tumbling=1.0, noise=0.5, perturbation=0.01)

Return the kinetic theory of flocking, whose polar order grows above a transition density.

Particles stream along their orientation, tumble at a rate by a wrapped normal noise, and align in binary collisions at their mean angle with the same noise. A homogeneous disordered state is unstable above the density pi tumbling (1 - exp(-noise^2 / 2)) / (8 interaction_range speed (exp(-noise^2 / 2) - 2 / 3)).

Parameters:

Name	Type	Description	Default
positions	int	The number of positions along each of the two periodic position axes of unit spacing.	4
angles	int	The number of cells of the orientation axis.	16
density	float	The density of the homogeneous state.	1.0
speed	float	The speed of the particles.	1.0
interaction_range	float	The interaction range.	0.5
tumbling	float	The rate of the tumbling events.	1.0
noise	float	The standard deviation of the noise of tumbling and alignment.	0.5
perturbation	float	The relative amplitude of the initial polar perturbation.	0.01

Returns:

Name	Type	Description
system	System	The system, with the density flock.f.

Examples:

>>> from skcomplex.systems.physics.active_matter import bertin_droz_gregoire
>>> bertin_droz_gregoire(2, 8).elements[0].prepare().size_
32

Source code in src/skcomplex/systems/physics/active_matter/_kinetic_flocks.py

@register_system(
    field='physics/active_matter',
    tags=('field', 'kinetic', 'phase space', 'flocking', 'pde'),
    references=('Bertin, Droz, and Grégoire 2006', 'Bertin, Droz, and Grégoire 2009'),
    defaults={'until': 0.5, 'dt': 0.1, 'solver': RK4Solver(dt=0.05)},
)
def bertin_droz_gregoire(
    positions: int = 4,
    angles: int = 16,
    *,
    density: float = 1.0,
    speed: float = 1.0,
    interaction_range: float = 0.5,
    tumbling: float = 1.0,
    noise: float = 0.5,
    perturbation: float = 0.01,
) -> System:
    """Return the kinetic theory of flocking, whose polar order grows above a transition density.

    Particles stream along their orientation, tumble at a rate by a wrapped normal noise, and align in binary collisions
    at their mean angle with the same noise. A homogeneous disordered state is unstable above the density
    `pi tumbling (1 - exp(-noise^2 / 2)) / (8 interaction_range speed (exp(-noise^2 / 2) - 2 / 3))`.

    Args:
        positions:
            The number of positions along each of the two periodic position axes of unit spacing.

        angles:
            The number of cells of the orientation axis.

        density:
            The density of the homogeneous state.

        speed:
            The speed of the particles.

        interaction_range:
            The interaction range.

        tumbling:
            The rate of the tumbling events.

        noise:
            The standard deviation of the noise of tumbling and alignment.

        perturbation:
            The relative amplitude of the initial polar perturbation.

    Returns:
        system:
            The system, with the density `flock.f`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import bertin_droz_gregoire
        >>> bertin_droz_gregoire(2, 8).elements[0].prepare().size_
        32
    """
    axis = AxisSpec('theta', angles, 0.0, 2 * math.pi, periodic=True)
    space = PhaseSpace(LatticeSpace((positions, positions)), [axis]).prepare()
    theta = space.coordinates_['theta']
    profile = density * (1.0 + perturbation * np.cos(theta)) / (2 * math.pi)
    values = np.broadcast_to(profile, space.shape_).copy()
    flock = FieldElement(
        'flock',
        element=DensityElement,
        space=space,
        values={'f': values},
        discretization='finite-volume',
    )

    def drift(_: SimpleNamespace) -> list[Expression]:
        return [speed * cos(coordinate('theta')), speed * sin(coordinate('theta')), 0]

    return System(
        elements=[flock],
        interactions=[
            AdvectionInteraction('stream', target='flock', state='f', drift=drift),
            AngularDiffusionInteraction('tumble', candidate='flock', angles=angles, tumbling=tumbling, noise=noise),
            AngularAlignmentInteraction(
                'align',
                candidate='flock',
                field='flock',
                angles=angles,
                frequency=interaction_range * speed,
                noise=noise,
            ),
        ],
    )

binary_collision_flocking(n_particles=400, *, size=20.0, speed=1.0, lam=1.0, sigma=0.3, d0=0.5, substeps=10, random_state=0)

Return self-propelled particles that flock through aligning binary collisions, the Bertin-Droz-Grégoire system.

Particles fly straight at a constant speed. Their headings jump by Gaussian angles in self-diffusion events at the rate lam, and two particles that come closer than d0 both take their mean heading plus Gaussian noise. The hydrodynamic equations derived from the Boltzmann equation predict collective motion above the density rho_t = pi lam (1 - exp(-sigma^2 / 2)) / (4 (exp(-sigma^2 / 2) - 2/3)), for speed = 1 and d0 = 1/2, so collective motion needs sigma < sqrt(2 ln(3/2)), about 0.90, at any density.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	400
size	float	The side of the periodic square.	20.0
speed	float	The self-propulsion speed.	1.0
lam	float	The rate of self-diffusion events.	1.0
sigma	float	The standard deviation of the angular noise.	0.3
d0	float	The collision distance.	0.5
substeps	int	The number of substeps per unit of time.	10
random_state	int | Generator | None	The seed or generator of the uniform initial positions and headings.	0

Returns:

Name	Type	Description
system	System	The system, with the agent group birds and the observables mean_cos, mean_sin, and polarization.

Examples:

>>> from skcomplex.systems.physics.active_matter import binary_collision_flocking
>>> binary_collision_flocking(10).interactions[0].sigma
0.3

Source code in src/skcomplex/systems/physics/active_matter/_collisions.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'agents', 'flocking', 'collisions', 'boltzmann equation'),
    references=('Bertin, Droz, and Grégoire 2006', 'Bertin, Droz, and Grégoire 2009', 'Bertin 2021'),
    defaults={'until': 200, 'dt': 1},
    observables=_flocking_observables,
)
def binary_collision_flocking(
    n_particles: int = 400,
    *,
    size: float = 20.0,
    speed: float = 1.0,
    lam: float = 1.0,
    sigma: float = 0.3,
    d0: float = 0.5,
    substeps: int = 10,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return self-propelled particles that flock through aligning binary collisions, the Bertin-Droz-Grégoire system.

    Particles fly straight at a constant speed. Their headings jump by Gaussian angles in self-diffusion events at the
    rate `lam`, and two particles that come closer than `d0` both take their mean heading plus Gaussian noise. The
    hydrodynamic equations derived from the Boltzmann equation predict collective motion above the density
    `rho_t = pi lam (1 - exp(-sigma^2 / 2)) / (4 (exp(-sigma^2 / 2) - 2/3))`, for `speed = 1` and `d0 = 1/2`, so
    collective motion needs `sigma < sqrt(2 ln(3/2))`, about 0.90, at any density.

    Args:
        n_particles:
            The number of particles.

        size:
            The side of the periodic square.

        speed:
            The self-propulsion speed.

        lam:
            The rate of self-diffusion events.

        sigma:
            The standard deviation of the angular noise.

        d0:
            The collision distance.

        substeps:
            The number of substeps per unit of time.

        random_state:
            The seed or generator of the uniform initial positions and headings.

    Returns:
        system:
            The system, with the agent group `birds` and the observables `mean_cos`, `mean_sin`, and `polarization`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import binary_collision_flocking
        >>> binary_collision_flocking(10).interactions[0].sigma
        0.3
    """
    rng = np.random.default_rng(random_state)
    birds = GroupElement(
        'birds',
        element=FlockerElement,
        size=n_particles,
        values={'theta': rng.uniform(0.0, 2 * math.pi, n_particles)},
        location=rng.uniform(0.0, size, (n_particles, 2)),
    )
    collisions = BinaryCollisionsInteraction(
        'collisions',
        target='birds',
        speed=speed,
        lam=lam,
        sigma=sigma,
        d0=d0,
        substeps=float(substeps),
    )
    return System(
        elements=[birds],
        interactions=[collisions],
        space=ContinuousSpace(extent=(size, size), boundary='periodic'),
    )

chiral_active_particles(n_particles=200, *, speed=3.0, omega=2.0, D_t=0.1, D_r=0.1, size=50.0, random_state=0)

Return chiral active particles, which swim along circles of radius v / omega that drift by diffusion.

Their effective diffusion coefficient at long times is D_t + v^2 D_r / (2 (D_r^2 + omega^2)).

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	200
speed	float	The self-propulsion speed.	3.0
omega	float	The angular velocity.	2.0
D_t	float	The translational diffusion coefficient.	0.1
D_r	float	The rotational diffusion coefficient.	0.1
size	float	The side of the square in which the particles start.	50.0
random_state	int | Generator | None	The seed or generator of the initial positions and orientations.	0

Returns:

Name	Type	Description
system	System	The system, with the positions swimmers.x and orientations swimmers.theta.

Examples:

>>> from skcomplex.systems.physics.active_matter import chiral_active_particles
>>> chiral_active_particles(10, omega=1.0).elements[0].values['omega']
1.0

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'self-propulsion', 'chirality'),
    references=('Bechinger et al. 2016', 'van Teeffelen and Löwen 2008'),
    defaults={'until': 20.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def chiral_active_particles(
    n_particles: int = 200,
    *,
    speed: float = 3.0,
    omega: float = 2.0,
    D_t: float = 0.1,
    D_r: float = 0.1,
    size: float = 50.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return chiral active particles, which swim along circles of radius `v / omega` that drift by diffusion.

    Their effective diffusion coefficient at long times is `D_t + v^2 D_r / (2 (D_r^2 + omega^2))`.

    Args:
        n_particles:
            The number of particles.

        speed:
            The self-propulsion speed.

        omega:
            The angular velocity.

        D_t:
            The translational diffusion coefficient.

        D_r:
            The rotational diffusion coefficient.

        size:
            The side of the square in which the particles start.

        random_state:
            The seed or generator of the initial positions and orientations.

    Returns:
        system:
            The system, with the positions `swimmers.x` and orientations `swimmers.theta`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import chiral_active_particles
        >>> chiral_active_particles(10, omega=1.0).elements[0].values['omega']
        1.0
    """
    positions, angles = _swimmers(n_particles, size, random_state)
    group = GroupElement(
        'swimmers',
        element=ActiveParticleElement,
        size=n_particles,
        values={'x': positions, 'theta': angles, 'omega': omega, 'D_t': D_t, 'D_r': D_r},
    )
    propulsion = SelfPropulsionInteraction('propulsion', target='swimmers', speed=speed)
    return System(elements=[group], interactions=[propulsion], space=ContinuousSpace(n_dim=2))

cucker_smale_flocking(n_agents=50, *, K=1.0, beta=0.25, size=10.0, random_state=0)

Return Cucker-Smale flocking, where every agent aligns its velocity with all others, weighted by distance.

The mean velocity is conserved, and for beta <= 1/2 the velocities converge to it from any start.

Parameters:

Name	Type	Description	Default
n_agents	int	The number of agents.	50
K	float	The coupling strength.	1.0
beta	float	The decay exponent of the communication weight.	0.25
size	float	The side of the square in which the agents start.	10.0
random_state	int | Generator | None	The seed or generator of the initial positions and velocities.	0

Returns:

Name	Type	Description
system	System	The system, with the positions birds.x and velocities birds.x_d1.

Examples:

>>> from skcomplex.systems.physics.active_matter import cucker_smale_flocking
>>> cucker_smale_flocking(5).interactions[0].n_agents
5

Source code in src/skcomplex/systems/physics/active_matter/_collective.py

@register_system(
    field='physics/active_matter',
    tags=('flocking', 'ode', 'alignment', 'pair interaction'),
    references=('Cucker and Smale 2007', 'Vicsek and Zafeiris 2012', 'Bellomo et al. 2017'),
    defaults={'until': 20.0, 'dt': 0.5, 'solver': RK4Solver(dt=0.05)},
)
def cucker_smale_flocking(
    n_agents: int = 50,
    *,
    K: float = 1.0,
    beta: float = 0.25,
    size: float = 10.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return Cucker-Smale flocking, where every agent aligns its velocity with all others, weighted by distance.

    The mean velocity is conserved, and for `beta <= 1/2` the velocities converge to it from any start.

    Args:
        n_agents:
            The number of agents.

        K:
            The coupling strength.

        beta:
            The decay exponent of the communication weight.

        size:
            The side of the square in which the agents start.

        random_state:
            The seed or generator of the initial positions and velocities.

    Returns:
        system:
            The system, with the positions `birds.x` and velocities `birds.x_d1`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import cucker_smale_flocking
        >>> cucker_smale_flocking(5).interactions[0].n_agents
        5
    """
    rng = np.random.default_rng(random_state)
    birds = GroupElement(
        'birds',
        element=PointElement,
        size=n_agents,
        values={'x': rng.uniform(0.0, size, (n_agents, 2)), 'x_d1': rng.normal(0.0, 1.0, (n_agents, 2))},
        visual=VisualSpec('points', {'position': 'x', 'color_by': 'x_d1', 'states': ()}),
    )
    alignment = CuckerSmaleAlignmentInteraction(
        'align',
        first='birds',
        second='birds',
        cutoff=1e6,
        K=K,
        beta=beta,
        n_agents=n_agents,
    )
    return System(elements=[birds], interactions=[alignment], space=ContinuousSpace(n_dim=2))

deformable_active_particles(n_particles=20, *, gamma=1.0, a=1.0, b=1.0, kappa=0.5, noise=0.0, random_state=0)

Return deformable self-propelled particles of Ohta and Ohkuma, which move straight or along circles.

Every particle follows x' = v, v' = gamma v - |v|² v - a S v, and S' = -kappa S + b (v v - |v|² I / 2), with the traceless symmetric deformation tensor S. A fast relaxing shape, large kappa, gives straight motion at the speed sqrt(gamma / (1 + a b / (2 kappa))). For a b (gamma - kappa / 2) > kappa² a circular motion exists, at the speed sqrt(gamma - kappa / 2) with the angular velocity sqrt(a b (gamma - kappa / 2) - kappa²) / 2, its elongation turning with its velocity. For the defaults it exists below kappa ≈ 0.78, and particles started from small velocities reach it for kappa <= 0.5, while closer to the threshold straight motion stays stable too. The particles do not interact, and they start at random positions, small random velocities, and no deformation.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	20
gamma	float	The distance from the onset of self-propulsion.	1.0
a	float	The coupling of the velocity to the deformation.	1.0
b	float	The coupling of the deformation to the velocity.	1.0
kappa	float	The relaxation rate of the deformation.	0.5
noise	float	The amplitude of the noise of the velocity.	0.0
random_state	int | Generator | None	The seed or generator of the initial positions and velocities.	0

Returns:

Name	Type	Description
system	System	The system, with the positions cells.x, the velocities cells.v, and the deformations cells.S.

Examples:

>>> from skcomplex.systems.physics.active_matter import deformable_active_particles
>>> deformable_active_particles(5).elements[0].prepare().size_
5

Source code in src/skcomplex/systems/physics/active_matter/_deformable.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'deformable particles', 'tensor', 'circular motion', 'bifurcation'),
    references=('Ohta and Ohkuma 2009', 'Ohta 2017', 'Tarama 2017'),
    defaults={'until': 100.0, 'dt': 0.25, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def deformable_active_particles(
    n_particles: int = 20,
    *,
    gamma: float = 1.0,
    a: float = 1.0,
    b: float = 1.0,
    kappa: float = 0.5,
    noise: float = 0.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return deformable self-propelled particles of Ohta and Ohkuma, which move straight or along circles.

    Every particle follows `x' = v`, `v' = gamma v - |v|² v - a S v`, and `S' = -kappa S + b (v v - |v|² I / 2)`, with
    the traceless symmetric deformation tensor `S`. A fast relaxing shape, large `kappa`, gives straight motion at the
    speed `sqrt(gamma / (1 + a b / (2 kappa)))`. For `a b (gamma - kappa / 2) > kappa²` a circular motion exists, at the
    speed `sqrt(gamma - kappa / 2)` with the angular velocity `sqrt(a b (gamma - kappa / 2) - kappa²) / 2`, its
    elongation turning with its velocity. For the defaults it exists below `kappa ≈ 0.78`, and particles started from
    small velocities reach it for `kappa <= 0.5`, while closer to the threshold straight motion stays stable too. The
    particles do not interact, and they start at random positions, small random velocities, and no deformation.

    Args:
        n_particles:
            The number of particles.

        gamma:
            The distance from the onset of self-propulsion.

        a:
            The coupling of the velocity to the deformation.

        b:
            The coupling of the deformation to the velocity.

        kappa:
            The relaxation rate of the deformation.

        noise:
            The amplitude of the noise of the velocity.

        random_state:
            The seed or generator of the initial positions and velocities.

    Returns:
        system:
            The system, with the positions `cells.x`, the velocities `cells.v`, and the deformations `cells.S`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import deformable_active_particles
        >>> deformable_active_particles(5).elements[0].prepare().size_
        5
    """
    rng = np.random.default_rng(random_state)
    values = {
        'x': rng.uniform(0.0, 20.0, (n_particles, 2)),
        'v': rng.normal(0.0, 0.1, (n_particles, 2)),
        'S': np.zeros((n_particles, 2)),
        'gamma': gamma,
        'a': a,
        'b': b,
        'kappa': kappa,
        'noise': noise,
    }
    cells = GroupElement('cells', element=DeformableParticleElement, size=n_particles, values=values)
    return System(elements=[cells], space=ContinuousSpace(n_dim=2))

depot_brownian_particles(n_particles=50, *, gamma=20.0, q=10.0, c=1.0, d=10.0, D=0.01, random_state=0)

Return Brownian particles with energy depots, which start at rest and take up the speed sqrt(q/gamma - c/d).

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	50
gamma	float	The friction coefficient.	20.0
q	float	The rate of energy intake.	10.0
c	float	The rate of internal energy loss.	1.0
d	float	The rate of conversion of depot energy into motion.	10.0
D	float	The intensity of the velocity noise.	0.01
random_state	int | Generator | None	The seed or generator of the small initial velocities.	0

Returns:

Name	Type	Description
system	System	The system, with the positions walkers.x, velocities walkers.v, and depots walkers.e.

Examples:

>>> from skcomplex.systems.physics.active_matter import depot_brownian_particles
>>> depot_brownian_particles(5).elements[0].values['q']
10.0

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'energy depot', 'self-propulsion'),
    references=('Schweitzer, Ebeling, and Tilch 1998', 'Helbing 2001'),
    defaults={'until': 10.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.001)},
)
def depot_brownian_particles(
    n_particles: int = 50,
    *,
    gamma: float = 20.0,
    q: float = 10.0,
    c: float = 1.0,
    d: float = 10.0,
    D: float = 0.01,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return Brownian particles with energy depots, which start at rest and take up the speed `sqrt(q/gamma - c/d)`.

    Args:
        n_particles:
            The number of particles.

        gamma:
            The friction coefficient.

        q:
            The rate of energy intake.

        c:
            The rate of internal energy loss.

        d:
            The rate of conversion of depot energy into motion.

        D:
            The intensity of the velocity noise.

        random_state:
            The seed or generator of the small initial velocities.

    Returns:
        system:
            The system, with the positions `walkers.x`, velocities `walkers.v`, and depots `walkers.e`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import depot_brownian_particles
        >>> depot_brownian_particles(5).elements[0].values['q']
        10.0
    """
    velocities = np.random.default_rng(random_state).normal(0.0, 0.01, (n_particles, 2))
    group = GroupElement(
        'walkers',
        element=DepotParticleElement,
        size=n_particles,
        values={'v': velocities, 'gamma': gamma, 'q': q, 'c': c, 'd': d, 'D': D},
    )
    return System(elements=[group], space=ContinuousSpace(n_dim=2))

dry_active_nematics(shape=(64, 64), *, density=1.5, critical_density=1.0, a=1.0, c=2.0, K=1.0, gamma=1.0, D=0.5, current=0.5, beta=0.5, perturbation=0.01, random_state=0)

Return a dry active nematic, a conserved density coupled to a symmetric traceless alignment tensor.

On a periodic lattice the density and the tensor follow Section II.2 of Marchetti et al.,

rho_t = -div(J),    J = -D grad(rho) - current div(Q),
Q_t = H / gamma + beta (partial_i partial_j rho - delta_ij laplacian(rho) / 2),

with the molecular field H of the Landau-de Gennes energy, whose coefficient a (1 - rho / rho_c) turns negative above the critical density rho_c. The current along div(Q) is the active current, which no free energy gives.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of sites along every axis of the periodic lattice.	(64, 64)
density	float	The mean density, which the dynamics conserve.	1.5
critical_density	float	The density at which the isotropic state loses stability.	1.0
a	float	The scale of the coefficient of tr(Q^2) / 2 of the energy.	1.0
c	float	The coefficient of tr(Q^2)^2 / 4 of the energy.	2.0
K	float	The elastic constant.	1.0
gamma	float	The rotational viscosity.	1.0
D	float	The diffusion coefficient of the density.	0.5
current	float	The coefficient of the active current along div(Q).	0.5
beta	float	The coupling of the tensor to the curvature of the density.	0.5
perturbation	float	The standard deviation of the initial noise of the density and the tensor.	0.01
random_state	int | Generator | None	The seed or generator of the initial noise.	0

Returns:

Name	Type	Description
system	System	The system, with the fields nematic.rho and nematic.Q.

Examples:

>>> from skcomplex.systems.physics.active_matter import dry_active_nematics
>>> dry_active_nematics((8, 8)).elements[0].values['beta']
0.5

Source code in src/skcomplex/systems/physics/active_matter/_nematics.py

@register_system(
    field='physics/active_matter',
    tags=('field', 'active matter', 'nematic', 'tensor', 'liquid crystal', 'giant number fluctuations', 'pde'),
    references=('Marchetti et al. 2013', 'arXiv:1207.2929'),
    defaults={'until': 100.0, 'dt': 5.0, 'solver': RK4Solver(dt=0.02), 'backend': NumbaBackend()},
)
def dry_active_nematics(
    shape: Sequence[int] = (64, 64),
    *,
    density: float = 1.5,
    critical_density: float = 1.0,
    a: float = 1.0,
    c: float = 2.0,
    K: float = 1.0,
    gamma: float = 1.0,
    D: float = 0.5,
    current: float = 0.5,
    beta: float = 0.5,
    perturbation: float = 0.01,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return a dry active nematic, a conserved density coupled to a symmetric traceless alignment tensor.

    On a periodic lattice the density and the tensor follow Section II.2 of Marchetti et al.,

        rho_t = -div(J),    J = -D grad(rho) - current div(Q),
        Q_t = H / gamma + beta (partial_i partial_j rho - delta_ij laplacian(rho) / 2),

    with the molecular field `H` of the Landau-de Gennes energy, whose coefficient `a (1 - rho / rho_c)` turns negative
    above the critical density `rho_c`. The current along `div(Q)` is the active current, which no free energy gives.

    Args:
        shape:
            The number of sites along every axis of the periodic lattice.

        density:
            The mean density, which the dynamics conserve.

        critical_density:
            The density at which the isotropic state loses stability.

        a:
            The scale of the coefficient of `tr(Q^2) / 2` of the energy.

        c:
            The coefficient of `tr(Q^2)^2 / 4` of the energy.

        K:
            The elastic constant.

        gamma:
            The rotational viscosity.

        D:
            The diffusion coefficient of the density.

        current:
            The coefficient of the active current along `div(Q)`.

        beta:
            The coupling of the tensor to the curvature of the density.

        perturbation:
            The standard deviation of the initial noise of the density and the tensor.

        random_state:
            The seed or generator of the initial noise.

    Returns:
        system:
            The system, with the fields `nematic.rho` and `nematic.Q`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import dry_active_nematics
        >>> dry_active_nematics((8, 8)).elements[0].values['beta']
        0.5
    """
    shape = tuple(shape)
    rng = np.random.default_rng(random_state)
    sites = int(np.prod(shape))
    values = {
        'rho': density + perturbation * rng.standard_normal(sites),
        'Q': perturbation * rng.standard_normal((sites, 2)),
        'beta': beta,
    }
    nematic = FieldElement('nematic', element=NematicDensitySiteElement, space=LatticeSpace(shape), values=values)
    energy = LandauDeGennesEnergyInteraction(
        'energy',
        target='nematic',
        a=a,
        c=c,
        K=K,
        gamma=gamma,
        quadratic=lambda s: s.a * (1 - s.target.rho / critical_density),
    )

    def flux(s: SimpleNamespace) -> list[Expression]:
        """Return the current of the density, its diffusion and its active current."""
        active = tensor_divergence(s.target.Q)
        return [-D * item - current * other for item, other in zip(gradient(s.target.rho), active, strict=True)]

    transport = ConservationInteraction(
        'transport',
        target='nematic',
        state='rho',
        flux=flux,
        limiter=None,
        riemann=None,
    )
    return System(elements=[nematic], interactions=[energy, transport])

ginelli_nematics(n_particles=400, *, size=10.0, v0=0.3, eta=0.08, seed=0)

Return polar particles with nematic alignment, which order along a common axis but move both ways along it.

At every step particle i takes the direction of sum_j sign(cos(theta_j - theta_i)) exp(i theta_j) over the particles within the unit distance, itself included, so it aligns with the axis of every neighbour on the side nearer its own heading, plus a uniform angle in [-eta pi, eta pi], and moves v0 along it. At low noise the particles share an axis, with as many moving one way as the other: the nematic order parameter |<exp(2 i theta)>| grows towards 1 while the polar order parameter stays small.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	400
size	float	The side of the periodic box, in units of the interaction radius.	10.0
v0	float	The distance a particle moves in one step.	0.3
eta	float	The amplitude of the noise, as a fraction of pi.	0.08
seed	int	The seed of the uniform initial positions and headings.	0

Returns:

Name	Type	Description
system	System	The system, with the positions flock.x, the headings flock.theta, and the observables polarization and nematic_order.

Examples:

>>> from skcomplex.systems.physics.active_matter import ginelli_nematics
>>> ginelli_nematics(10).elements[0].values['eta']
0.08

Source code in src/skcomplex/systems/physics/active_matter/_flocks.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'flocking', 'nematic', 'discrete time', 'neighborhood'),
    references=('Ginelli et al. 2010', 'Chaté 2020'),
    defaults={'until': 300.0, 'dt': 1.0},
    observables=lambda system: _orders('flock'),
)
def ginelli_nematics(
    n_particles: int = 400,
    *,
    size: float = 10.0,
    v0: float = 0.3,
    eta: float = 0.08,
    seed: int = 0,
) -> System:
    """Return polar particles with nematic alignment, which order along a common axis but move both ways along it.

    At every step particle `i` takes the direction of `sum_j sign(cos(theta_j - theta_i)) exp(i theta_j)` over the
    particles within the unit distance, itself included, so it aligns with the axis of every neighbour on the side
    nearer its own heading, plus a uniform angle in `[-eta pi, eta pi]`, and moves `v0` along it. At low noise the
    particles share an axis, with as many moving one way as the other: the nematic order parameter
    `|<exp(2 i theta)>|` grows towards 1 while the polar order parameter stays small.

    Args:
        n_particles:
            The number of particles.

        size:
            The side of the periodic box, in units of the interaction radius.

        v0:
            The distance a particle moves in one step.

        eta:
            The amplitude of the noise, as a fraction of `pi`.

        seed:
            The seed of the uniform initial positions and headings.

    Returns:
        system:
            The system, with the positions `flock.x`, the headings `flock.theta`, and the observables `polarization`
            and `nematic_order`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import ginelli_nematics
        >>> ginelli_nematics(10).elements[0].values['eta']
        0.08
    """
    alignment = HeadingAlignmentInteraction('alignment', target='flock', nematic=True)
    return _flock(NematicParticleElement, n_particles, size, seed, v0=v0, alignment=alignment, eta=eta)

gregoire_chate(n_particles=200, *, size=5.0, v0=0.05, eta=0.1, alpha=1.0, beta=0.5, seed=0)

Return the Grégoire-Chaté system: aligning particles that also attract and repel, forming cohesive flocks.

At every step particle i takes the direction of alpha sum_j exp(i theta_j) + beta sum_j f(r_ij) e_ij + eta n_i xi, over the n_i particles within the unit distance, where e_ij points from i to j, xi is a random unit vector, and the force f is a strong repulsion below r_c = 0.2, zero at the preferred distance r_e = 0.5, and attractive up to r_a = 0.8. The particle then moves v0 along the new heading. At low noise and moderate cohesion the particles gather into cohesive, aligned flocks, and the polar order parameter grows towards 1. A strong cohesion locks them in place, and the order is lost.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	200
size	float	The side of the periodic box, in units of the interaction radius.	5.0
v0	float	The distance a particle moves in one step.	0.05
eta	float	The strength of the vectorial noise.	0.1
alpha	float	The strength of the alignment.	1.0
beta	float	The strength of the cohesion.	0.5
seed	int	The seed of the uniform initial positions and headings.	0

Returns:

Name	Type	Description
system	System	The system, with the positions flock.x, the headings flock.theta, and the observables polarization and nematic_order.

Examples:

>>> from skcomplex.systems.physics.active_matter import gregoire_chate
>>> gregoire_chate(10, beta=2.0).elements[0].values['beta']
2.0

Source code in src/skcomplex/systems/physics/active_matter/_flocks.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'flocking', 'cohesion', 'discrete time', 'neighborhood'),
    references=('Grégoire and Chaté 2004', 'Vicsek and Zafeiris 2012'),
    defaults={'until': 300.0, 'dt': 1.0},
    observables=lambda system: _orders('flock'),
)
def gregoire_chate(
    n_particles: int = 200,
    *,
    size: float = 5.0,
    v0: float = 0.05,
    eta: float = 0.1,
    alpha: float = 1.0,
    beta: float = 0.5,
    seed: int = 0,
) -> System:
    """Return the Grégoire-Chaté system: aligning particles that also attract and repel, forming cohesive flocks.

    At every step particle `i` takes the direction of `alpha sum_j exp(i theta_j) + beta sum_j f(r_ij) e_ij +
    eta n_i xi`, over the `n_i` particles within the unit distance, where `e_ij` points from `i` to `j`, `xi` is a
    random unit vector, and the force `f` is a strong repulsion below `r_c = 0.2`, zero at the preferred distance
    `r_e = 0.5`, and attractive up to `r_a = 0.8`. The particle then moves `v0` along the new heading. At low noise and
    moderate cohesion the particles gather into cohesive, aligned flocks, and the polar order parameter grows towards
    1. A strong cohesion locks them in place, and the order is lost.

    Args:
        n_particles:
            The number of particles.

        size:
            The side of the periodic box, in units of the interaction radius.

        v0:
            The distance a particle moves in one step.

        eta:
            The strength of the vectorial noise.

        alpha:
            The strength of the alignment.

        beta:
            The strength of the cohesion.

        seed:
            The seed of the uniform initial positions and headings.

    Returns:
        system:
            The system, with the positions `flock.x`, the headings `flock.theta`, and the observables `polarization`
            and `nematic_order`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import gregoire_chate
        >>> gregoire_chate(10, beta=2.0).elements[0].values['beta']
        2.0
    """
    return _flock(CohesiveParticleElement, n_particles, size, seed, v0=v0, eta=eta, alpha=alpha, beta=beta)

inertial_active_brownian_particles(n_particles=200, *, speed=3.0, mass=0.5, moment=0.1, D_t=0.1, D_r=1.0, size=50.0, random_state=0)

Return active Brownian particles with inertia, which lag behind their orientation as they speed up and turn.

With unit friction, every particle follows m x'' = -x' + v0 n + sqrt(2 D_t) xi and J theta'' = -theta' + sqrt(2 D_r) xi_r, the inertial active Brownian particles of macroscopic swimmers, flyers, and vibrated grains. The velocity relaxes over the time m, so it trails the orientation, and the mean projection of the velocity on the orientation falls below v0. With rotational inertia alone, the orientation correlation is exp(-D_r (t - J (1 - exp(-t / J)))), which gives the steady projection v0 / m ∫ exp(-t / m) exp(-D_r (t - J (1 - exp(-t / J)))) dt, and v0 / (1 + m D_r) without rotational inertia. The particles start at rest at uniform positions and orientations.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	200
speed	float	The self-propulsion speed v0.	3.0
mass	float	The mass m, in units of the friction.	0.5
moment	float	The moment of inertia J, in units of the rotational friction.	0.1
D_t	float	The translational diffusion coefficient.	0.1
D_r	float	The rotational diffusion coefficient.	1.0
size	float	The side of the square in which the particles start.	50.0
random_state	int | Generator | None	The seed or generator of the initial positions and orientations.	0

Returns:

Name	Type	Description
system	System	The system, with the positions swimmers.x, the velocities swimmers.x_d1, the orientations swimmers.theta, and the angular velocities swimmers.theta_d1.

Examples:

>>> from skcomplex.systems.physics.active_matter import inertial_active_brownian_particles
>>> inertial_active_brownian_particles(10).elements[0].prepare().size_
10

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'self-propulsion', 'inertia', 'second order'),
    references=('Löwen 2020', 'Scholz et al. 2018', 'Bechinger et al. 2016'),
    defaults={'until': 20.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.005)},
)
def inertial_active_brownian_particles(
    n_particles: int = 200,
    *,
    speed: float = 3.0,
    mass: float = 0.5,
    moment: float = 0.1,
    D_t: float = 0.1,
    D_r: float = 1.0,
    size: float = 50.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return active Brownian particles with inertia, which lag behind their orientation as they speed up and turn.

    With unit friction, every particle follows `m x'' = -x' + v0 n + sqrt(2 D_t) xi` and
    `J theta'' = -theta' + sqrt(2 D_r) xi_r`, the inertial active Brownian particles of macroscopic swimmers, flyers,
    and vibrated grains. The velocity relaxes over the time `m`, so it trails the orientation, and the mean projection
    of the velocity on the orientation falls below `v0`. With rotational inertia alone, the orientation correlation is
    `exp(-D_r (t - J (1 - exp(-t / J))))`, which gives the steady projection
    `v0 / m ∫ exp(-t / m) exp(-D_r (t - J (1 - exp(-t / J)))) dt`, and `v0 / (1 + m D_r)` without rotational inertia.
    The particles start at rest at uniform positions and orientations.

    Args:
        n_particles:
            The number of particles.

        speed:
            The self-propulsion speed `v0`.

        mass:
            The mass `m`, in units of the friction.

        moment:
            The moment of inertia `J`, in units of the rotational friction.

        D_t:
            The translational diffusion coefficient.

        D_r:
            The rotational diffusion coefficient.

        size:
            The side of the square in which the particles start.

        random_state:
            The seed or generator of the initial positions and orientations.

    Returns:
        system:
            The system, with the positions `swimmers.x`, the velocities `swimmers.x_d1`, the orientations
            `swimmers.theta`, and the angular velocities `swimmers.theta_d1`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import inertial_active_brownian_particles
        >>> inertial_active_brownian_particles(10).elements[0].prepare().size_
        10
    """
    positions, angles = _swimmers(n_particles, size, random_state)
    group = GroupElement(
        'swimmers',
        element=InertialActiveParticleElement,
        size=n_particles,
        values={
            'x': positions,
            'theta': angles,
            'speed': speed,
            'mass': mass,
            'moment': moment,
            'D_t': D_t,
            'D_r': D_r,
        },
    )
    return System(elements=[group], space=ContinuousSpace(n_dim=2))

keller_segel(shape=(32, 32), *, length=32.0, D=1.0, chi=5.0, Dc=1.0, k=1.0, density=1.0, random_state=0)

Return the Keller-Segel system of cells that diffuse and climb the gradient of a chemical they secrete.

The cell density follows rho_t = D laplacian(rho) - chi div(rho grad(c)) and the chemical c_t = Dc laplacian(c) + rho - k c, on a periodic square. The total number of cells is conserved. The uniform state of density rho_0 is unstable when the chemotactic drive chi rho_0 exceeds D (k + Dc q^2) for the longest wave q = 2 pi / length. The cells then aggregate into dense spots, which in two dimensions collapse in finite time, so a run above the threshold ends before the density blows up.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of sites along both axes of the periodic lattice.	(32, 32)
length	float	The side of the square.	32.0
D	float	The diffusivity of the cells.	1.0
chi	float	The chemotactic sensitivity.	5.0
Dc	float	The diffusivity of the chemical.	1.0
k	float	The decay rate of the chemical.	1.0
density	float	The mean cell density rho_0, perturbed by 1% noise, with the chemical at its uniform value rho_0 / k.	1.0
random_state	int | Generator | None	The seed or generator of the noise.	0

Returns:

Name	Type	Description
system	System	The system, with the cell density colony.rho and the chemical colony.c.

Examples:

>>> from skcomplex.systems.physics.active_matter import keller_segel
>>> keller_segel((8, 8), chi=0.5).elements[0].values['chi']
0.5

Source code in src/skcomplex/systems/physics/active_matter/_chemotaxis.py

@register_system(
    field='physics/active_matter',
    tags=('field', 'chemotaxis', 'aggregation', 'collapse', 'pde'),
    references=('te Vrugt and Wittkowski 2024', 'Keller and Segel 1970', 'Horstmann 2003'),
    defaults={'until': 10.5, 'dt': 0.25},
)
def keller_segel(
    shape: Sequence[int] = (32, 32),
    *,
    length: float = 32.0,
    D: float = 1.0,
    chi: float = 5.0,
    Dc: float = 1.0,
    k: float = 1.0,
    density: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Keller-Segel system of cells that diffuse and climb the gradient of a chemical they secrete.

    The cell density follows `rho_t = D laplacian(rho) - chi div(rho grad(c))` and the chemical
    `c_t = Dc laplacian(c) + rho - k c`, on a periodic square. The total number of cells is conserved. The uniform state
    of density `rho_0` is unstable when the chemotactic drive `chi rho_0` exceeds `D (k + Dc q^2)` for the longest
    wave `q = 2 pi / length`. The cells then aggregate into dense spots, which in two dimensions collapse in finite
    time, so a run above the threshold ends before the density blows up.

    Args:
        shape:
            The number of sites along both axes of the periodic lattice.

        length:
            The side of the square.

        D:
            The diffusivity of the cells.

        chi:
            The chemotactic sensitivity.

        Dc:
            The diffusivity of the chemical.

        k:
            The decay rate of the chemical.

        density:
            The mean cell density `rho_0`, perturbed by 1% noise, with the chemical at its uniform value `rho_0 / k`.

        random_state:
            The seed or generator of the noise.

    Returns:
        system:
            The system, with the cell density `colony.rho` and the chemical `colony.c`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import keller_segel
        >>> keller_segel((8, 8), chi=0.5).elements[0].values['chi']
        0.5
    """
    shape = tuple(shape)
    rng = np.random.default_rng(random_state)
    values = {
        'rho': density * (1 + 0.01 * rng.uniform(-1.0, 1.0, shape)),
        'c': np.full(shape, density / k),
        'D': D,
        'chi': chi,
        'Dc': Dc,
        'k': k,
    }
    space = LatticeSpace(shape, spacing=length / shape[0])
    return System(elements=[FieldElement('colony', element=ChemotaxisElement, space=space, values=values)])

living_crystals(n_particles=100, *, size=50.0, speed=3.0, D_r=1.0, v0=20.0, cutoff=10.0, k=200.0, random_state=0)

Return active colloids of radius 1 with phoretic attraction, which gather into clusters that slowly coarsen.

The colloids swim as active Brownian particles, drift towards each other at the speed v0 / r^2 within the cutoff, and a stiff soft repulsion stands in for the hard-sphere correction at contact.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of colloids.	100
size	float	The side of the periodic box.	50.0
speed	float	The self-propulsion speed.	3.0
D_r	float	The rotational diffusion coefficient.	1.0
v0	float	The phoretic drift speed at the distance of one radius.	20.0
cutoff	float	The distance beyond which the attraction is neglected.	10.0
k	float	The stiffness of the repulsion times the mobility.	200.0
random_state	int | Generator | None	The seed or generator of the initial positions and orientations.	0

Returns:

Name	Type	Description
system	System	The system, with the positions swimmers.x and orientations swimmers.theta.

Examples:

>>> from skcomplex.systems.physics.active_matter import living_crystals
>>> living_crystals(16).interactions[2].v0
20.0

Source code in src/skcomplex/systems/physics/active_matter/_collective.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'clusters', 'phoresis', 'pair interaction'),
    references=('Argun et al. 2021', 'Palacci et al. 2013'),
    defaults={'until': 20.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.002)},
)
def living_crystals(
    n_particles: int = 100,
    *,
    size: float = 50.0,
    speed: float = 3.0,
    D_r: float = 1.0,
    v0: float = 20.0,
    cutoff: float = 10.0,
    k: float = 200.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return active colloids of radius 1 with phoretic attraction, which gather into clusters that slowly coarsen.

    The colloids swim as active Brownian particles, drift towards each other at the speed `v0 / r^2` within the
    cutoff, and a stiff soft repulsion stands in for the hard-sphere correction at contact.

    Args:
        n_particles:
            The number of colloids.

        size:
            The side of the periodic box.

        speed:
            The self-propulsion speed.

        D_r:
            The rotational diffusion coefficient.

        v0:
            The phoretic drift speed at the distance of one radius.

        cutoff:
            The distance beyond which the attraction is neglected.

        k:
            The stiffness of the repulsion times the mobility.

        random_state:
            The seed or generator of the initial positions and orientations.

    Returns:
        system:
            The system, with the positions `swimmers.x` and orientations `swimmers.theta`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import living_crystals
        >>> living_crystals(16).interactions[2].v0
        20.0
    """
    system = _active_disks(n_particles, size, speed, D_r, random_state)
    repulsion = SoftSphereInteraction('contact', first='swimmers', second='swimmers', cutoff=2.0, period=size, k=k)
    attraction = PhoreticAttractionInteraction(
        'phoresis',
        first='swimmers',
        second='swimmers',
        cutoff=cutoff,
        period=size,
        v0=v0,
        R=1.0,
    )
    interactions = [*system.interactions, repulsion, attraction]
    return System(elements=system.elements, interactions=interactions, space=system.space)

motility_induced_phase_separation(n_particles=400, *, packing_fraction=0.5, speed=1.0, D_r=0.01, k=50.0, random_state=0)

Return repulsive active Brownian disks of unit diameter, which phase separate into clusters at high speed.

Disks that collide head on block each other until their orientations diffuse away, so fast swimmers with slow rotational diffusion, a large Péclet number v / D_r, accumulate where they are slow and form dense clusters in a dilute gas.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of disks.	400
packing_fraction	float	The fraction of the periodic box the disks cover.	0.5
speed	float	The self-propulsion speed.	1.0
D_r	float	The rotational diffusion coefficient.	0.01
k	float	The stiffness of the repulsion times the mobility.	50.0
random_state	int | Generator | None	The seed or generator of the initial positions and orientations.	0

Returns:

Name	Type	Description
system	System	The system, with the positions swimmers.x and orientations swimmers.theta.

Examples:

>>> from skcomplex.systems.physics.active_matter import motility_induced_phase_separation
>>> motility_induced_phase_separation(16).interactions[1].k
50.0

Source code in src/skcomplex/systems/physics/active_matter/_collective.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'phase separation', 'pair interaction'),
    references=('Cates and Tailleur 2015', 'Bechinger et al. 2016', 'Fily and Marchetti 2012'),
    defaults={'until': 20.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.005)},
)
def motility_induced_phase_separation(
    n_particles: int = 400,
    *,
    packing_fraction: float = 0.5,
    speed: float = 1.0,
    D_r: float = 0.01,
    k: float = 50.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return repulsive active Brownian disks of unit diameter, which phase separate into clusters at high speed.

    Disks that collide head on block each other until their orientations diffuse away, so fast swimmers with slow
    rotational diffusion, a large Péclet number `v / D_r`, accumulate where they are slow and form dense clusters in a
    dilute gas.

    Args:
        n_particles:
            The number of disks.

        packing_fraction:
            The fraction of the periodic box the disks cover.

        speed:
            The self-propulsion speed.

        D_r:
            The rotational diffusion coefficient.

        k:
            The stiffness of the repulsion times the mobility.

        random_state:
            The seed or generator of the initial positions and orientations.

    Returns:
        system:
            The system, with the positions `swimmers.x` and orientations `swimmers.theta`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import motility_induced_phase_separation
        >>> motility_induced_phase_separation(16).interactions[1].k
        50.0
    """
    size = math.sqrt(n_particles * math.pi / 4 / packing_fraction)
    system = _active_disks(n_particles, size, speed, D_r, random_state)
    repulsion = SoftSphereInteraction('contact', first='swimmers', second='swimmers', cutoff=1.0, period=size, k=k)
    return System(elements=system.elements, interactions=[*system.interactions, repulsion], space=system.space)

nonreciprocal_vicsek(n_movers=400, *, size=5.0, J_ab=1.0, J_ba=-1.0, eta=0.3, random_state=1)

Return two species of Vicsek movers in a periodic box whose alignment with each other is nonreciprocal.

Half the movers are species A and half species B, each aligning with its own species. With J_ab = J_ba the two species flock together. With J_ab = -J_ba, A aligns with B while B antialigns with A, so neither settles: the mean headings of the species lock at a quarter turn from each other and rotate together, the chiral phase of Fruchart et al. (2021), which breaks the time-reversal symmetry that reciprocal alignment keeps.

Parameters:

Name	Type	Description	Default
n_movers	int	The number of movers, an even number.	400
size	float	The side of the periodic box.	5.0
J_ab	float	The alignment of A with B.	1.0
J_ba	float	The alignment of B with A.	-1.0
eta	float	The width of the uniform noise of the headings.	0.3
random_state	int | Generator | None	The seed or generator of the initial positions and headings.	1

Returns:

Name	Type	Description
system	System	The system, with the mean headings heading_a and heading_b of the two species.

Examples:

>>> from skcomplex.systems.physics.active_matter import nonreciprocal_vicsek
>>> nonreciprocal_vicsek(20).elements[0].prepare().size_
20

Source code in src/skcomplex/systems/physics/active_matter/_nonreciprocal.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'flocking', 'nonreciprocal', 'vicsek', 'chiral phase'),
    references=('Fruchart, Hanai, Littlewood, and Vitelli 2021', 'te Vrugt and Wittkowski 2024'),
    defaults={'until': 300.0, 'dt': 1.0},
    observables=lambda system: [
        ObservableSpec('heading_a', _heading_a, group='movers'),
        ObservableSpec('heading_b', _heading_b, group='movers'),
    ],
)
def nonreciprocal_vicsek(
    n_movers: int = 400,
    *,
    size: float = 5.0,
    J_ab: float = 1.0,
    J_ba: float = -1.0,
    eta: float = 0.3,
    random_state: int | np.random.Generator | None = 1,
) -> System:
    """Return two species of Vicsek movers in a periodic box whose alignment with each other is nonreciprocal.

    Half the movers are species A and half species B, each aligning with its own species. With `J_ab = J_ba` the two
    species flock together. With `J_ab = -J_ba`, A aligns with B while B antialigns with A, so neither settles: the
    mean headings of the species lock at a quarter turn from each other and rotate together, the chiral phase of
    Fruchart et al. (2021), which breaks the time-reversal symmetry that reciprocal alignment keeps.

    Args:
        n_movers:
            The number of movers, an even number.

        size:
            The side of the periodic box.

        J_ab:
            The alignment of A with B.

        J_ba:
            The alignment of B with A.

        eta:
            The width of the uniform noise of the headings.

        random_state:
            The seed or generator of the initial positions and headings.

    Returns:
        system:
            The system, with the mean headings `heading_a` and `heading_b` of the two species.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import nonreciprocal_vicsek
        >>> nonreciprocal_vicsek(20).elements[0].prepare().size_
        20
    """
    rng = np.random.default_rng(random_state)
    species = (np.arange(n_movers) % 2).astype(float)
    movers = GroupElement(
        'movers',
        element=NonreciprocalMoverElement,
        size=n_movers,
        values={
            'x': rng.uniform(0, size, (n_movers, 2)),
            'theta': rng.uniform(-math.pi, math.pi, n_movers),
            'species': species,
            'J_ab': J_ab,
            'J_ba': J_ba,
            'eta': eta,
        },
    )
    return System(
        elements=[movers],
        interactions=[PropulsionStepInteraction('propulsion', target='movers', speed=0.05)],
        space=ContinuousSpace(extent=(size, size), boundary='periodic'),
        time=DiscreteTime(),
    )

run_and_tumble(n_agents=100, *, speed=0.1, tumble_rate=0.1, size=10.0, random_state=0)

Return run-and-tumble agents in a periodic square, which diffuse with v^2 (1 + q) / (4 (1 - q)) per step.

The persistence q = exp(-tumble_rate) is the probability that an agent keeps its direction for one step.

Parameters:

Name	Type	Description	Default
n_agents	int	The number of agents.	100
speed	float	The distance run per step.	0.1
tumble_rate	float	The rate of tumbles per step.	0.1
size	float	The side of the periodic square.	10.0
random_state	int | Generator | None	The seed or generator of the initial positions and directions.	0

Returns:

Name	Type	Description
system	System	The system, with the agent group bacteria.

Examples:

>>> from skcomplex.systems.physics.active_matter import run_and_tumble
>>> run_and_tumble(10).elements[0].prepare().size_
10

Source code in src/skcomplex/systems/physics/active_matter/_particles.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'agents', 'bacteria', 'self-propulsion'),
    references=('Berg 2004', 'Bechinger et al. 2016'),
    defaults={'until': 100, 'dt': 1},
)
def run_and_tumble(
    n_agents: int = 100,
    *,
    speed: float = 0.1,
    tumble_rate: float = 0.1,
    size: float = 10.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return run-and-tumble agents in a periodic square, which diffuse with `v^2 (1 + q) / (4 (1 - q))` per step.

    The persistence `q = exp(-tumble_rate)` is the probability that an agent keeps its direction for one step.

    Args:
        n_agents:
            The number of agents.

        speed:
            The distance run per step.

        tumble_rate:
            The rate of tumbles per step.

        size:
            The side of the periodic square.

        random_state:
            The seed or generator of the initial positions and directions.

    Returns:
        system:
            The system, with the agent group `bacteria`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import run_and_tumble
        >>> run_and_tumble(10).elements[0].prepare().size_
        10
    """
    positions, angles = _swimmers(n_agents, size, random_state)
    group = GroupElement(
        'bacteria',
        element=TumblerElement,
        size=n_agents,
        values={'theta': angles, 'speed': speed, 'tumble_rate': tumble_rate, 'width': size},
        location=positions,
    )
    return System(elements=[group], space=ContinuousSpace(extent=(size, size), boundary='periodic'))

squirmers(n_squirmers=27, *, B1=1.5, beta=-3.0, radius=1.0, viscosity=1.0, spacing=5.0, random_state=0)

Return spherical squirmers that swim at 2 B1 / 3 and stir each other as pushers or pullers.

The squirmers start on a jittered cubic lattice of the spacing in an unbounded fluid, with random orientations. beta = B2 / B1 sets the stresslet: pushers (beta < 0) draw their neighbours in from the sides and pullers (beta > 0) push them away.

Parameters:

Name	Type	Description	Default
n_squirmers	int	The number of squirmers.	27
B1	float	The first slip mode.	1.5
beta	float	The ratio B2 / B1 of the slip modes.	-3.0
radius	float	The radius of the squirmers.	1.0
viscosity	float	The viscosity of the fluid.	1.0
spacing	float	The spacing of the starting lattice.	5.0
random_state	int | Generator | None	The seed or generator of the starting jitter and orientations.	0

Returns:

Name	Type	Description
system	System	The system, with the positions squirmers.x and orientations squirmers.p.

Examples:

>>> from skcomplex.systems.physics.active_matter import squirmers
>>> squirmers(8, beta=2.0).interactions[0].B2
3.0

Source code in src/skcomplex/systems/physics/active_matter/_squirmers.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'microswimmers', 'squirmers', 'hydrodynamics', 'stokes flow', 'pushers', 'pullers'),
    references=(
        'Lighthill 1952',
        'Blake 1971',
        'Ishikawa, Simmonds, and Pedley 2006',
        'te Vrugt and Wittkowski 2024',
    ),
    defaults={'until': 40.0, 'dt': 0.5},
)
def squirmers(
    n_squirmers: int = 27,
    *,
    B1: float = 1.5,
    beta: float = -3.0,
    radius: float = 1.0,
    viscosity: float = 1.0,
    spacing: float = 5.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return spherical squirmers that swim at `2 B1 / 3` and stir each other as pushers or pullers.

    The squirmers start on a jittered cubic lattice of the spacing in an unbounded fluid, with random orientations.
    `beta = B2 / B1` sets the stresslet: pushers (`beta < 0`) draw their neighbours in from the sides and pullers
    (`beta > 0`) push them away.

    Args:
        n_squirmers:
            The number of squirmers.

        B1:
            The first slip mode.

        beta:
            The ratio `B2 / B1` of the slip modes.

        radius:
            The radius of the squirmers.

        viscosity:
            The viscosity of the fluid.

        spacing:
            The spacing of the starting lattice.

        random_state:
            The seed or generator of the starting jitter and orientations.

    Returns:
        system:
            The system, with the positions `squirmers.x` and orientations `squirmers.p`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import squirmers
        >>> squirmers(8, beta=2.0).interactions[0].B2
        3.0
    """
    rng = np.random.default_rng(random_state)
    side = math.ceil(n_squirmers ** (1 / 3) - 1e-9)
    grid = np.stack(np.meshgrid(*[np.arange(side)] * 3, indexing='ij'), axis=-1).reshape(-1, 3)[:n_squirmers]
    positions = spacing * grid + rng.uniform(-0.1, 0.1, size=(n_squirmers, 3)) * spacing
    orientations = rng.standard_normal((n_squirmers, 3))
    orientations /= np.linalg.norm(orientations, axis=1, keepdims=True)
    group = GroupElement(
        'squirmers',
        element=SquirmerElement,
        size=n_squirmers,
        values={'x': positions, 'p': orientations},
    )
    slip = SquirmerSlipInteraction(
        'slip',
        target='squirmers',
        B1=B1,
        B2=beta * B1,
        radius=radius,
        viscosity=viscosity,
    )
    return System(elements=[group], interactions=[slip], space=ContinuousSpace(n_dim=3))

szabo_cells(n_cells=200, *, density=1.0, eta=0.6, v0=1.0, mu=1.0, tau=1.0, R_eq=5.0 / 6.0, R_0=1.0, F_rep=30.0, F_adh=0.75, random_state=0)

Return the cells of Szabó et al. (2006) in a periodic box, which move together above a critical density.

A cell moves at x' = v0 n + mu F, where F sums the pair forces of the neighbours within R_0, the diameter of a cell, and its heading turns towards its velocity, theta' = arcsin((n x v / |v|) . e_z) / tau + xi. No cell reads the headings of the others, and adhesion alone orders them. With eta = 0.6 the transition lies at the number density rho_c = 0.18 / R_0^2.

Parameters:

Name	Type	Description	Default
n_cells	int	The number of cells.	200
density	float	The number of cells per unit area.	1.0
eta	float	The strength of the angular noise.	0.6
v0	float	The speed of the motility.	1.0
mu	float	The mobility.	1.0
tau	float	The relaxation time of the heading.	1.0
R_eq	float	The equilibrium distance of a pair.	5.0 / 6.0
R_0	float	The distance beyond which a pair exerts no force.	1.0
F_rep	float	The largest repulsive force.	30.0
F_adh	float	The largest adhesive force.	0.75
random_state	int | Generator | None	The seed or generator of the initial positions and headings.	0

Returns:

Name	Type	Description
system	System	The system, with the positions cells.x, the headings cells.theta, and the observables mean_heading and polarization.

Examples:

>>> from skcomplex.systems.physics.active_matter import szabo_cells
>>> szabo_cells(16, density=0.5).space.extent
(5.656854249492381, 5.656854249492381)

Source code in src/skcomplex/systems/physics/active_matter/_cells.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'flocking', 'tissue', 'neighbours'),
    references=('Szabó et al. 2006', 'Vicsek and Zafeiris 2012'),
    defaults={'until': 60.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.01), 'backend': NumbaBackend()},
    observables=lambda system: [
        ObservableSpec('mean_heading', _heading, group='cells'),
        ObservableSpec('polarization', _polarization),
    ],
)
def szabo_cells(
    n_cells: int = 200,
    *,
    density: float = 1.0,
    eta: float = 0.6,
    v0: float = 1.0,
    mu: float = 1.0,
    tau: float = 1.0,
    R_eq: float = 5.0 / 6.0,
    R_0: float = 1.0,
    F_rep: float = 30.0,
    F_adh: float = 0.75,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the cells of Szabó et al. (2006) in a periodic box, which move together above a critical density.

    A cell moves at `x' = v0 n + mu F`, where `F` sums the pair forces of the neighbours within `R_0`, the diameter of a
    cell, and its heading turns towards its velocity, `theta' = arcsin((n x v / |v|) . e_z) / tau + xi`. No cell reads
    the headings of the others, and adhesion alone orders them. With `eta = 0.6` the transition lies at the number
    density `rho_c = 0.18 / R_0^2`.

    Args:
        n_cells:
            The number of cells.

        density:
            The number of cells per unit area.

        eta:
            The strength of the angular noise.

        v0:
            The speed of the motility.

        mu:
            The mobility.

        tau:
            The relaxation time of the heading.

        R_eq:
            The equilibrium distance of a pair.

        R_0:
            The distance beyond which a pair exerts no force.

        F_rep:
            The largest repulsive force.

        F_adh:
            The largest adhesive force.

        random_state:
            The seed or generator of the initial positions and headings.

    Returns:
        system:
            The system, with the positions `cells.x`, the headings `cells.theta`, and the observables `mean_heading` and
            `polarization`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import szabo_cells
        >>> szabo_cells(16, density=0.5).space.extent
        (5.656854249492381, 5.656854249492381)
    """
    size = math.sqrt(n_cells / density)
    rng = np.random.default_rng(random_state)
    cells = GroupElement(
        'cells',
        element=MotileCellElement,
        size=n_cells,
        values={
            'x': rng.uniform(0.0, size, (n_cells, 2)),
            'theta': rng.uniform(0.0, 2 * math.pi, n_cells),
            'v0': v0,
            'mu': mu,
            'tau': tau,
            'eta': eta,
            'R_eq': R_eq,
            'R_0': R_0,
            'F_rep': F_rep,
            'F_adh': F_adh,
        },
        visual=VisualSpec('points', {'position': 'x', 'radius': R_0 / 2, 'color_by': 'theta'}),
    )
    motion = RateInteraction('motion', target='cells', state='x', rate=functools.partial(_velocity, reach=R_0))
    turning = RateInteraction(
        'turning',
        target='cells',
        state='theta',
        rate=functools.partial(_turning, reach=R_0),
        noise=_angular_noise,
    )
    return System(
        elements=[cells],
        interactions=[motion, turning],
        space=ContinuousSpace(extent=(size, size), boundary='periodic'),
    )

toner_tu(shape=(32, 32), *, rho0=1.0, rho_c=0.5, a=1.0, beta=1.0, lam=1.0, sigma=0.5, D=0.5, spacing=1.0, random_state=0)

Return the Toner-Tu hydrodynamics of a polar flock on a periodic lattice, from a disordered velocity field.

Above the critical density the velocity grows and aligns until the whole fluid moves together with the speed sqrt(a (rho0 - rho_c) / beta).

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of sites along each axis.	(32, 32)
rho0	float	The mean density.	1.0
rho_c	float	The critical density of the ordering transition.	0.5
a	float	The rate of the ordering instability per unit density.	1.0
beta	float	The saturation of the speed.	1.0
lam	float	The strength of the self-advection.	1.0
sigma	float	The pressure per unit density.	0.5
D	float	The viscosity of the velocity.	0.5
spacing	float	The distance between sites.	1.0
random_state	int | Generator | None	The seed or generator of the initial velocities.	0

Returns:

Name	Type	Description
system	System	The system, with the fields flock.rho, flock.vx, and flock.vy.

Examples:

>>> from skcomplex.systems.physics.active_matter import toner_tu
>>> toner_tu((8, 8)).elements[0].prepare().size_
64

Source code in src/skcomplex/systems/physics/active_matter/_collective.py

@register_system(
    field='physics/active_matter',
    tags=('field', 'flocking', 'hydrodynamics', 'pde'),
    references=('Toner and Tu 1995', 'Marchetti et al. 2013'),
    defaults={'until': 20.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.02)},
)
def toner_tu(
    shape: Sequence[int] = (32, 32),
    *,
    rho0: float = 1.0,
    rho_c: float = 0.5,
    a: float = 1.0,
    beta: float = 1.0,
    lam: float = 1.0,
    sigma: float = 0.5,
    D: float = 0.5,
    spacing: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Toner-Tu hydrodynamics of a polar flock on a periodic lattice, from a disordered velocity field.

    Above the critical density the velocity grows and aligns until the whole fluid moves together with the speed
    `sqrt(a (rho0 - rho_c) / beta)`.

    Args:
        shape:
            The number of sites along each axis.

        rho0:
            The mean density.

        rho_c:
            The critical density of the ordering transition.

        a:
            The rate of the ordering instability per unit density.

        beta:
            The saturation of the speed.

        lam:
            The strength of the self-advection.

        sigma:
            The pressure per unit density.

        D:
            The viscosity of the velocity.

        spacing:
            The distance between sites.

        random_state:
            The seed or generator of the initial velocities.

    Returns:
        system:
            The system, with the fields `flock.rho`, `flock.vx`, and `flock.vy`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import toner_tu
        >>> toner_tu((8, 8)).elements[0].prepare().size_
        64
    """
    shape = tuple(shape)
    rng = np.random.default_rng(random_state)
    values = {
        'rho': np.full(shape, rho0),
        'vx': rng.normal(0.0, 0.1, shape),
        'vy': rng.normal(0.0, 0.1, shape) + 0.05,
    }
    parameters = {'a': a, 'rho_c': rho_c, 'beta': beta, 'lam': lam, 'sigma': sigma, 'D': D}
    space = LatticeSpace(shape, spacing=spacing)
    return System(
        elements=[FieldElement('flock', element=PolarFluidElement, space=space, values={**values, **parameters})],
    )

vicsek(n_particles=400, *, size=10.0, v0=0.5, eta=0.1, seed=0, neighbors=None, vision=None, delay=0, memory=5)

Return the Vicsek system: particles that move at a constant speed along the mean heading of their neighbours.

At every step particle i takes the direction of the sum of the headings of the particles within the unit distance, itself included, plus a uniform angle in [-eta pi, eta pi], and moves v0 along it in a periodic box. At low noise and high density the particles align into a flock, and the polar order parameter |<exp(i theta)>| grows from about 1 / sqrt(N) to near 1. Raising the noise or lowering the density destroys the order.

Two variants of Argun et al. (2021, Section 8.3) make the interaction non-metric or non-reciprocal. With neighbors, a particle aligns with its neighbors nearest particles and itself at any distance, a topological interaction that need not be reciprocal, and the clusters grow with neighbors. With vision, a particle sees only the particles within the unit distance whose direction from it makes an angle of at most vision with its heading, a vision cone that breaks reciprocity. The two variants do not combine, since the nearest neighbourhood of the core has no cone.

With a positive delay h, the sensory delay of Section 8.2, a particle aligns with the headings its neighbours had h steps ago, theta_j(n + 1) = <theta_k(n - h)> + noise. Argun et al. report that a short delay stabilizes the clusters, and a long delay hinders them. A negative delay is a prediction: the headings of the neighbours |h| steps ahead are extrapolated by a least-squares line through their last memory headings, and it hinders clustering. Every particle then keeps its last 25 headings, so the delay lies in [-25, 25], and before the start they equal the initial heading.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	400
size	float	The side of the periodic box, in units of the interaction radius.	10.0
v0	float	The distance a particle moves in one step.	0.5
eta	float	The amplitude of the noise, as a fraction of pi.	0.1
seed	int	The seed of the uniform initial positions and headings.	0
neighbors	int | None	The number of nearest particles a particle aligns with, or None for the particles within the unit distance.	None
vision	float | None	The half-angle of the vision cone in (0, pi], or None for the whole disk.	None
delay	int	The delay of the alignment in steps, negative for a prediction.	0
memory	int	The number of past headings a prediction fits, from 2 to 26.	5

Returns:

Name	Type	Description
system	System	The system, with the positions flock.x, the headings flock.theta, and the observables polarization and nematic_order.

Raises:

Type	Description
ParameterError	Both neighbors and vision are given, or the delay or the memory exceeds the 25 kept headings.

Examples:

>>> from skcomplex.systems.physics.active_matter import vicsek
>>> vicsek(10).elements[0].size
10
>>> type(vicsek(10, neighbors=4).interactions[1].neighborhood).__name__
'NearestNeighborhood'

Source code in src/skcomplex/systems/physics/active_matter/_flocks.py

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'flocking', 'alignment', 'discrete time', 'neighborhood'),
    references=('Vicsek et al. 1995', 'Argun et al. 2021', 'Jensen 2023', 'Vicsek and Zafeiris 2012'),
    defaults={'until': 200.0, 'dt': 1.0},
    observables=lambda system: _orders('flock'),
)
def vicsek(
    n_particles: int = 400,
    *,
    size: float = 10.0,
    v0: float = 0.5,
    eta: float = 0.1,
    seed: int = 0,
    neighbors: int | None = None,
    vision: float | None = None,
    delay: int = 0,
    memory: int = 5,
) -> System:
    """Return the Vicsek system: particles that move at a constant speed along the mean heading of their neighbours.

    At every step particle `i` takes the direction of the sum of the headings of the particles within the unit
    distance, itself included, plus a uniform angle in `[-eta pi, eta pi]`, and moves `v0` along it in a periodic box.
    At low noise and high density the particles align into a flock, and the polar order parameter
    `|<exp(i theta)>|` grows from about `1 / sqrt(N)` to near 1. Raising the noise or lowering the density destroys
    the order.

    Two variants of Argun et al. (2021, Section 8.3) make the interaction non-metric or non-reciprocal. With
    `neighbors`, a particle aligns with its `neighbors` nearest particles and itself at any distance, a topological
    interaction that need not be reciprocal, and the clusters grow with `neighbors`. With `vision`, a particle sees only
    the particles within the unit distance whose direction from it makes an angle of at most `vision` with its heading,
    a vision cone that breaks reciprocity. The two variants do not combine, since the nearest neighbourhood of the core
    has no cone.

    With a positive `delay` `h`, the sensory delay of Section 8.2, a particle aligns with the headings its neighbours
    had `h` steps ago, `theta_j(n + 1) = <theta_k(n - h)> + noise`. Argun et al. report that a short delay stabilizes
    the clusters, and a long delay hinders them. A negative delay is a prediction: the headings of the neighbours `|h|`
    steps ahead are extrapolated by a least-squares line through their last `memory` headings, and it hinders
    clustering. Every particle then keeps its last 25 headings, so the delay lies in `[-25, 25]`, and before the start
    they equal the initial heading.

    Args:
        n_particles:
            The number of particles.

        size:
            The side of the periodic box, in units of the interaction radius.

        v0:
            The distance a particle moves in one step.

        eta:
            The amplitude of the noise, as a fraction of `pi`.

        seed:
            The seed of the uniform initial positions and headings.

        neighbors:
            The number of nearest particles a particle aligns with, or `None` for the particles within the unit
            distance.

        vision:
            The half-angle of the vision cone in `(0, pi]`, or `None` for the whole disk.

        delay:
            The delay of the alignment in steps, negative for a prediction.

        memory:
            The number of past headings a prediction fits, from 2 to 26.

    Returns:
        system:
            The system, with the positions `flock.x`, the headings `flock.theta`, and the observables `polarization`
            and `nematic_order`.

    Raises:
        ParameterError:
            Both `neighbors` and `vision` are given, or the delay or the memory exceeds the 25 kept headings.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import vicsek
        >>> vicsek(10).elements[0].size
        10
        >>> type(vicsek(10, neighbors=4).interactions[1].neighborhood).__name__
        'NearestNeighborhood'
    """
    if neighbors is not None and vision is not None:
        error_msg = 'A Vicsek flock aligns with its nearest neighbours or within a vision cone, not both.'
        raise ParameterError(error_msg)
    if abs(delay) > MEMORY or not FIT_POINTS <= memory <= MEMORY + 1:
        error_msg = (
            f'The delay {delay!r} and the memory {memory!r} need |delay| <= {MEMORY} and 2 <= memory <= {MEMORY + 1}.'
        )
        raise ParameterError(error_msg)
    neighborhood: BaseNeighborhood | None = None
    if neighbors is not None:
        neighborhood = NearestNeighborhood(neighbors, center=True)
    if vision is not None:
        neighborhood = ConeNeighborhood(1.0, vision, direction='theta', center=True)
    alignment = HeadingAlignmentInteraction(
        'alignment',
        target='flock',
        neighborhood=neighborhood,
        latency=delay,
        memory=memory,
    )
    element = VicsekParticleElement if delay == 0 else DelayedVicsekParticleElement
    return _flock(element, n_particles, size, seed, v0=v0, alignment=alignment, eta=eta)