Active matter

Active matter: self-propelled particles, flocking, polar fluids, chemotaxis, Active Model B+, and their orders.

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.models.physics.active_matter import ActiveParticleElement
>>> ActiveParticleElement('swimmer', D_r=2.0).D_r
2.0

Source code in src/skcomplex/models/physics/active_matter/_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 Model 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.models.physics.active_matter import ActiveScalarElement
>>> ActiveScalarElement('site').zeta
4.0

Source code in src/skcomplex/models/physics/active_matter/_model_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.models.physics.active_matter import AngularAlignmentInteraction
>>> align = AngularAlignmentInteraction('align', candidate='flock', field='flock', frequency=0.5, noise=0.5)
>>> align.check().outcome
'conservative'

Source code in src/skcomplex/models/physics/active_matter/_angular_alignment.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.models.physics.active_matter import AngularDiffusionInteraction
>>> AngularDiffusionInteraction('tumble', candidate='flock', angles=8, tumbling=1.0, noise=0.5).check().field

Source code in src/skcomplex/models/physics/active_matter/_angular_diffusion.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.models.physics.active_matter import BinaryCollisionsInteraction
>>> BinaryCollisionsInteraction('collisions', target='birds').d0
0.5

Source code in src/skcomplex/models/physics/active_matter/_binary_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(model)

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

Parameters:

Name	Type	Description	Default
model	ModelView	The model.	required

