Active Ornstein-Uhlenbeck particles

Active Ornstein-Uhlenbeck particles in the plane, driven by exponentially correlated propulsion.

• Name: active_ornstein_uhlenbeck_particles
• Field: physics/active_matter
• Formalism: stochastic differential equations
• References: Bechinger et al. 2016; Fodor et al. 2016

Factory

@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))

Components

class OrnsteinUhlenbeckParticleElement(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.

    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.

        u:
            The initial propulsion velocity.

        tau:
            The persistence time.

        D:
            The active diffusion coefficient.

        D_t:
            The translational diffusion coefficient.

    Examples:
        >>> from skcomplex.models.physics.active_matter import OrnsteinUhlenbeckParticleElement
        >>> OrnsteinUhlenbeckParticleElement('swimmer').tau
        1.0
    """

    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

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the propulsion, its relaxation and noise, and the translational noise.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the component.
        """
        return [
            RateLaw(s.x, s.u),
            NoiseLaw(s.x, sqrt(2 * s.D_t)),
            RateLaw(s.u, -s.u / s.tau),
            NoiseLaw(s.u, sqrt(2 * s.D) / s.tau),
        ]

Equations

swimmers.x[i, 0]: diff(x_swimmers_0(t), t) == u_swimmers_0(t)
swimmers.x[i, 1]: diff(x_swimmers_1(t), t) == u_swimmers_1(t)
swimmers.u[i, 0]: diff(u_swimmers_0(t), t) == -u_swimmers_0(t)/tau_swimmers
swimmers.u[i, 1]: diff(u_swimmers_1(t), t) == -u_swimmers_1(t)/tau_swimmers