Chiral active particles
Chiral active particles, which swim along circles of radius v / omega that drift by diffusion.
- Name:
chiral_active_particles - Field:
physics/active_matter - Formalism: stochastic differential equations
- References: Bechinger et al. 2016; van Teeffelen and Löwen 2008
Factory
@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))
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 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(t))
swimmers.x[i, 1]: diff(x_swimmers_1(t), t) == speed_propulsion*sin(theta_swimmers(t))
swimmers.theta[i]: diff(theta_swimmers(t), t) == omega_swimmers