Dissipative zero-range process on a tree
The dissipative zero-range process on a tree, with energy injected at the root and lost at the leaves.
- Name:
dissipative_zero_range - Field:
physics/statistical_mechanics - Formalism: jump process
- References: Bertin 2021; Bertin 2005
Factory
@register_model(
field='physics/statistical_mechanics',
tags=('zero-range process', 'energy cascade', 'nonequilibrium', 'tree', 'jump process'),
references=('Bertin 2021', 'Bertin 2005'),
defaults={'until': 200.0, 'dt': 1.0, 'backend': NumbaBackend()},
)
def dissipative_zero_range(
n_levels: int = 5,
*,
branches: int = 2,
nu: float = 1.0,
alpha: float = 0.5,
beta_ext: float = 1.0,
epsilon0: float = 1.0,
dissipation: float = 1.0,
) -> System:
"""Return the dissipative zero-range process on a tree, with energy injected at the root and lost at the leaves.
Level `j` of the tree has `m^(j - 1)` nodes, each with `m` branches to the next level. A quantum `epsilon0` of
energy crosses a branch between levels `j` and `j + 1` at rate `nu_j = nu k_j^alpha`, `k_j = m^(j - 1)`, in either
direction from an occupied node. A reservoir at the inverse temperature `beta_ext` exchanges quanta with the root,
adding one at rate `nu exp(-beta_ext epsilon0)` and removing one at rate `nu`, and every leaf loses a quantum at the
rate `dissipation`. The stationary distribution is a product of geometric laws `z_j^n`, with `z_j` the solution of
the linear balance of the currents `nu_j (z_j - z_(j + 1))` at every level, and the injected flux is
`nu (exp(-beta_ext epsilon0) - z_1)`.
Args:
n_levels:
The number of levels `M`.
branches:
The number `m` of branches of every node above the leaves.
nu:
The rate scale of transfers and of the exchange with the reservoir.
alpha:
The exponent of the growth of the transfer rate with the pseudo wave number `k_j`.
beta_ext:
The inverse temperature of the reservoir.
epsilon0:
The quantum of energy.
dissipation:
The rate at which every leaf loses a quantum.
Returns:
system:
The model, with the quanta `nodes.n` of the nodes in breadth-first order, from the root to the leaves.
Examples:
>>> from skcomplex.models.physics.statistical_mechanics import dissipative_zero_range
>>> system = dissipative_zero_range(3, branches=3)
>>> system.elements[0].size, len(system.interactions)
(13, 14)
"""
n_nodes = _level_start(n_levels + 1, branches)
nodes = GroupElement('nodes', element=TreeNodeElement, size=n_nodes, values={'x': _tree_layout(n_levels, branches)})
interactions: list[BaseInteraction] = []
for level in range(1, n_levels):
rate = nu * float(branches ** (level - 1)) ** alpha
start, below = _level_start(level, branches), _level_start(level + 1, branches)
for index in range(branches ** (level - 1)):
for branch in range(branches):
child = below + branches * index + branch
link = [f'nodes[{start + index}]', f'nodes[{child}]']
interactions.append(
ZeroRangeHopInteraction(f'transfer_{start + index}_{child}', link=link, p=rate, q=rate, b=0.0),
)
leaves = [f'nodes[{index}]' for index in range(_level_start(n_levels, branches), n_nodes)]
interactions += [
ReservoirInteraction('injection', sites=['nodes[0]'], inflow=nu * np.exp(-beta_ext * epsilon0), outflow=nu),
ReservoirInteraction('dissipation', sites=leaves, outflow=dissipation),
]
return System(elements=[nodes], interactions=interactions, space=ContinuousSpace(n_dim=2))
Components
class TreeNodeElement(ParticleSiteElement):
"""Node of a tree that holds a number `n` of energy quanta, drawn at a fixed position.
Args:
label:
The label of the node or the group.
visual:
How the element is drawn, or `None` for `default_visual_`.
n:
The initial number of quanta.
x:
The position of the node in the drawing of the tree.
Examples:
>>> from skcomplex.models.physics.statistical_mechanics import TreeNodeElement
>>> TreeNodeElement('root').n
0.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = VisualSpec(
'points',
{'position': 'x', 'radius': 0.04, 'color_by': 'n', 'color_range': (0.0, 2.0), 'states': (), 'trail': 0},
),
*,
n: float = 0.0,
x: Annotated[npt.ArrayLike, StateSpec(vector=True)] = 0.0,
) -> None:
super().__init__(label, visual=visual, n=n)
self.x = x
class ZeroRangeHopInteraction(TransferInteraction):
"""Hop of one particle along a link, forward at rate `p u(n)` and backward at rate `q u(n')`.
The rate depends only on the occupation of the departure site, through `u(n) = 1 + b / n`. It is a transfer of one
particle from the first site to the second, with jumps back at the reverse rate.
Args:
label:
The label of the link.
link:
The departure and the arrival site of a forward hop.
p:
The forward rate scale.
q:
The backward rate scale.
b:
The strength of the decrease of `u(n)` with the occupation, 0 for a constant rate.
Examples:
>>> from skcomplex.models.physics.statistical_mechanics import ZeroRangeHopInteraction
>>> ZeroRangeHopInteraction('hop', link=['sites[0]', 'sites[1]']).b
4.0
"""
def __init__(
self,
label: str,
*,
link: Annotated[list[str], RoleSpec(variadic=True)],
p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
q: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
b: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
) -> None:
super().__init__(
label,
source='',
target='',
state='n',
amount=1.0,
rate=self._get_forward_rate,
reverse_rate=self._get_backward_rate,
)
self.link = link
self.p = p
self.q = q
self.b = b
def _get_holders(self, namespace: SimpleNamespace) -> list[tuple[SimpleNamespace, SimpleNamespace]]:
"""Return the two sites of the link, the departure and the arrival of a forward hop."""
first, second = namespace.link
return [(first, second)]
@staticmethod
def _get_forward_rate(s: SimpleNamespace) -> Expression:
"""Return the rate `p u(n)` of a hop from the first site."""
return s.p * ZeroRangeHopInteraction._departure_rate(s.source.n, s.b)
@staticmethod
def _get_backward_rate(s: SimpleNamespace) -> Expression:
"""Return the rate `q u(n')` of a hop from the second site."""
return s.q * ZeroRangeHopInteraction._departure_rate(s.target.n, s.b)
@staticmethod
def _departure_rate(n: Expression, b: Expression) -> Expression:
"""Return the zero-range rate `u(n) = 1 + b / n` of a site with `n` particles, or 0 for an empty site."""
return where(n > 0, 1 + b / where(n > 0, n, 1), 0)
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)