Active Brownian particles with local alignment
Active Brownian particles in a periodic box whose headings align with the neighbours within d0.
- Name:
aligning_active_brownian_particles - Field:
physics/active_matter - Formalism: stochastic differential equations
- References: Bertin 2021; Peruani, Deutsch, and Bär 2008; Chaté 2020
Factory
@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'),
)
Components
class ActiveParticleElement(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.
Args:
label:
The label of the particle or the group.
visual:
How the element is drawn, or `None` for `default_visual_`.
x:
The initial position.
theta:
The initial orientation.
D_t:
The translational diffusion coefficient.
D_r:
The rotational diffusion coefficient.
omega:
The angular velocity of a chiral particle.
Examples:
>>> from skcomplex.models.physics.active_matter import ActiveParticleElement
>>> ActiveParticleElement('swimmer', D_r=2.0).D_r
2.0
"""
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
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the translational noise, and the turning and the rotational noise.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
return [
NoiseLaw(s.x, sqrt(2 * s.D_t)),
RateLaw(s.theta, s.omega),
NoiseLaw(s.theta, sqrt(2 * s.D_r)),
]
class PolarAlignmentInteraction(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.
Args:
label:
The label of the interaction, unique among the parts of its system.
first:
The path of the group or particle of the first members.
second:
The path of the group or particle of the second members.
cutoff:
The interaction range `d0`.
position:
The name of the vector state that holds the positions of the members.
gamma:
The strength of the alignment.
period:
The side of a periodic box, for minimum image distances, or `None` for the period of the space.
Examples:
>>> from skcomplex.models.physics.active_matter import PolarAlignmentInteraction
>>> PolarAlignmentInteraction('align', first='swimmers', second='swimmers', gamma=0.5).gamma
0.5
"""
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
@staticmethod
def _torque(namespace: SimpleNamespace) -> Expression:
"""Return the torque `gamma sin(theta_j - theta_i)` on the first member of a pair."""
return namespace.gamma * sin(namespace.second.theta - namespace.first.theta)
class SelfPropulsionInteraction(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.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the particle or the group whose positions move.
speed:
The self-propulsion speed.
state:
The name of the vector state that holds the position.
heading:
The name of the state that holds the heading.
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))]
"""
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
def _get_rate(self, namespace: SimpleNamespace) -> FreeModuleElement:
"""Return the velocity `speed (cos theta, sin theta)`.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The rate of the position.
"""
heading = getattr(namespace.target, self.heading)
return namespace.speed * vector(SR, [cos(heading), sin(heading)])
Equations
swimmers.x[i, 0]: diff(x_swimmers_0(t), t) == speed_propulsion*cos(theta_swimmers)
swimmers.x[i, 1]: diff(x_swimmers_1(t), t) == speed_propulsion*sin(theta_swimmers)