Source code in src/skcomplex/models/physics/active_matter/_binary_collisions.py

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

    Args:
        model:
            The model.
    """
    agents = model.agents(self.target, FlockerView)
    size = np.asarray(model.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 = (model.parameter(name) for name in ('speed', 'lam', 'sigma', 'd0'))
    substeps = max(1, round(model.parameter('substeps')))
    tau, rng = 1.0 / substeps, model.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 model, 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.models.physics.active_matter import ChemotaxisElement
>>> ChemotaxisElement('site').chi
5.0

Source code in src/skcomplex/models/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.models.physics.active_matter import CohesiveParticleElement
>>> CohesiveParticleElement('bird').r_e
0.5

Source code in src/skcomplex/models/physics/active_matter/_cohesive_particle.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.models.physics.active_matter import ContractileGelSiteElement
>>> ContractileGelSiteElement('site').T0
4.0

Source code in src/skcomplex/models/physics/active_matter/_contractile_gel_site.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 model.	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.models.physics.active_matter import CuckerSmaleAlignmentInteraction
>>> CuckerSmaleAlignmentInteraction('align', first='birds', second='birds', beta=0.25).beta
0.25

Source code in src/skcomplex/models/physics/active_matter/_cucker_smale_alignment.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

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 model 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.models.physics.active_matter import DepotParticleElement
>>> DepotParticleElement('walker').q
10.0

Source code in src/skcomplex/models/physics/active_matter/_depot_particle.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.models.physics.active_matter import FlockerElement
>>> FlockerElement('bird', theta=1.0).theta
1.0

Source code in src/skcomplex/models/physics/active_matter/_flocker.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 model.	required

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

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

HeadingAlignmentInteraction(label, *, target, nematic=False)

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.

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

Examples:

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

Source code in src/skcomplex/models/physics/active_matter/_heading_alignment.py

def __init__(self, label: str, *, target: Annotated[str, RoleSpec()], nematic: bool = False) -> 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),
    )
    self.nematic = nematic

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.models.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/models/physics/active_matter/_landau_de_gennes_energy.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

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 model.

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.models.physics.active_matter import NematicDensitySiteElement
>>> NematicDensitySiteElement('site').beta
0.5

Source code in src/skcomplex/models/physics/active_matter/_nematic_density_site.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.models.physics.active_matter import NematicGelSiteElement
>>> NematicGelSiteElement('site').zeta
0.5

Source code in src/skcomplex/models/physics/active_matter/_nematic_gel_site.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.models.physics.active_matter import NematicParticleElement
>>> NematicParticleElement('rod', eta=0.2).eta
0.2

Source code in src/skcomplex/models/physics/active_matter/_nematic_particle.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.models.physics.active_matter import NonreciprocalMoverElement
>>> NonreciprocalMoverElement('mover', J_ba=-1.0).J_ba
-1.0

Source code in src/skcomplex/models/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.models.physics.active_matter import OrnsteinUhlenbeckParticleElement
>>> OrnsteinUhlenbeckParticleElement('swimmer').tau
1.0

Source code in src/skcomplex/models/physics/active_matter/_ornstein_uhlenbeck_particle.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.models.physics.active_matter import PhoreticAttractionInteraction
>>> PhoreticAttractionInteraction('phoresis', first='colloids', second='colloids', v0=20.0).v0
20.0

Source code in src/skcomplex/models/physics/active_matter/_phoretic_attraction.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.models.physics.active_matter import PolarAlignmentInteraction
>>> PolarAlignmentInteraction('align', first='swimmers', second='swimmers', gamma=0.5).gamma
0.5

Source code in src/skcomplex/models/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.models.physics.active_matter import PolarFluidElement
>>> PolarFluidElement('site').rho_c
0.5

Source code in src/skcomplex/models/physics/active_matter/_polar_fluid.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.models.physics.active_matter import PolarGelSiteElement
>>> PolarGelSiteElement('site').xi
0.0

Source code in src/skcomplex/models/physics/active_matter/_polar_gel_site.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.models.physics.active_matter import PropulsionStepInteraction
>>> step = PropulsionStepInteraction('propulsion', target='flock').check().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/models/physics/active_matter/_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.models.physics.active_matter import SelfPropulsionInteraction
>>> propulsion = SelfPropulsionInteraction('propulsion', target='swimmers').check().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/models/physics/active_matter/_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

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.models.physics.active_matter import TumblerElement
>>> TumblerElement('bacterium').speed
0.1

Source code in src/skcomplex/models/physics/active_matter/_tumbler.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, model)

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

Parameters:

Name	Type	Description	Default
agent	TumblerView	The agent.	required
model	ModelView	The model.	required

Source code in src/skcomplex/models/physics/active_matter/_tumbler.py

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

    Args:
        agent:
            The agent.

        model:
            The model.
    """
    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 model.rng.random() < -math.expm1(-agent.tumble_rate):
        agent.theta = float(model.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 model.	required

Source code in src/skcomplex/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.models.physics.active_matter import VicsekParticleElement
>>> VicsekParticleElement('bird').eta
0.1

Source code in src/skcomplex/models/physics/active_matter/_vicsek_particle.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 model, with the positions swimmers.x and orientations swimmers.theta.

Examples:

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

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

@register_model(
    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 model, with the positions `swimmers.x` and orientations `swimmers.theta`.

    Examples:
        >>> from skcomplex.models.physics.active_matter import active_brownian_particles
        >>> active_brownian_particles(10).elements[0].check().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 model, 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.models.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/models/physics/active_matter/_nematics.py

@register_model(
    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 model, 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.models.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_model_b(shape=(64, 64), *, mean=0.4, kappa=1.0, lam=1.0, zeta=4.0, D=0.05, random_state=0)

Return Active Model 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 Model 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 model, with the field fluid.phi.

Examples:

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

Source code in src/skcomplex/models/physics/active_matter/_model_b.py

@register_model(
    field='physics/active_matter',
    tags=('field', 'active matter', 'phase separation', 'microphase separation', 'conserved noise', 'model 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_model_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 Model 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 Model 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 model, with the field `fluid.phi`.

    Examples:
        >>> from skcomplex.models.physics.active_matter import active_model_b
        >>> active_model_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_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 model, with the positions swimmers.x and propulsions swimmers.u.

Examples:

>>> from skcomplex.models.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/models/physics/active_matter/_particles.py

@register_model(
    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 model, with the positions `swimmers.x` and propulsions `swimmers.u`.

    Examples:
        >>> from skcomplex.models.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_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 model, with the density gel.c and the displacement gel.u.

Examples:

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

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

@register_model(
    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 model, with the density `gel.c` and the displacement `gel.u`.

    Examples:
        >>> from skcomplex.models.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 model, with the positions swimmers.x, the headings swimmers.theta, and the observables mean_heading and polarization.

Examples:

>>> from skcomplex.models.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/models/physics/active_matter/_alignment.py

@register_model(
    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 model, with the positions `swimmers.x`, the headings `swimmers.theta`, and the observables
            `mean_heading` and `polarization`.

    Examples:
        >>> from skcomplex.models.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 model, with the density flock.f.

Examples:

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

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

@register_model(
    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 model, with the density `flock.f`.

    Examples:
        >>> from skcomplex.models.physics.active_matter import bertin_droz_gregoire
        >>> bertin_droz_gregoire(2, 8).elements[0].check().prepare().size_
        32
    """
    axis = AxisSpec('theta', angles, 0.0, 2 * math.pi, periodic=True)
    space = PhaseSpace(LatticeSpace((positions, positions)), [axis]).check().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 model.

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 model, with the agent group birds and the observables mean_cos, mean_sin, and polarization.

Examples:

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

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

@register_model(
    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 model.

    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 model, with the agent group `birds` and the observables `mean_cos`, `mean_sin`, and `polarization`.

    Examples:
        >>> from skcomplex.models.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 model, with the positions swimmers.x and orientations swimmers.theta.

Examples:

>>> from skcomplex.models.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/models/physics/active_matter/_particles.py

@register_model(
    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 model, with the positions `swimmers.x` and orientations `swimmers.theta`.

    Examples:
        >>> from skcomplex.models.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 model, with the positions birds.x and velocities birds.x_d1.

Examples:

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

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

@register_model(
    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 model, with the positions `birds.x` and velocities `birds.x_d1`.

    Examples:
        >>> from skcomplex.models.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))

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 model, with the positions walkers.x, velocities walkers.v, and depots walkers.e.

Examples:

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

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

@register_model(
    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 model, with the positions `walkers.x`, velocities `walkers.v`, and depots `walkers.e`.

    Examples:
        >>> from skcomplex.models.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 model, with the fields nematic.rho and nematic.Q.

Examples:

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

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

@register_model(
    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 model, with the fields `nematic.rho` and `nematic.Q`.

    Examples:
        >>> from skcomplex.models.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 model, with the positions flock.x, the headings flock.theta, and the observables polarization and nematic_order.

Examples:

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

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

@register_model(
    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 model, with the positions `flock.x`, the headings `flock.theta`, and the observables `polarization`
            and `nematic_order`.

    Examples:
        >>> from skcomplex.models.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é model: 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 model, with the positions flock.x, the headings flock.theta, and the observables polarization and nematic_order.

Examples:

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

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

@register_model(
    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é model: 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 model, with the positions `flock.x`, the headings `flock.theta`, and the observables `polarization`
            and `nematic_order`.

    Examples:
        >>> from skcomplex.models.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)

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 model 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 model, with the cell density colony.rho and the chemical colony.c.

Examples:

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

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

@register_model(
    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 model 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 model, with the cell density `colony.rho` and the chemical `colony.c`.

    Examples:
        >>> from skcomplex.models.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 model, with the positions swimmers.x and orientations swimmers.theta.

Examples:

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

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

@register_model(
    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 model, with the positions `swimmers.x` and orientations `swimmers.theta`.

    Examples:
        >>> from skcomplex.models.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 model, with the positions swimmers.x and orientations swimmers.theta.

Examples:

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

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

@register_model(
    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 model, with the positions `swimmers.x` and orientations `swimmers.theta`.

    Examples:
        >>> from skcomplex.models.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 model, with the mean headings heading_a and heading_b of the two species.

Examples:

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

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

@register_model(
    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 model, with the mean headings `heading_a` and `heading_b` of the two species.

    Examples:
        >>> from skcomplex.models.physics.active_matter import nonreciprocal_vicsek
        >>> nonreciprocal_vicsek(20).elements[0].check().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 model, with the agent group bacteria.

Examples:

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

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

@register_model(
    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 model, with the agent group `bacteria`.

    Examples:
        >>> from skcomplex.models.physics.active_matter import run_and_tumble
        >>> run_and_tumble(10).elements[0].check().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'))

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 model, with the fields flock.rho, flock.vx, and flock.vy.

Examples:

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

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

@register_model(
    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 model, with the fields `flock.rho`, `flock.vx`, and `flock.vy`.

    Examples:
        >>> from skcomplex.models.physics.active_matter import toner_tu
        >>> toner_tu((8, 8)).elements[0].check().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)

Return the Vicsek model: 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.

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

Returns:

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

Examples:

>>> from skcomplex.models.physics.active_matter import vicsek
>>> vicsek(10).elements[0].size
10

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

@register_model(
    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,
) -> System:
    """Return the Vicsek model: 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.

    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 model, with the positions `flock.x`, the headings `flock.theta`, and the observables `polarization`
            and `nematic_order`.

    Examples:
        >>> from skcomplex.models.physics.active_matter import vicsek
        >>> vicsek(10).elements[0].size
        10
    """
    alignment = HeadingAlignmentInteraction('alignment', target='flock')
    return _flock(VicsekParticleElement, n_particles, size, seed, v0=v0, alignment=alignment, eta=eta)