Nonreciprocal Vicsek model
Two species of Vicsek movers in a periodic box whose alignment with each other is nonreciprocal.
- Name:
nonreciprocal_vicsek - Field:
physics/active_matter - Formalism: discrete map
- References: Fruchart, Hanai, Littlewood, and Vitelli 2021; te Vrugt and Wittkowski 2024
Factory
@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(),
)
Components
class NonreciprocalMoverElement(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.
Args:
label:
The label of the mover.
visual:
How the element is drawn, or `None` for `default_visual_`.
x:
The position.
theta:
The heading.
species:
The species, 0 for A or 1 for B.
eta:
The width of the uniform noise of the heading.
J_aa:
The alignment of A with A.
J_ab:
The alignment of A with B.
J_ba:
The alignment of B with A.
J_bb:
The alignment of B with B.
radius:
The radius of the neighbourhood.
Examples:
>>> from skcomplex.models.physics.active_matter import NonreciprocalMoverElement
>>> NonreciprocalMoverElement('mover', J_ba=-1.0).J_ba
-1.0
"""
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
def _coupling(self, namespace: SimpleNamespace) -> Expression:
"""Return the alignment of a member with a neighbour, by the species of both."""
s = namespace
first, second = s.species, s.other.species
return (
s.J_aa * (1 - first) * (1 - second)
+ s.J_ab * (1 - first) * second
+ s.J_ba * first * (1 - second)
+ s.J_bb * first * second
)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the alignment of the heading with the species-weighted neighbours, plus noise.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
The update of the heading.
"""
s = namespace
over = RadiusNeighborhood(self.radius, center=True)
direction = neighbors(s.other.theta, over=over, reduce='circular-mean', weight=self._coupling(s))
return [UpdateLaw(s.theta, direction + uniform('noise', -s.eta / 2, s.eta / 2))]
class PropulsionStepInteraction(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.
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 distance moved in one step.
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 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)]
"""
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
def _step(self, namespace: SimpleNamespace) -> FreeModuleElement:
"""Return the step `speed (cos theta, sin theta)`."""
heading = getattr(namespace.target, self.heading)
return namespace.speed * vector(SR, [cos(heading), sin(heading)])
Equations
movers.x[i, 0]: x_movers_0(t + 1) == speed_propulsion*cos(theta_movers(t)) + x_movers_0(t)
movers.x[i, 1]: x_movers_1(t + 1) == speed_propulsion*sin(theta_movers(t)) + x_movers_1(t)
movers.theta[i]: theta_movers(t + 1) == eta_movers*uniform_noise_movers - 1/2*eta_movers + arctan2(_neighbours_1, _neighbours_0)
movers.species[i]: species_movers(t + 1) == species_movers(t)