Asymmetric exclusion process with open boundaries

The asymmetric simple exclusion process on an open line, in continuous time.

• Name: asep
• Field: physics/statistical_mechanics
• Formalism: jump process
• References: Derrida, Evans, Hakim, and Pasquier 1993; Bertin 2021

Factory

@register_model(
    field='physics/statistical_mechanics',
    tags=('exclusion', 'nonequilibrium', 'phase transition', 'traffic', 'jump process'),
    references=('Derrida, Evans, Hakim, and Pasquier 1993', 'Bertin 2021'),
    defaults={'until': 500.0, 'dt': 2.0, 'backend': NumbaBackend()},
)
def asep(
    n_sites: int = 100,
    *,
    alpha: float = 0.3,
    gamma: float = 0.7,
    p: float = 1.0,
    q: float = 0.0,
    density: float = 0.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the asymmetric simple exclusion process on an open line, in continuous time.

    A particle hops to the empty right neighbour at rate `p` and to the empty left neighbour at rate `q`, enters the
    empty first site at rate `alpha`, and leaves the last site at rate `gamma`, the exit rate that Bertin writes
    `gamma` and others write `beta`. For `p = 1` and `q = 0` the current is `alpha (1 - alpha)` in the low-density
    phase `alpha < 1/2, gamma > alpha`, `gamma (1 - gamma)` in the high-density phase `gamma < 1/2, gamma < alpha`, and
    `1/4` in the maximal-current phase `alpha, gamma > 1/2`.

    Args:
        n_sites:
            The number of sites.

        alpha:
            The entry rate.

        gamma:
            The exit rate.

        p:
            The rate of hops to the right.

        q:
            The rate of hops to the left.

        density:
            The initial fraction of occupied sites.

        random_state:
            The seed or generator of the initial occupations.

    Returns:
        system:
            The model, with the occupations `lane.n`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import asep
        >>> len(asep(10).interactions)
        3
    """
    occupied = (np.random.default_rng(random_state).random(n_sites) < density).astype(float)
    lane = LatticeElement(
        'lane',
        element=ParticleSiteElement,
        space=LatticeSpace(n_sites, boundary='fixed'),
        values={'n': occupied},
    )
    interactions = [
        ExclusionHopsInteraction('hops', sites='lane', p=p, q=q),
        ReservoirInteraction('entry', sites=['lane[0]'], inflow=alpha, exclusion=True),
        ReservoirInteraction('exit', sites=[f'lane[{n_sites - 1}]'], outflow=gamma, exclusion=True),
    ]
    return System(elements=[lane], interactions=interactions)

Components

class ParticleSiteElement(Element):
    """Site of a lattice of particles, holding the number `n` of particles or quanta on it.

    Args:
        label:
            The label of the site or the group.

        visual:
            How the element is drawn, or `None` for `default_visual_`.

        n:
            The initial number of particles.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import ParticleSiteElement
        >>> ParticleSiteElement('site', n=3.0).n
        3.0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        n: Annotated[float, StateSpec(domain=NN)] = 0.0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.n = n

class ExclusionHopsInteraction(TransferInteraction):
    """Hops of a particle from every site of a line to an empty neighbour, to the right at rate `p` and left at `q`.

    No particle hops past the ends of the line. It is a transfer of one particle along every edge of the line, to the
    right and back, from an occupied site to an empty one.

    Args:
        label:
            The label of the hops.

        sites:
            The group of sites, in line order.

        p:
            The rate of hops to the right.

        q:
            The rate of hops to the left.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import ExclusionHopsInteraction
        >>> ExclusionHopsInteraction('hops', sites='lane').p
        1.0
    """

    def __init__(
        self,
        label: str,
        *,
        sites: Annotated[str, RoleSpec()],
        p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
        q: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    ) -> None:
        super().__init__(
            label,
            source=sites,
            target=sites,
            state='n',
            amount=1.0,
            rate=self._get_right_rate,
            reverse_rate=self._get_left_rate,
            offset=1,
            boundary='fixed',
        )
        self.sites = sites
        self.p = p
        self.q = q

    def _get_holders(self, namespace: SimpleNamespace) -> list[tuple[SimpleNamespace, SimpleNamespace]]:
        """Return the sites, whose every member hops to the next one and back."""
        return [(namespace.sites, namespace.sites)]

    @staticmethod
    def _get_right_rate(s: SimpleNamespace) -> Expression:
        """Return the rate of a hop to the right, `p` from an occupied site to an empty one."""
        return s.p * s.source.n * (1 - s.target.n)

    @staticmethod
    def _get_left_rate(s: SimpleNamespace) -> Expression:
        """Return the rate of a hop to the left, `q` from an occupied site to an empty one."""
        return s.q * s.target.n * (1 - s.source.n)

class ReservoirInteraction(TransferInteraction):
    """Reservoir that adds one particle to each of its sites at rate `inflow` and removes one at rate `outflow`.

    A particle leaves only an occupied site. With exclusion a particle enters only an empty site. It is a transfer of
    one particle from the environment into every site, with jumps back to the environment at the outflow rate.

    Args:
        label:
            The label of the reservoir.

        sites:
            The sites in contact with the reservoir.

        inflow:
            The rate at which a particle enters a site.

        outflow:
            The rate at which a particle leaves an occupied site.

        exclusion:
            Whether a site holds at most one particle.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import ReservoirInteraction
        >>> ReservoirInteraction('bath', sites=['lane[0]'], inflow=0.5).outflow
        0.0
    """

    def __init__(
        self,
        label: str,
        *,
        sites: Annotated[list[str], RoleSpec(variadic=True)],
        inflow: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
        outflow: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
        exclusion: bool = False,
    ) -> None:
        super().__init__(
            label,
            source='',
            target='',
            state=Source(),
            into='n',
            amount=1.0,
            rate=self._get_entry_rate,
            reverse_rate=self._get_exit_rate,
        )
        self.sites = sites
        self.inflow = inflow
        self.outflow = outflow
        self.exclusion = exclusion

    def _get_holders(self, namespace: SimpleNamespace) -> list[tuple[SimpleNamespace, SimpleNamespace]]:
        """Return every site as the target of its own entries."""
        return [(site, site) for site in namespace.sites]

    def _get_entry_rate(self, s: SimpleNamespace) -> Expression:
        """Return the rate at which a particle enters a site, only an empty one under exclusion."""
        return s.inflow * (1 - s.target.n if self.exclusion else 1)

    @staticmethod
    def _get_exit_rate(s: SimpleNamespace) -> Expression:
        """Return the rate at which a particle leaves an occupied site."""
        return s.outflow * where(s.target.n > 0, 1, 0)