Motility-induced phase separation
Repulsive active Brownian disks of unit diameter, which phase separate into clusters at high speed.
- Name:
motility_induced_phase_separation - Field:
physics/active_matter - Formalism: stochastic differential equations
- References: Cates and Tailleur 2015; Bechinger et al. 2016; Fily and Marchetti 2012
Factory
@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)
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)])
class SoftSphereInteraction(PairCouplingInteraction):
"""Soft repulsion between overlapping spheres, `k * (cutoff - r) ** 2 / 2` closer than the cutoff.
Args:
label:
The label of the interaction.
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 contact distance, the sum of the two radii.
position:
The name of the vector state that holds the positions of the members.
write:
The name of the state whose rate the force changes, the position of order 2 by default.
divisor:
The name of the parameter that holds the masses of the members, or `None` for unit masses.
period:
The side of a periodic box, for minimum image distances, or `None` for the period of the space.
k:
The stiffness.
Examples:
>>> from skcomplex.models.physics.common.particles import SoftSphereInteraction
>>> SoftSphereInteraction('contact', first='gas', second='gas', cutoff=0.1, k=50.0).k
50.0
"""
def __init__(
self,
label: str,
first: str,
second: str,
cutoff: float = 1.0,
position: str = 'x',
*,
write: str = 'x',
divisor: str | None = None,
period: float | None = None,
k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
) -> None:
super().__init__(
label,
first,
second,
cutoff=cutoff,
position=position,
write=write,
potential=self._get_soft_sphere_energy,
divisor=divisor,
period=period,
)
self.k = k
@staticmethod
def _get_soft_sphere_energy(s: SimpleNamespace, r: Expression) -> Expression:
"""Return `k * (cutoff - r) ** 2 / 2`."""
return s.k * (s.cutoff - r) ** 2 / 2
Equations
swimmers.theta[i]: diff(theta_swimmers(t), t) == omega_swimmers