Depot-pumped Brownian particles

Brownian particles with energy depots, which start at rest and take up the speed sqrt(q/gamma - c/d).

• Name: depot_brownian_particles
• Field: physics/active_matter
• Formalism: stochastic differential equations
• References: Schweitzer, Ebeling, and Tilch 1998; Helbing 2001

Factory

@register_model(
    field='physics/active_matter',
    tags=('active matter', 'sde', 'energy depot', 'self-propulsion'),
    references=('Schweitzer, Ebeling, and Tilch 1998', 'Helbing 2001'),
    defaults={'until': 10.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.001)},
)
def depot_brownian_particles(
    n_particles: int = 50,
    *,
    gamma: float = 20.0,
    q: float = 10.0,
    c: float = 1.0,
    d: float = 10.0,
    D: float = 0.01,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return Brownian particles with energy depots, which start at rest and take up the speed `sqrt(q/gamma - c/d)`.

    Args:
        n_particles:
            The number of particles.

        gamma:
            The friction coefficient.

        q:
            The rate of energy intake.

        c:
            The rate of internal energy loss.

        d:
            The rate of conversion of depot energy into motion.

        D:
            The intensity of the velocity noise.

        random_state:
            The seed or generator of the small initial velocities.

    Returns:
        system:
            The model, with the positions `walkers.x`, velocities `walkers.v`, and depots `walkers.e`.

    Examples:
        >>> from skcomplex.models.physics.active_matter import depot_brownian_particles
        >>> depot_brownian_particles(5).elements[0].values['q']
        10.0
    """
    velocities = np.random.default_rng(random_state).normal(0.0, 0.01, (n_particles, 2))
    group = GroupElement(
        'walkers',
        element=DepotParticleElement,
        size=n_particles,
        values={'v': velocities, 'gamma': gamma, 'q': q, 'c': c, 'd': d, 'D': D},
    )
    return System(elements=[group], space=ContinuousSpace(n_dim=2))

Components

class DepotParticleElement(Element):
    """Brownian particle pumped by an internal energy depot, as in the model of Schweitzer, Ebeling, and Tilch.

    The depot takes energy at the rate `q`, loses it at the rate `c e`, and converts it into motion at the rate
    `d e |v|^2`, so `de/dt = q - c e - d e |v|^2` and `dv = (-gamma v + d e v) dt + sqrt(2 D) dW`. Above
    `q d = c gamma` the particle moves at the speed `sqrt(q / gamma - c / d)`.

    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.

        v:
            The initial velocity.

        e:
            The initial energy of the depot.

        gamma:
            The friction coefficient.

        q:
            The rate of energy intake.

        c:
            The rate of internal energy loss.

        d:
            The rate of conversion of depot energy into motion.

        D:
            The intensity of the velocity noise.

    Examples:
        >>> from skcomplex.models.physics.active_matter import DepotParticleElement
        >>> DepotParticleElement('walker').q
        10.0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = POINTS,
        *,
        x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
        v: Annotated[npt.ArrayLike, StateSpec(vector=True, role='velocity')] = 0.0,
        e: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
        gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 20.0,
        q: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 10.0,
        c: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
        d: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 10.0,
        D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
    ) -> None:
        super().__init__(label, visual=visual)
        self.x = x
        self.v = v
        self.e = e
        self.gamma = gamma
        self.q = q
        self.c = c
        self.d = d
        self.D = D

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the motion, the friction and pumping of the velocity, its noise, and the depot balance.

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

        Returns:
            laws:
                The laws of the component.
        """
        return [
            RateLaw(s.x, s.v),
            RateLaw(s.v, (-s.gamma + s.d * s.e) * s.v),
            NoiseLaw(s.v, sqrt(2 * s.D)),
            RateLaw(s.e, s.q - s.c * s.e - s.d * s.e * s.v.dot_product(s.v)),
        ]

Equations

walkers.x[i, 0]: diff(x_walkers_0(t), t) == v_walkers_0(t)
walkers.x[i, 1]: diff(x_walkers_1(t), t) == v_walkers_1(t)
walkers.v[i, 0]: diff(v_walkers_0(t), t) == (d_walkers*e_walkers(t) - gamma_walkers)*v_walkers_0(t)
walkers.v[i, 1]: diff(v_walkers_1(t), t) == (d_walkers*e_walkers(t) - gamma_walkers)*v_walkers_1(t)
walkers.e[i]: diff(e_walkers(t), t) == -(v_walkers_0(t)^2 + v_walkers_1(t)^2)*d_walkers*e_walkers(t) - c_walkers*e_walkers(t) + q_walkers