Soft matter

Soft matter: Brownian particles, optical and hydrodynamically coupled traps, anomalous diffusion, and transitions.

AnisotropicParticleElement(label, visual=None, *, x=0.0, y=0.0, z=0.0, ux=0.0, uy=0.0, uz=1.0, psi=0.0, D_par=1.0, D_perp=1.0, D_r=0.75, D_spin=0.75, D_c=0.0)

Bases: Element

Axisymmetric Brownian particle in three dimensions, whose diffusion tensor turns with its axis.

The position r = (x, y, z) diffuses with D_par along the unit axis u = (ux, uy, uz) and with D_perp across it, dr = (sqrt(2 D_par) u u + sqrt(2 D_perp) (I - u u)) dW. The axis diffuses on the unit sphere, du = -2 D_r u dt + sqrt(2 D_r) (I - u u) dW_r in the Itô form, and a boundary rescales it to unit length after every step against the drift of the discrete steps. The spin psi about the axis diffuses with D_spin, and its noise is correlated with the displacement along the axis through the coupling D_c of a chiral body, so that <dr_par dpsi> = 2 D_c dt, which needs D_c² <= D_par D_spin.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the particle is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial first coordinate.	0.0
y	Annotated[float, StateSpec()]	The initial second coordinate.	0.0
z	Annotated[float, StateSpec()]	The initial third coordinate.	0.0
ux	Annotated[float, StateSpec()]	The initial first component of the axis.	0.0
uy	Annotated[float, StateSpec()]	The initial second component of the axis.	0.0
uz	Annotated[float, StateSpec()]	The initial third component of the axis.	1.0
psi	Annotated[float, StateSpec()]	The initial spin angle about the axis.	0.0
D_par	Annotated[float, ParameterSpec(domain=open(0, oo))]	The translational diffusion coefficient along the axis.	1.0
D_perp	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The translational diffusion coefficient across the axis.	1.0
D_r	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rotational diffusion coefficient of the axis.	0.75
D_spin	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rotational diffusion coefficient about the axis.	0.75
D_c	Annotated[float, ParameterSpec()]	The coupling of the translation along the axis and the rotation about it.	0.0

Attributes:

Name	Type	Description
default_visual_		The visual of the particle when it sets none.

Examples:

>>> from skcomplex.systems.physics.soft_matter import AnisotropicParticleElement
>>> AnisotropicParticleElement('janus').uz
1.0

Source code in src/skcomplex/systems/physics/soft_matter/_rigid_bodies.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    y: Annotated[float, StateSpec()] = 0.0,
    z: Annotated[float, StateSpec()] = 0.0,
    ux: Annotated[float, StateSpec()] = 0.0,
    uy: Annotated[float, StateSpec()] = 0.0,
    uz: Annotated[float, StateSpec()] = 1.0,
    psi: Annotated[float, StateSpec()] = 0.0,
    D_par: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    D_perp: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    D_r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.75,
    D_spin: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.75,
    D_c: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.z = z
    self.ux = ux
    self.uy = uy
    self.uz = uz
    self.psi = psi
    self.D_par = D_par
    self.D_perp = D_perp
    self.D_r = D_r
    self.D_spin = D_spin
    self.D_c = D_c

BlinkingTrapsInteraction(label, *, target, k=1.0, on_time=1.0, off_time=0.0, tensor='rotne-prager-yamakawa', root='cholesky', radius=1.0, viscosity=1.0, temperature=1.0, dt=0.01)

Bases: MobilityInteraction

Blinking optical traps on spheres that move through a shared fluid, in steps of Brownian dynamics.

Every sphere feels the harmonic force -k (x - center) of its own trap while the traps are on, for on_time of every period on_time + off_time, and no force while they are off. The mobility turns the forces on all spheres into their velocities, and the thermal energy temperature adds Brownian displacements whose correlations between spheres follow the same mobility, so spheres released from their traps diffuse together (Batchelor 1976). Every step lasts dt, and the time of the system counts steps.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the group of spheres.	required
k	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The stiffness of every trap.	1.0
on_time	float	The time the traps stay on in every period.	1.0
off_time	float	The time the traps stay off in every period.	0.0
tensor	str	'free-draining', 'oseen', or 'rotne-prager-yamakawa'.	'rotne-prager-yamakawa'
root	str	'cholesky' or 'lanczos'.	'cholesky'
radius	Annotated[float, ParameterSpec(domain=open(0, oo))]	The radius of the spheres.	1.0
viscosity	Annotated[float, ParameterSpec(domain=open(0, oo))]	The viscosity of the fluid.	1.0
temperature	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The thermal energy kT.	1.0
dt	Annotated[float, ParameterSpec(domain=open(0, oo))]	The length of a step.	0.01

Examples:

>>> from skcomplex.systems.physics.soft_matter import BlinkingTrapsInteraction
>>> BlinkingTrapsInteraction('tweezers', target='beads', off_time=0.2).form
'map'

Source code in src/skcomplex/systems/physics/soft_matter/_coupled_traps.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    on_time: float = 1.0,
    off_time: float = 0.0,
    tensor: str = 'rotne-prager-yamakawa',
    root: str = 'cholesky',
    radius: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    viscosity: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    temperature: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    dt: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
) -> None:
    super().__init__(
        label,
        target=target,
        tensor=tensor,
        form='map',
        force=self._get_trap_force,
        root=root,
        radius=radius,
        viscosity=viscosity,
        temperature=temperature,
        dt=dt,
    )
    self.k = k
    self.on_time = on_time
    self.off_time = off_time

BrownianParticleElement(label, visual=POINTS, *, x=0.0, gamma=1.0, temperature=1.0)

Bases: Element

Overdamped Brownian particle at the position x, with the friction gamma and the thermal energy temperature.

The particle moves by dx = F / gamma dt + sqrt(2 temperature / gamma) dW, with the thermal noise of ThermalNoiseInteraction and the forces of other interactions, which divide by gamma.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The friction coefficient.	1.0
temperature	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The thermal energy.	1.0

Examples:

>>> from skcomplex.systems.physics.common.particles import BrownianParticleElement
>>> BrownianParticleElement('bead', temperature=2.0).temperature
2.0

Source code in src/skcomplex/systems/physics/common/particles/_brownian_particle.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    temperature: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.gamma = gamma
    self.temperature = temperature

ColoredNoiseParticleElement(label, visual=None, *, v=0.0, xi=0.0, gamma=1.0, D=1.0, tau=1.0)

Bases: Element

Particle with the velocity v, driven by the noise xi of correlation time tau.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle.	required
visual	VisualSpec | None	How the particle is drawn, or None for none.	None
v	Annotated[float, StateSpec()]	The initial velocity.	0.0
xi	Annotated[float, StateSpec()]	The initial noise.	0.0
gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The friction coefficient.	1.0
D	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The noise intensity.	1.0
tau	Annotated[float, ParameterSpec(domain=open(0, oo))]	The correlation time of the noise.	1.0

Examples:

>>> from skcomplex.systems.physics.soft_matter import ColoredNoiseParticleElement
>>> ColoredNoiseParticleElement('particle', tau=0.5).tau
0.5

Source code in src/skcomplex/systems/physics/soft_matter/_colored_noise.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    v: Annotated[float, StateSpec()] = 0.0,
    xi: Annotated[float, StateSpec()] = 0.0,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.v = v
    self.xi = xi
    self.gamma = gamma
    self.D = D
    self.tau = tau

ContinuousTimeWalkerElement(label, visual=POINTS, *, x=0.0, alpha=0.5, tau0=1.0, sigma=1.0)

Bases: Element

Walker that waits a power-law time tau0 r^(-1 / alpha) for r uniform, then jumps by a Gaussian step.

The jumps are a jump interaction of the system, which resets the position after every waiting time.

Parameters:

Name	Type	Description	Default
label	str	The label of the walker or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
alpha	Annotated[float, ParameterSpec(domain=open(0, oo))]	The tail exponent of the waiting times, below 1 for waiting times without a mean.	0.5
tau0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The shortest waiting time.	1.0
sigma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The standard deviation of the jumps along every axis.	1.0

Examples:

>>> from skcomplex.systems.physics.soft_matter import ContinuousTimeWalkerElement
>>> ContinuousTimeWalkerElement('walker', alpha=0.7).alpha
0.7

Source code in src/skcomplex/systems/physics/soft_matter/_anomalous.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    tau0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.alpha = alpha
    self.tau0 = tau0
    self.sigma = sigma

DimerElement(label, visual=None, *, c=0.0, phi=0.0, theta1=0.0, theta2=0.0, L=1.0, v0=0.0, D_t=0.1, D_r=1.0)

Bases: Element

Two Brownian or active spheres joined by a rigid rod, described by the centre and the angle of the rod.

Each sphere diffuses with D_t and swims at the speed v0 along its own orientation, theta1 or theta2, which diffuse with D_r, and the rod keeps them at the distance L. Moving both spheres freely and then restoring their distance moves the centre c by the mean of their displacements and turns the rod by the part of their relative displacement across it, so in the limit of small steps dc = v0 (n1 + n2) / 2 dt + sqrt(D_t) dW and dphi = v0 (n1 - n2) . e / L dt + 2 sqrt(D_t) / L dW, with e = (-sin phi, cos phi) across the rod.

Parameters:

Name	Type	Description	Default
label	str	The label of the dimer or the group.	required
visual	VisualSpec | None	How the dimer is drawn, or None for default_visual_.	None
c	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial centre.	0.0
phi	Annotated[float, StateSpec(role=direction)]	The initial angle of the rod.	0.0
theta1	Annotated[float, StateSpec(period=2 * pi)]	The initial orientation of the first sphere.	0.0
theta2	Annotated[float, StateSpec(period=2 * pi)]	The initial orientation of the second sphere.	0.0
L	Annotated[float, ParameterSpec(domain=open(0, oo))]	The length of the rod.	1.0
v0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The swimming speed of each sphere, 0 for Brownian spheres.	0.0
D_t	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The translational diffusion coefficient of each sphere.	0.1
D_r	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rotational diffusion coefficient of the orientation of each sphere.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the dimer when it sets none.

Examples:

>>> from skcomplex.systems.physics.soft_matter import DimerElement
>>> DimerElement('molecule').L
1.0

Source code in src/skcomplex/systems/physics/soft_matter/_rigid_bodies.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    c: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    phi: Annotated[float, StateSpec(role='direction')] = 0.0,
    theta1: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
    theta2: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
    L: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    v0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    D_t: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    D_r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.c = c
    self.phi = phi
    self.theta1 = theta1
    self.theta2 = theta2
    self.L = L
    self.v0 = v0
    self.D_t = D_t
    self.D_r = D_r

DoubleWellElement(label, visual=None, *, x=-1.0, amplitude=0.1, frequency=0.01, D=0.12)

Bases: Element

Particle at the position x in the double well x^4 / 4 - x^2 / 2, forced periodically and by noise.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle.	required
visual	VisualSpec | None	How the particle is drawn, or None for none.	None
x	Annotated[float, StateSpec(role=position)]	The initial position.	-1.0
amplitude	Annotated[float, ParameterSpec()]	The amplitude of the periodic forcing.	0.1
frequency	Annotated[float, ParameterSpec()]	The frequency of the forcing.	0.01
D	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The noise intensity.	0.12

Examples:

>>> from skcomplex.systems.physics.common.particles import DoubleWellElement
>>> DoubleWellElement('particle', D=0.2).D
0.2

Source code in src/skcomplex/systems/physics/common/particles/_double_well.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(role='position')] = -1.0,
    amplitude: Annotated[float, ParameterSpec()] = 0.1,
    frequency: Annotated[float, ParameterSpec()] = 0.01,
    D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.12,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.amplitude = amplitude
    self.frequency = frequency
    self.D = D

FlashingSawtoothInteraction(label, *, target, barrier=1.0, L=1.0, asymmetry=0.1, switches=())

Bases: RateInteraction

Sawtooth potential on overdamped particles on a line, turned on and off at given times.

Its force is the rate of the position, divided by the friction gamma of the target.

The potential of period L has a minimum at every multiple of L, rises linearly by barrier over the fraction asymmetry of the period, and falls back over the rest. It is on from time 0 to the first of the switches, off until the next one, and so on, and it keeps its state after the last one.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the particle or the group.	required
barrier	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The height of the sawtooth.	1.0
L	Annotated[float, ParameterSpec(domain=open(0, oo))]	The period of the sawtooth.	1.0
asymmetry	Annotated[float, ParameterSpec(domain=open(0, oo))]	The fraction of the period over which the potential rises, 0.5 for a symmetric potential.	0.1
switches	tuple[float, ...]	The increasing times at which the potential turns off and on in turn.	()

Examples:

>>> from skcomplex.systems.physics.soft_matter import FlashingSawtoothInteraction
>>> FlashingSawtoothInteraction('ratchet', target='beads', switches=(1.0, 2.5)).switches
(1.0, 2.5)

Source code in src/skcomplex/systems/physics/soft_matter/_ratchets.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    barrier: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    L: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    asymmetry: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
    switches: tuple[float, ...] = (),
) -> None:
    super().__init__(label, target=target, state='x', divisor='gamma')
    self.barrier = barrier
    self.L = L
    self.asymmetry = asymmetry
    self.switches = switches

FractionalWalkerElement(label, visual=POINTS, *, x=0.0, D=1.0, H=0.3)

Bases: Element

Walker driven by fractional Brownian motion, dx = sqrt(2 D) dB_H, whose increments are correlated in time.

Parameters:

Name	Type	Description	Default
label	str	The label of the walker or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
D	Annotated[float, ParameterSpec(domain=open(0, oo))]	The diffusion coefficient.	1.0
H	Annotated[float, ParameterSpec(domain=open(0, 1))]	The Hurst index of the driving fractional Brownian motion.	0.3

Examples:

>>> from skcomplex.systems.physics.soft_matter import FractionalWalkerElement
>>> FractionalWalkerElement('walker', H=0.3).H
0.3

Source code in src/skcomplex/systems/physics/soft_matter/_anomalous.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    D: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    H: Annotated[float, ParameterSpec(domain=RealSet.open(0, 1))] = 0.3,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.D = D
    self.H = H

FrictionalBlockElement(label, visual=None, *, x=0.0, v=0.0, m=1.0, k=1.0, anchor=0.0, mu_s=1.0, mu_d=0.8, amplitude=2.0, shake_time=2.0, relax_time=8.0, v_stick=0.02)

Bases: Element

Block of a chain on a rough floor, at the position x with the velocity v, joined to its neighbours.

The springs of stiffness k and an anchor spring pull the block, Coulomb friction holds it while the force stays below mu_s and slows it by mu_d while it slides, and periodic shaking of amplitude amplitude kicks it during the first shake_time of every shake_time + relax_time.

Parameters:

Name	Type	Description	Default
label	str	The label of the block.	required
visual	VisualSpec | None	How the block is drawn, or None for none.	None
x	Annotated[float, StateSpec()]	The initial position.	0.0
v	Annotated[float, StateSpec()]	The initial velocity.	0.0
m	Annotated[float, ParameterSpec(domain=open(0, oo))]	The mass.	1.0
k	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The stiffness of the springs between neighbours.	1.0
anchor	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The stiffness of the spring to the rest position.	0.0
mu_s	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The static friction force.	1.0
mu_d	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The dynamic friction force.	0.8
amplitude	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the shaking force.	2.0
shake_time	Annotated[float, ParameterSpec(domain=open(0, oo))]	The duration of every shaking.	2.0
relax_time	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time between shakings.	8.0
v_stick	Annotated[float, ParameterSpec(domain=open(0, oo))]	The speed below which the block can stick.	0.02

Examples:

>>> from skcomplex.systems.physics.soft_matter import FrictionalBlockElement
>>> FrictionalBlockElement('block', mu_s=2.0).mu_s
2.0

Source code in src/skcomplex/systems/physics/soft_matter/_friction.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    v: Annotated[float, StateSpec()] = 0.0,
    m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    anchor: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    mu_s: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    mu_d: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.8,
    amplitude: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
    shake_time: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
    relax_time: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 8.0,
    v_stick: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.02,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.v = v
    self.m = m
    self.k = k
    self.anchor = anchor
    self.mu_s = mu_s
    self.mu_d = mu_d
    self.amplitude = amplitude
    self.shake_time = shake_time
    self.relax_time = relax_time
    self.v_stick = v_stick

LevyFlierElement(label, visual=POINTS, *, x=0.0, mu=1.2, l0=1.0)

Bases: Element

Walker that jumps at every step by a length of power-law tail l^(-1 - mu) in a random direction.

Parameters:

Name	Type	Description	Default
label	str	The label of the walker or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
mu	Annotated[float, ParameterSpec(domain=open(0, oo))]	The tail exponent of the jump lengths, below 2 for jumps of infinite variance.	1.2
l0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The shortest jump length.	1.0

Examples:

>>> from skcomplex.systems.physics.soft_matter import LevyFlierElement
>>> LevyFlierElement('walker', mu=1.2).mu
1.2

Source code in src/skcomplex/systems/physics/soft_matter/_anomalous.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    mu: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.2,
    l0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.mu = mu
    self.l0 = l0

LevyWalkerElement(label, visual=POINTS, *, x=0.0, v=0.0, sigma=1.5, tau0=1.0, speed=1.0)

Bases: Element

Walker that flies at constant speed in a random direction for a power-law flight time, then turns.

The flight is a rate interaction of the position, and the turns a jump interaction that resets the velocity after every flight time.

Parameters:

Name	Type	Description	Default
label	str	The label of the walker or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	POINTS
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
v	Annotated[ArrayLike, StateSpec(vector=True, role=velocity)]	The initial velocity.	0.0
sigma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The tail exponent of the flight times, between 1 and 2 for superdiffusion.	1.5
tau0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The shortest flight time.	1.0
speed	Annotated[float, ParameterSpec(domain=open(0, oo))]	The speed of every flight.	1.0

Examples:

>>> from skcomplex.systems.physics.soft_matter import LevyWalkerElement
>>> LevyWalkerElement('walker', speed=2.0).speed
2.0

Source code in src/skcomplex/systems/physics/soft_matter/_anomalous.py

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,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.5,
    tau0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    speed: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.v = v
    self.sigma = sigma
    self.tau0 = tau0
    self.speed = speed

ModulatedParticleElement(label, visual=None, *, x=50.0, s0=1.0, ds=1.8, L=100.0)

Bases: Element

Particle at the position x in a box [0, L] whose walls reflect it.

The noise amplitude sigma = s0 + ds (x / L - 1/2) grows across the box. The walls reflect the particle at every frame, so frames must be short against the time to cross the box.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(role=position)]	The initial position.	50.0
s0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The noise amplitude at the middle of the box.	1.0
ds	Annotated[float, ParameterSpec()]	The change of the noise amplitude across the box.	1.8
L	Annotated[float, ParameterSpec(domain=open(0, oo))]	The length of the box.	100.0

Examples:

>>> from skcomplex.systems.physics.soft_matter import ModulatedParticleElement
>>> ModulatedParticleElement('bead').ds
1.8

Source code in src/skcomplex/systems/physics/soft_matter/_diffusion.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(role='position')] = 50.0,
    s0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    ds: Annotated[float, ParameterSpec()] = 1.8,
    L: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.s0 = s0
    self.ds = ds
    self.L = L

OpticalTrapInteraction(label, *, target, k=1.0, amplitude=0.0, frequency=0.0, blink_period=None, duty=0.5)

Bases: GradientInteraction

Harmonic optical trap on overdamped particles, the potential sum_d k_d (x_d - c_d(t))^2 / 2 of every axis.

The center oscillates as c(t) = amplitude sin(2 pi frequency t), and a blinking trap is on for the first duty fraction of every blink_period and off for the rest. The force divides by the friction gamma of the target.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the particle or the group.	required
k	Annotated[ArrayLike, ParameterSpec(vector=True)]	The stiffness along every axis.	1.0
amplitude	Annotated[ArrayLike, ParameterSpec(vector=True)]	The amplitude of the oscillation of the center along every axis.	0.0
frequency	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The frequency of the oscillation of the center.	0.0
blink_period	float | None	The period of a blinking trap, or None for a trap that is always on.	None
duty	float	The fraction of every blinking period during which the trap is on.	0.5

Examples:

>>> from skcomplex.systems.physics.soft_matter import OpticalTrapInteraction
>>> OpticalTrapInteraction('tweezers', target='bead', k=[1.0, 0.25]).k
[1.0, 0.25]

Source code in src/skcomplex/systems/physics/soft_matter/_diffusion.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    k: Annotated[npt.ArrayLike, ParameterSpec(vector=True)] = 1.0,
    amplitude: Annotated[npt.ArrayLike, ParameterSpec(vector=True)] = 0.0,
    frequency: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    blink_period: float | None = None,
    duty: float = 0.5,
) -> None:
    super().__init__(label, target=target, states='x', divisor='gamma')
    self.k = k
    self.amplitude = amplitude
    self.frequency = frequency
    self.blink_period = blink_period
    self.duty = duty

ScaledWalkerElement(label, visual=None, *, x=0.0, K=1.0, alpha=0.5, t0=0.01)

Bases: Element

Walker at the position x whose diffusion coefficient K alpha (t + t0)^(alpha - 1) changes in time.

Parameters:

Name	Type	Description	Default
label	str	The label of the walker or the group.	required
visual	VisualSpec | None	How the walker is drawn, or None for none.	None
x	Annotated[float, StateSpec(role=position)]	The initial position.	0.0
K	Annotated[float, ParameterSpec(domain=open(0, oo))]	The generalized diffusion coefficient.	1.0
alpha	Annotated[float, ParameterSpec(domain=open(0, oo))]	The anomalous diffusion exponent.	0.5
t0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time offset of the diffusion coefficient.	0.01

Examples:

>>> from skcomplex.systems.physics.soft_matter import ScaledWalkerElement
>>> ScaledWalkerElement('walker', alpha=1.5).alpha
1.5

Source code in src/skcomplex/systems/physics/soft_matter/_diffusion.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(role='position')] = 0.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    t0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.K = K
    self.alpha = alpha
    self.t0 = t0

SpeckleFieldInteraction(label, *, target, strength=1.0, centers=(), intensities=(), lengths=())

Bases: GradientInteraction

Random optical speckle field of Gaussian hotspots, which drags particles up its gradient at u grad(I).

The potential is -u I(x), with the intensity I(x) = sum_i I_i exp(-|x - c_i|^2 / d_i^2) of hotspots whose centers, intensities relative to strength, and decay lengths are fixed when the field is created.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the particle or the group.	required
strength	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The product u I0 of the polarizability constant and the reference intensity.	1.0
centers	tuple[tuple[float, float], ...]	The centers of the hotspots.	()
intensities	tuple[float, ...]	The intensities of the hotspots in units of the reference intensity.	()
lengths	tuple[float, ...]	The decay lengths of the hotspots.	()

Examples:

>>> from skcomplex.systems.physics.soft_matter import SpeckleFieldInteraction
>>> SpeckleFieldInteraction(
...     'speckle', target='bead', centers=((0.0, 0.0),), intensities=(1.0,), lengths=(1.0,)
... ).strength
1.0

Source code in src/skcomplex/systems/physics/soft_matter/_diffusion.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    strength: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    centers: tuple[tuple[float, float], ...] = (),
    intensities: tuple[float, ...] = (),
    lengths: tuple[float, ...] = (),
) -> None:
    super().__init__(label, target=target, states='x')
    self.strength = strength
    self.centers = centers
    self.intensities = intensities
    self.lengths = lengths

ThermalNoiseInteraction(label, *, target, state='x')

Bases: RateInteraction

Thermal noise of overdamped particles, the random force sqrt(2 gamma temperature) dW divided by gamma.

The friction gamma and the thermal energy temperature are parameters of the target, so every component of the position diffuses with the coefficient temperature / gamma.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the particle or the group.	required
state	str	The name of the state that holds the positions.	'x'

Examples:

>>> from skcomplex.systems.physics.common.particles import ThermalNoiseInteraction
>>> ThermalNoiseInteraction('bath', target='beads').divisor
'gamma'

Source code in src/skcomplex/systems/physics/common/particles/_thermal_noise.py

def __init__(self, label: str, *, target: Annotated[str, RoleSpec()], state: str = 'x') -> None:
    super().__init__(label, target=target, state=state, divisor='gamma')

TiltedWashboardInteraction(label, *, target, U0=1.0, L=1.0, f=0.0)

Bases: GradientInteraction

Tilted periodic potential U(x) - f x, with U(x) = U0 cos(2 pi x / L), on overdamped particles on a line.

The force divides by the friction gamma of the target.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the particle or the group.	required
U0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the periodic potential, half its barrier height.	1.0
L	Annotated[float, ParameterSpec(domain=open(0, oo))]	The period of the potential.	1.0
f	Annotated[float, ParameterSpec()]	The constant non-conservative force that tilts the potential.	0.0

Examples:

>>> from skcomplex.systems.physics.soft_matter import TiltedWashboardInteraction
>>> TiltedWashboardInteraction('washboard', target='beads', f=0.5).f
0.5

Source code in src/skcomplex/systems/physics/soft_matter/_ratchets.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    U0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    L: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    f: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, target=target, states='x', divisor='gamma')
    self.U0 = U0
    self.L = L
    self.f = f

TransientWalkerElement(label, visual=VisualSpec('points', {'position': 'x', 'color_by': 'D', 'color_range': (0.0, 1.0)}), *, x=0.0, D=0.5, sigma=0.5, gamma=1.0)

Bases: Element

Brownian walker whose diffusion coefficient D is redrawn after waiting times of mean D^(-gamma).

The coefficient takes values of density sigma D^(sigma - 1) on (0, 1), the beta distribution B(sigma, 1), so slow walkers wait long before they change.

Parameters:

Name	Type	Description	Default
label	str	The label of the walker or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	VisualSpec('points', {'position': 'x', 'color_by': 'D', 'color_range': (0.0, 1.0)})
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
D	Annotated[float, StateSpec(domain=open(0, oo))]	The initial diffusion coefficient.	0.5
sigma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The exponent of the density of diffusion coefficients near zero.	0.5
gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The exponent of the mean waiting time D^(-gamma).	1.0

Examples:

>>> from skcomplex.systems.physics.soft_matter import TransientWalkerElement
>>> TransientWalkerElement('walker', gamma=2.0).gamma
2.0

Source code in src/skcomplex/systems/physics/soft_matter/_anomalous.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('points', {'position': 'x', 'color_by': 'D', 'color_range': (0.0, 1.0)}),
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    D: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 0.5,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.D = D
    self.sigma = sigma
    self.gamma = gamma

TrapSiteElement(label, visual=None, *, occupied=0.0, visited=0.0, clock=0.0, tau=1.0)

Bases: Element

Site of a line of traps, which holds a walker for its trapping time tau before the walker hops to the right.

The walker on a site counts the steps it has spent there with clock and hops to the next site after tau steps, and every site keeps a record of whether a walker has visited it. A walker hops only onto a site of positive trapping time, so on a lattice that reads -1 outside its edges a walker on the last site stays there.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
occupied	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	Whether a walker is on the site.	0.0
visited	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	Whether a walker has been on the site.	0.0
clock	Annotated[float, StateSpec(domain=NN)]	The number of steps the walker has spent on the site.	0.0
tau	Annotated[float, StateSpec(domain=open(0, oo))]	The trapping time of the site.	1.0

Examples:

>>> from skcomplex.systems.physics.soft_matter import TrapSiteElement
>>> TrapSiteElement('site', tau=3.0).tau
3.0

Source code in src/skcomplex/systems/physics/soft_matter/_traps.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    occupied: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
    visited: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
    clock: Annotated[float, StateSpec(domain=NN)] = 0.0,
    tau: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.occupied = occupied
    self.visited = visited
    self.clock = clock
    self.tau = tau

TrappedSphereElement(label, visual=SPHERES, *, x=0.0, center=0.0)

Bases: Element

Colloidal sphere at the position x, held by an optical trap centred at center.

The sphere has no friction of its own. A mobility interaction moves it through the fluid it shares with the other spheres, and its trap is one of the forces the mobility turns into velocities.

Parameters:

Name	Type	Description	Default
label	str	The label of the sphere or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	SPHERES
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The initial position.	0.0
center	Annotated[ArrayLike, ParameterSpec(vector=True)]	The centre of the trap of the sphere.	0.0

Examples:

>>> from skcomplex.systems.physics.soft_matter import TrappedSphereElement
>>> TrappedSphereElement('bead', center=[3.0, 0.0]).center
[3.0, 0.0]

Source code in src/skcomplex/systems/physics/soft_matter/_coupled_traps.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = SPHERES,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    center: Annotated[npt.ArrayLike, ParameterSpec(vector=True)] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.center = center

WallParticleElement(label, visual=None, *, z=1.0, D0=1.0, a=1.0, drift=0.2)

Bases: Element

Brownian particle at the distance z from a wall, which reflects it, and whose diffusion slows near the wall.

Parameters:

Name	Type	Description	Default
label	str	The label of the particle.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
z	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The distance from the wall.	1.0
D0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The diffusion coefficient far from the wall.	1.0
a	Annotated[float, ParameterSpec(domain=open(0, oo))]	The distance over which the wall slows diffusion.	1.0
drift	Annotated[float, ParameterSpec()]	The settling velocity far from the wall.	0.2

Examples:

>>> from skcomplex.systems.physics.soft_matter import WallParticleElement
>>> WallParticleElement('bead').a
1.0

Source code in src/skcomplex/systems/physics/soft_matter/_brownian.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    z: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    D0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    a: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    drift: Annotated[float, ParameterSpec()] = 0.2,
) -> None:
    super().__init__(label, visual=visual)
    self.z = z
    self.D0 = D0
    self.a = a
    self.drift = drift

anisotropic_brownian_particles(n_particles=200, *, D_par=1.0, D_perp=1.0, D_r=0.75, D_spin=0.75, D_c=0.0, random_state=0)

Return overdamped axisymmetric Brownian particles in three dimensions that translate and rotate.

Every particle diffuses with D_par along its axis u and with D_perp across it, its axis diffuses on the unit sphere with D_r, so that <u(t) . u(0)> = exp(-2 D_r t), and it spins about its axis with D_spin. The defaults are a Janus sphere of unit radius, D_par = D_perp = D_t and D_r = 3 D_t / 4, whose orientation is visible although its diffusion is isotropic. A prolate ellipsoid has D_par > D_perp from Perrin's factors, so its displacement over times shorter than 1 / D_r is elongated along its axis. A chiral body, such as a spiral, has a coupling D_c of the translation along its axis and the rotation about it, <dr_par dpsi> = 2 D_c dt. These are Challenges 5.1 to 5.3 of Argun et al. The axis is held as a unit vector and rescaled to unit length after every step, and the particles start at the origin with axes drawn uniformly on the sphere.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	200
D_par	float	The translational diffusion coefficient along the axis.	1.0
D_perp	float	The translational diffusion coefficient across the axis.	1.0
D_r	float	The rotational diffusion coefficient of the axis.	0.75
D_spin	float	The rotational diffusion coefficient about the axis.	0.75
D_c	float	The coupling of the translation along the axis and the rotation about it.	0.0
random_state	int | Generator | None	The seed or generator of the initial axes.	0

Returns:

Name	Type	Description
system	System	The system, with the positions particles.x, particles.y, particles.z, the axes particles.ux, particles.uy, particles.uz, and the spins particles.psi.

Raises:

Type	Description
ParameterError	The coupling exceeds sqrt(D_par D_spin).

Examples:

>>> from skcomplex.systems.physics.soft_matter import anisotropic_brownian_particles
>>> anisotropic_brownian_particles(5, D_par=1.5, D_perp=0.8).elements[0].values['D_perp']
0.8

Source code in src/skcomplex/systems/physics/soft_matter/_rigid_bodies.py

@register_system(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'rotational diffusion', 'three dimensions', 'anisotropy', 'chirality'),
    references=('Perrin 1934', 'Han et al. 2006', 'Argun et al. 2021'),
    defaults={'until': 10.0, 'dt': 0.05, 'solver': EulerMaruyamaSolver(dt=0.001)},
)
def anisotropic_brownian_particles(
    n_particles: int = 200,
    *,
    D_par: float = 1.0,
    D_perp: float = 1.0,
    D_r: float = 0.75,
    D_spin: float = 0.75,
    D_c: float = 0.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return overdamped axisymmetric Brownian particles in three dimensions that translate and rotate.

    Every particle diffuses with `D_par` along its axis `u` and with `D_perp` across it, its axis diffuses on the unit
    sphere with `D_r`, so that `<u(t) . u(0)> = exp(-2 D_r t)`, and it spins about its axis with `D_spin`. The defaults
    are a Janus sphere of unit radius, `D_par = D_perp = D_t` and `D_r = 3 D_t / 4`, whose orientation is visible
    although its diffusion is isotropic. A prolate ellipsoid has `D_par > D_perp` from Perrin's factors, so its
    displacement over times shorter than `1 / D_r` is elongated along its axis. A chiral body, such as a spiral, has a
    coupling `D_c` of the translation along its axis and the rotation about it, `<dr_par dpsi> = 2 D_c dt`. These are
    Challenges 5.1 to 5.3 of Argun et al. The axis is held as a unit vector and rescaled to unit length after every
    step, and the particles start at the origin with axes drawn uniformly on the sphere.

    Args:
        n_particles:
            The number of particles.

        D_par:
            The translational diffusion coefficient along the axis.

        D_perp:
            The translational diffusion coefficient across the axis.

        D_r:
            The rotational diffusion coefficient of the axis.

        D_spin:
            The rotational diffusion coefficient about the axis.

        D_c:
            The coupling of the translation along the axis and the rotation about it.

        random_state:
            The seed or generator of the initial axes.

    Returns:
        system:
            The system, with the positions `particles.x`, `particles.y`, `particles.z`, the axes `particles.ux`,
            `particles.uy`, `particles.uz`, and the spins `particles.psi`.

    Raises:
        ParameterError:
            The coupling exceeds `sqrt(D_par D_spin)`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import anisotropic_brownian_particles
        >>> anisotropic_brownian_particles(5, D_par=1.5, D_perp=0.8).elements[0].values['D_perp']
        0.8
    """
    if D_c**2 > D_par * D_spin:
        error_msg = f'The coupling {D_c!r} exceeds sqrt(D_par D_spin), so the diffusion tensor is not positive.'
        raise ParameterError(error_msg)
    axes = np.random.default_rng(random_state).normal(size=(n_particles, 3))
    axes /= np.linalg.norm(axes, axis=1, keepdims=True)
    values = {
        'ux': axes[:, 0],
        'uy': axes[:, 1],
        'uz': axes[:, 2],
        'D_par': D_par,
        'D_perp': D_perp,
        'D_r': D_r,
        'D_spin': D_spin,
        'D_c': D_c,
    }
    particles = GroupElement('particles', element=AnisotropicParticleElement, size=n_particles, values=values)
    return System(elements=[particles])

annealed_transient_time_motion(n_walkers=500, *, sigma=0.5, gamma=1.0, n_dim=2, random_state=0)

Return annealed transient time motion, Brownian walkers whose diffusion coefficient changes at random times.

Every walker diffuses with a coefficient D of density sigma D^(sigma - 1) on (0, 1) and keeps it for an exponential time of mean D^(-gamma), so slow walkers stay slow for long. For sigma < gamma the ensemble mean squared displacement grows as t^(sigma / gamma).

Parameters:

Name	Type	Description	Default
n_walkers	int	The number of independent walkers.	500
sigma	float	The exponent of the density of diffusion coefficients near zero.	0.5
gamma	float	The exponent of the mean waiting time D^(-gamma).	1.0
n_dim	int	The number of dimensions.	2
random_state	int | Generator | None	The seed or generator of the first diffusion coefficients.	0

Returns:

Name	Type	Description
system	System	The system, with the positions walkers.x, the coefficients walkers.D, and the change count changes.renewals.

Examples:

>>> from skcomplex.systems.physics.soft_matter import annealed_transient_time_motion
>>> annealed_transient_time_motion(10, gamma=2.0).elements[0].values['gamma']
2.0

Source code in src/skcomplex/systems/physics/soft_matter/_anomalous.py

@register_system(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'renewal process', 'heterogeneity', 'sde'),
    references=('Argun et al. 2021', 'Massignan et al. 2014'),
    defaults={'until': 500.0, 'dt': 2.5, 'solver': EulerMaruyamaSolver(dt=0.5)},
)
def annealed_transient_time_motion(
    n_walkers: int = 500,
    *,
    sigma: float = 0.5,
    gamma: float = 1.0,
    n_dim: int = 2,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return annealed transient time motion, Brownian walkers whose diffusion coefficient changes at random times.

    Every walker diffuses with a coefficient `D` of density `sigma D^(sigma - 1)` on `(0, 1)` and keeps it for an
    exponential time of mean `D^(-gamma)`, so slow walkers stay slow for long. For `sigma < gamma` the ensemble mean
    squared displacement grows as `t^(sigma / gamma)`.

    Args:
        n_walkers:
            The number of independent walkers.

        sigma:
            The exponent of the density of diffusion coefficients near zero.

        gamma:
            The exponent of the mean waiting time `D^(-gamma)`.

        n_dim:
            The number of dimensions.

        random_state:
            The seed or generator of the first diffusion coefficients.

    Returns:
        system:
            The system, with the positions `walkers.x`, the coefficients `walkers.D`, and the change count
            `changes.renewals`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import annealed_transient_time_motion
        >>> annealed_transient_time_motion(10, gamma=2.0).elements[0].values['gamma']
        2.0
    """
    coefficients = np.random.default_rng(random_state).beta(sigma, 1.0, n_walkers)
    walkers = GroupElement(
        'walkers',
        element=TransientWalkerElement,
        size=n_walkers,
        values={'x': np.zeros((n_walkers, n_dim)), 'D': coefficients, 'sigma': sigma, 'gamma': gamma},
    )
    noise = RateInteraction('noise', target='walkers', noise=lambda s: sqrt(2 * s.target.D))
    changes = JumpInteraction(
        'changes',
        target='walkers',
        waiting=lambda s: draw('wait', 'expon', scale=s.target.D ** (-s.target.gamma)),
        reset={'D': lambda s: draw('coefficient', 'beta', a=s.target.sigma, b=1)},
    )
    return System(elements=[walkers], interactions=[noise, changes], space=ContinuousSpace(n_dim=n_dim))

biased_trap_walk(n_walkers=100, *, n_sites=300, alpha=0.5, tau0=1.0, random_state=0)

Return walkers that hop always to the right through quenched traps with broadly distributed trapping times.

Every walker moves along its own line of sites and stays tau_i steps on site i before it hops to site i + 1. The trapping times are drawn once, as ceil(tau0 u^(-1 / alpha)) for uniform numbers u, so their distribution decays as tau0^alpha / tau^(1 + alpha). For alpha < 1 their mean diverges, the time to cross N sites grows as N^(1 / alpha), and the displacement grows subdiffusively as x ~ (t / tau0)^alpha, with the mean <x> = (t / tau0)^alpha sin(pi alpha) / (pi alpha) at long times. The walkers start on the first site.

Parameters:

Name	Type	Description	Default
n_walkers	int	The number of walkers, each in its own landscape.	100
n_sites	int	The number of sites of every landscape.	300
alpha	float	The exponent of the distribution of trapping times.	0.5
tau0	float	The smallest trapping time.	1.0
random_state	int | Generator | None	The seed or generator of the trapping times.	0

Returns:

Name	Type	Description
system	System	The system, with the lattice traps of shape (n_walkers, n_sites) and the mean displacement position.

Examples:

>>> from skcomplex.systems.physics.soft_matter import biased_trap_walk
>>> biased_trap_walk(4, n_sites=10).elements[0].prepare().size_
40

Source code in src/skcomplex/systems/physics/soft_matter/_traps.py

@register_system(
    field='physics/soft_matter',
    tags=('lattice', 'anomalous diffusion', 'disorder', 'traps', 'subdiffusion'),
    references=('Bertin 2021', 'Bouchaud and Georges 1990'),
    defaults={'until': 2000, 'dt': 10},
    observables=_trap_observables,
)
def biased_trap_walk(
    n_walkers: int = 100,
    *,
    n_sites: int = 300,
    alpha: float = 0.5,
    tau0: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return walkers that hop always to the right through quenched traps with broadly distributed trapping times.

    Every walker moves along its own line of sites and stays `tau_i` steps on site `i` before it hops to site
    `i + 1`. The trapping times are drawn once, as `ceil(tau0 u^(-1 / alpha))` for uniform numbers `u`, so their
    distribution decays as `tau0^alpha / tau^(1 + alpha)`. For `alpha < 1` their mean diverges, the time to cross `N`
    sites grows as `N^(1 / alpha)`, and the displacement grows subdiffusively as `x ~ (t / tau0)^alpha`, with the
    mean `<x> = (t / tau0)^alpha sin(pi alpha) / (pi alpha)` at long times. The walkers start on the first site.

    Args:
        n_walkers:
            The number of walkers, each in its own landscape.

        n_sites:
            The number of sites of every landscape.

        alpha:
            The exponent of the distribution of trapping times.

        tau0:
            The smallest trapping time.

        random_state:
            The seed or generator of the trapping times.

    Returns:
        system:
            The system, with the lattice `traps` of shape `(n_walkers, n_sites)` and the mean displacement `position`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import biased_trap_walk
        >>> biased_trap_walk(4, n_sites=10).elements[0].prepare().size_
        40
    """
    shape = (n_walkers, n_sites)
    tau = np.ceil(tau0 * (1 - np.random.default_rng(random_state).random(shape)) ** (-1 / alpha))
    start = np.zeros(shape)
    start[:, 0] = 1.0
    traps = LatticeElement(
        'traps',
        element=TrapSiteElement,
        space=LatticeSpace(shape, boundary='fixed', fill=-1.0),
        values={'occupied': start, 'visited': start, 'tau': tau},
        visual=VISITED,
    )
    hops = LocalRuleInteraction(
        'hops',
        target='traps',
        state=('occupied', 'visited', 'clock'),
        rule=_hop,
        neighborhood=ALONG_ROW,
        reads=('occupied', 'visited', 'clock', 'tau'),
    )
    return System(elements=[traps], interactions=[hops])

brownian_motion(n_particles=500, *, inertial=True, m=1.0, gamma=1.0, kT=1.0, n_dim=1, random_state=0)

Return free Brownian particles with inertia, m x'' = -gamma x' + sqrt(2 kT gamma) W, or without it.

With inertia the motion is ballistic, <x^2> = kT t^2 / m, below the momentum relaxation time m / gamma, and diffusive, <x^2> = 2 kT t / gamma per axis, above it. Without inertia it is diffusive at every time.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	500
inertial	bool	Whether the particles have inertia.	True
m	float	The mass, used with inertia.	1.0
gamma	float	The friction coefficient.	1.0
kT	float	The thermal energy.	1.0
n_dim	int	The number of dimensions.	1
random_state	int | Generator | None	The seed or generator of the thermal initial velocities.	0

Returns:

Name	Type	Description
system	System	The system, with the positions beads.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import brownian_motion
>>> brownian_motion(10, inertial=False).elements[0].element.__name__
'BrownianParticleElement'

Source code in src/skcomplex/systems/physics/soft_matter/_diffusion.py

@register_system(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'diffusion', 'inertia'),
    references=('Argun et al. 2021', 'Li et al. 2010', 'Volpe and Volpe 2013'),
    defaults={'until': 10.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.005)},
)
def brownian_motion(
    n_particles: int = 500,
    *,
    inertial: bool = True,
    m: float = 1.0,
    gamma: float = 1.0,
    kT: float = 1.0,
    n_dim: int = 1,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return free Brownian particles with inertia, `m x'' = -gamma x' + sqrt(2 kT gamma) W`, or without it.

    With inertia the motion is ballistic, `<x^2> = kT t^2 / m`, below the momentum relaxation time `m / gamma`, and
    diffusive, `<x^2> = 2 kT t / gamma` per axis, above it. Without inertia it is diffusive at every time.

    Args:
        n_particles:
            The number of particles.

        inertial:
            Whether the particles have inertia.

        m:
            The mass, used with inertia.

        gamma:
            The friction coefficient.

        kT:
            The thermal energy.

        n_dim:
            The number of dimensions.

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

    Returns:
        system:
            The system, with the positions `beads.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import brownian_motion
        >>> brownian_motion(10, inertial=False).elements[0].element.__name__
        'BrownianParticleElement'
    """
    space = ContinuousSpace(n_dim=n_dim)
    if not inertial:
        beads = GroupElement(
            'beads',
            element=BrownianParticleElement,
            size=n_particles,
            values={'x': np.zeros((n_particles, n_dim)), 'gamma': gamma, 'temperature': kT},
        )
        return System(elements=[beads], interactions=[ThermalNoiseInteraction('bath', target='beads')], space=space)
    velocities = np.random.default_rng(random_state).normal(0.0, math.sqrt(kT / m), (n_particles, n_dim))
    beads = GroupElement(
        'beads',
        element=MassPointElement,
        size=n_particles,
        values={'x': np.zeros((n_particles, n_dim)), 'x_d1': velocities, 'm': m},
        visual=POINTS,
    )
    bath = LangevinThermostatInteraction('bath', target='beads', temperature=kT, gamma=gamma)
    return System(elements=[beads], interactions=[bath], space=space)

brownian_near_wall(z=1.0, *, D0=1.0, a=1.0, drift=0.2, convention='anti-ito')

Return a Brownian particle whose diffusion slows near a wall, D(z) = D0 z² / (z² + a²), as it settles on it.

The diffusion is a smooth form of the hydrodynamic hindrance that vanishes at the wall, and the drift is the settling velocity scaled by the same mobility. Diffusion that depends on the position makes the noise multiplicative, and the anti-Itô convention gives the equilibrium of physical particles. The laws are symmetric about the wall at z = 0, which reflects the particle at every frame.

Parameters:

Name	Type	Description	Default
z	float	The initial distance from the wall.	1.0
D0	float	The diffusion coefficient far from the wall.	1.0
a	float	The distance over which the wall slows diffusion, such as the particle radius.	1.0
drift	float	The settling velocity far from the wall.	0.2
convention	str	The convention of the noise, 'ito', 'stratonovich', or 'anti-ito'.	'anti-ito'

Returns:

Name	Type	Description
system	System	The system, with the distance particle.z.

Examples:

>>> from skcomplex.systems.physics.soft_matter import brownian_near_wall
>>> brownian_near_wall(convention='ito').interactions[0].convention
'ito'

Source code in src/skcomplex/systems/physics/soft_matter/_brownian.py

@register_system(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'multiplicative noise', 'hydrodynamics'),
    references=('Brenner 1961', 'Volpe and Wehr 2016'),
    defaults={'until': 10.0, 'dt': 0.1},
)
def brownian_near_wall(
    z: float = 1.0,
    *,
    D0: float = 1.0,
    a: float = 1.0,
    drift: float = 0.2,
    convention: str = 'anti-ito',
) -> System:
    """Return a Brownian particle whose diffusion slows near a wall, `D(z) = D0 z² / (z² + a²)`, as it settles on it.

    The diffusion is a smooth form of the hydrodynamic hindrance that vanishes at the wall, and the drift is the
    settling velocity scaled by the same mobility. Diffusion that depends on the position makes the noise
    multiplicative, and the anti-Itô convention gives the equilibrium of physical particles. The laws are symmetric
    about the wall at `z = 0`, which reflects the particle at every frame.

    Args:
        z:
            The initial distance from the wall.

        D0:
            The diffusion coefficient far from the wall.

        a:
            The distance over which the wall slows diffusion, such as the particle radius.

        drift:
            The settling velocity far from the wall.

        convention:
            The convention of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.

    Returns:
        system:
            The system, with the distance `particle.z`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import brownian_near_wall
        >>> brownian_near_wall(convention='ito').interactions[0].convention
        'ito'
    """
    mobility = RateInteraction(
        'mobility',
        target='particle',
        state='z',
        rate=lambda s: -s.target.drift * s.target.z * abs(s.target.z) / (s.target.z**2 + s.target.a**2),
        noise=lambda s: sqrt(2 * s.target.D0) * s.target.z / sqrt(s.target.z**2 + s.target.a**2),
        convention=convention,
    )
    particle = WallParticleElement('particle', z=z, D0=D0, a=a, drift=drift)
    return System(elements=[particle], interactions=[mobility])

colored_noise_langevin(n_particles=500, *, gamma=1.0, D=1.0, tau=1.0, random_state=0)

Return particles whose velocity dv/dt = -gamma v + xi is driven by Ornstein-Uhlenbeck noise xi.

The noise follows tau dxi = -xi dt + sqrt(2 D) dW, so its correlation is <xi(t) xi(t')> = D / tau exp(-|t - t'| / tau), and it tends to the white noise of intensity 2 D as tau goes to 0. The stationary velocity variance is <v^2> = D / (gamma (1 + gamma tau)), which is D / gamma in the white noise limit. The velocities start at 0 and the noise starts in its stationary state. Solver steps must be short against tau.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of independent particles.	500
gamma	float	The friction coefficient.	1.0
D	float	The noise intensity.	1.0
tau	float	The correlation time of the noise.	1.0
random_state	int | Generator | None	The seed or generator of the initial noise values.	0

Returns:

Name	Type	Description
system	System	The system, with the velocities particles.v and the noise particles.xi.

Examples:

>>> from skcomplex.systems.physics.soft_matter import colored_noise_langevin
>>> colored_noise_langevin(10, tau=0.5).elements[0].values['tau']
0.5

Source code in src/skcomplex/systems/physics/soft_matter/_colored_noise.py

@register_system(
    field='physics/soft_matter',
    tags=('noise', 'sde', 'colored noise', 'langevin equation'),
    references=('Bertin 2021', 'Hänggi and Jung 1995'),
    defaults={'until': 50.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def colored_noise_langevin(
    n_particles: int = 500,
    *,
    gamma: float = 1.0,
    D: float = 1.0,
    tau: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return particles whose velocity `dv/dt = -gamma v + xi` is driven by Ornstein-Uhlenbeck noise `xi`.

    The noise follows `tau dxi = -xi dt + sqrt(2 D) dW`, so its correlation is `<xi(t) xi(t')> = D / tau
    exp(-|t - t'| / tau)`, and it tends to the white noise of intensity `2 D` as `tau` goes to 0. The stationary
    velocity variance is `<v^2> = D / (gamma (1 + gamma tau))`, which is `D / gamma` in the white noise limit. The
    velocities start at 0 and the noise starts in its stationary state. Solver steps must be short against `tau`.

    Args:
        n_particles:
            The number of independent particles.

        gamma:
            The friction coefficient.

        D:
            The noise intensity.

        tau:
            The correlation time of the noise.

        random_state:
            The seed or generator of the initial noise values.

    Returns:
        system:
            The system, with the velocities `particles.v` and the noise `particles.xi`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import colored_noise_langevin
        >>> colored_noise_langevin(10, tau=0.5).elements[0].values['tau']
        0.5
    """
    noise = np.random.default_rng(random_state).normal(0.0, math.sqrt(D / tau), n_particles)
    particles = GroupElement(
        'particles',
        element=ColoredNoiseParticleElement,
        size=n_particles,
        values={'xi': noise, 'gamma': gamma, 'D': D, 'tau': tau},
    )
    drag = RateInteraction(
        'drag',
        target='particles',
        state='v',
        rate=lambda s: -s.target.gamma * s.target.v + s.target.xi,
    )
    noise = RelaxationInteraction(
        'noise',
        target='particles',
        state='xi',
        time=lambda s: s.target.tau,
        noise=lambda s: sqrt(2 * s.target.D) / s.target.tau,
    )
    return System(elements=[particles], interactions=[drag, noise])

continuous_time_random_walk(n_walkers=500, *, alpha=0.5, tau0=1.0, sigma=1.0, n_dim=2)

Return the continuous-time random walk, Gaussian jumps after waiting times with the tail tau^(-1 - alpha).

For alpha < 1 the waiting times have no mean, so the number of jumps grows as t^alpha and the ensemble mean squared displacement grows as t^alpha, slower than Brownian motion. The walkers renew on their own clocks.

Parameters:

Name	Type	Description	Default
n_walkers	int	The number of independent walkers.	500
alpha	float	The tail exponent of the waiting times.	0.5
tau0	float	The shortest waiting time.	1.0
sigma	float	The standard deviation of the jumps along every axis.	1.0
n_dim	int	The number of dimensions.	2

Returns:

Name	Type	Description
system	System	The system, with the positions walkers.x and the jump count jumps.renewals.

Examples:

>>> from skcomplex.systems.physics.soft_matter import continuous_time_random_walk
>>> continuous_time_random_walk(10, alpha=0.7).elements[0].values['alpha']
0.7

Source code in src/skcomplex/systems/physics/soft_matter/_anomalous.py

@register_system(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'renewal process', 'subdiffusion', 'random walk'),
    references=('Argun et al. 2021', 'Metzler and Klafter 2000', 'Montroll and Weiss 1965', 'Bertin 2021'),
    defaults={'until': 1000.0, 'dt': 5.0},
)
def continuous_time_random_walk(
    n_walkers: int = 500,
    *,
    alpha: float = 0.5,
    tau0: float = 1.0,
    sigma: float = 1.0,
    n_dim: int = 2,
) -> System:
    """Return the continuous-time random walk, Gaussian jumps after waiting times with the tail `tau^(-1 - alpha)`.

    For `alpha < 1` the waiting times have no mean, so the number of jumps grows as `t^alpha` and the ensemble mean
    squared displacement grows as `t^alpha`, slower than Brownian motion. The walkers renew on their own clocks.

    Args:
        n_walkers:
            The number of independent walkers.

        alpha:
            The tail exponent of the waiting times.

        tau0:
            The shortest waiting time.

        sigma:
            The standard deviation of the jumps along every axis.

        n_dim:
            The number of dimensions.

    Returns:
        system:
            The system, with the positions `walkers.x` and the jump count `jumps.renewals`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import continuous_time_random_walk
        >>> continuous_time_random_walk(10, alpha=0.7).elements[0].values['alpha']
        0.7
    """
    walkers = GroupElement(
        'walkers',
        element=ContinuousTimeWalkerElement,
        size=n_walkers,
        values={'x': np.zeros((n_walkers, n_dim)), 'alpha': alpha, 'tau0': tau0, 'sigma': sigma},
    )
    jumps = JumpInteraction(
        'jumps',
        target='walkers',
        waiting=lambda s: draw('wait', 'pareto', b=s.target.alpha, scale=s.target.tau0),
        reset={
            'x': lambda s: [item + normal(f'jump{index}', 0, s.target.sigma) for index, item in enumerate(s.target.x)],
        },
    )
    return System(elements=[walkers], interactions=[jumps], space=ContinuousSpace(n_dim=n_dim))

diffusion_in_a_speckle_field(n_particles=100, *, n_hotspots=50, size=100.0, strength=30.0, length=10.0, D=1.0, random_state=0)

Return Brownian particles in a random speckle field, trapped at hotspots for longer as the field grows stronger.

The hotspots sit at random places in a square, with intensities r^(-1/2) and decay lengths length r^(-1/2) for independent uniform numbers r. The particles start at random places in the middle half of the square.

The default strength 30 draws the particles into the hotspots over the first ten time units and lets them escape now and then, so from time 10 to 60 their mean squared displacement grows as t^0.36, a tenth of free diffusion. A strength of 100 or more traps every particle within the first frames and leaves only the rattling inside one hotspot, with a flat mean squared displacement.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of independent particles.	100
n_hotspots	int	The number of hotspots.	50
size	float	The side of the square of the hotspots.	100.0
strength	float	The product u I0 of the polarizability constant and the reference intensity.	30.0
length	float	The reference decay length of the hotspots.	10.0
D	float	The diffusion coefficient.	1.0
random_state	int | Generator | None	The seed or generator of the hotspots and the initial positions.	0

Returns:

Name	Type	Description
system	System	The system, with the positions beads.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import diffusion_in_a_speckle_field
>>> len(diffusion_in_a_speckle_field(5, n_hotspots=3).interactions[0].centers)
3

Source code in src/skcomplex/systems/physics/soft_matter/_diffusion.py

@register_system(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'sde', 'optical forces', 'disorder'),
    references=('Argun et al. 2021', 'Volpe, Volpe, and Gigan 2014'),
    defaults={'until': 60.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def diffusion_in_a_speckle_field(
    n_particles: int = 100,
    *,
    n_hotspots: int = 50,
    size: float = 100.0,
    strength: float = 30.0,
    length: float = 10.0,
    D: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return Brownian particles in a random speckle field, trapped at hotspots for longer as the field grows stronger.

    The hotspots sit at random places in a square, with intensities `r^(-1/2)` and decay lengths `length r^(-1/2)` for
    independent uniform numbers `r`. The particles start at random places in the middle half of the square.

    The default strength 30 draws the particles into the hotspots over the first ten time units and lets them escape
    now and then, so from time 10 to 60 their mean squared displacement grows as `t^0.36`, a tenth of free diffusion.
    A strength of 100 or more traps every particle within the first frames and leaves only the rattling inside one
    hotspot, with a flat mean squared displacement.

    Args:
        n_particles:
            The number of independent particles.

        n_hotspots:
            The number of hotspots.

        size:
            The side of the square of the hotspots.

        strength:
            The product `u I0` of the polarizability constant and the reference intensity.

        length:
            The reference decay length of the hotspots.

        D:
            The diffusion coefficient.

        random_state:
            The seed or generator of the hotspots and the initial positions.

    Returns:
        system:
            The system, with the positions `beads.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import diffusion_in_a_speckle_field
        >>> len(diffusion_in_a_speckle_field(5, n_hotspots=3).interactions[0].centers)
        3
    """
    rng = np.random.default_rng(random_state)
    centers = rng.uniform(0.0, size, (n_hotspots, 2))
    intensities = (1 - rng.random(n_hotspots)) ** -0.5
    lengths = length * (1 - rng.random(n_hotspots)) ** -0.5
    speckle = SpeckleFieldInteraction(
        'speckle',
        target='beads',
        strength=strength,
        centers=tuple((float(x), float(y)) for x, y in centers),
        intensities=tuple(float(item) for item in intensities),
        lengths=tuple(float(item) for item in lengths),
    )
    positions = rng.uniform(size / 4, 3 * size / 4, (n_particles, 2))
    beads = GroupElement(
        'beads',
        element=BrownianParticleElement,
        size=n_particles,
        values={'x': positions, 'gamma': 1.0, 'temperature': D},
    )
    bath = ThermalNoiseInteraction('bath', target='beads')
    return System(elements=[beads], interactions=[speckle, bath], space=ContinuousSpace(n_dim=2))

flashing_ratchet(n_particles=500, *, barrier=1.0, asymmetry=0.1, L=1.0, rate=1.0, duration=100.0, gamma=1.0, kT=0.05, random_state=0)

Return the on-off flashing ratchet, Brownian particles in a sawtooth potential switched on and off at random.

The potential switches at the rate rate, so on and off phases last 1 / rate on average. While it is off the particles diffuse freely, and while it is on they slide to the minima of the sawtooth. A potential that rises over a short fraction asymmetry of its period sends more particles over its near maximum than over its far one, so the particles drift with a mean velocity of the sign of 0.5 - asymmetry, and a symmetric potential gives no current. The switching times are drawn once, as a Poisson process over duration shared by every particle. The particles start at the minimum at the origin.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of independent particles.	500
barrier	float	The height of the sawtooth.	1.0
asymmetry	float	The fraction of the period over which the potential rises, 0.5 for a symmetric potential.	0.1
L	float	The period of the sawtooth.	1.0
rate	float	The rate at which the potential switches.	1.0
duration	float	The time span of the drawn switching times.	100.0
gamma	float	The friction coefficient.	1.0
kT	float	The thermal energy.	0.05
random_state	int | Generator | None	The seed or generator of the switching times.	0

Returns:

Name	Type	Description
system	System	The system, with the positions particles.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import flashing_ratchet
>>> len(flashing_ratchet(10, rate=2.0, duration=5.0).interactions[0].switches) > 0
True

Source code in src/skcomplex/systems/physics/soft_matter/_ratchets.py

@register_system(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'ratchet', 'nonequilibrium steady state', 'current'),
    references=('Bertin 2021', 'Ajdari and Prost 1992', 'Reimann 2002'),
    defaults={'until': 100.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.002)},
)
def flashing_ratchet(
    n_particles: int = 500,
    *,
    barrier: float = 1.0,
    asymmetry: float = 0.1,
    L: float = 1.0,
    rate: float = 1.0,
    duration: float = 100.0,
    gamma: float = 1.0,
    kT: float = 0.05,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the on-off flashing ratchet, Brownian particles in a sawtooth potential switched on and off at random.

    The potential switches at the rate `rate`, so on and off phases last `1 / rate` on average. While it is off the
    particles diffuse freely, and while it is on they slide to the minima of the sawtooth. A potential that rises over
    a short fraction `asymmetry` of its period sends more particles over its near maximum than over its far one, so
    the particles drift with a mean velocity of the sign of `0.5 - asymmetry`, and a symmetric potential gives no
    current. The switching times are drawn once, as a Poisson process over `duration` shared by every particle. The
    particles start at the minimum at the origin.

    Args:
        n_particles:
            The number of independent particles.

        barrier:
            The height of the sawtooth.

        asymmetry:
            The fraction of the period over which the potential rises, 0.5 for a symmetric potential.

        L:
            The period of the sawtooth.

        rate:
            The rate at which the potential switches.

        duration:
            The time span of the drawn switching times.

        gamma:
            The friction coefficient.

        kT:
            The thermal energy.

        random_state:
            The seed or generator of the switching times.

    Returns:
        system:
            The system, with the positions `particles.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import flashing_ratchet
        >>> len(flashing_ratchet(10, rate=2.0, duration=5.0).interactions[0].switches) > 0
        True
    """
    rng = np.random.default_rng(random_state)
    times = np.cumsum(rng.exponential(1 / rate, int(2 * rate * duration) + 10))
    particles = GroupElement(
        'particles',
        element=BrownianParticleElement,
        size=n_particles,
        values={'x': np.zeros((n_particles, 1)), 'gamma': gamma, 'temperature': kT},
    )
    sawtooth = FlashingSawtoothInteraction(
        'sawtooth',
        target='particles',
        barrier=barrier,
        L=L,
        asymmetry=asymmetry,
        switches=tuple(float(time) for time in times[times < duration]),
    )
    bath = ThermalNoiseInteraction('bath', target='particles')
    return System(elements=[particles], interactions=[sawtooth, bath], space=ContinuousSpace(n_dim=1))

fractional_brownian_motion(n_walkers=500, *, H=0.3, D=1.0, n_dim=2)

Return fractional Brownian motion, dx = sqrt(2 D) dB_H, with <x^2> = 2 D t^(2 H) along every axis.

The increments are correlated over the whole run: they anticorrelate for H < 1/2, which makes the walkers subdiffusive, and correlate for H > 1/2, which makes them superdiffusive.

Parameters:

Name	Type	Description	Default
n_walkers	int	The number of independent walkers.	500
H	float	The Hurst index.	0.3
D	float	The diffusion coefficient.	1.0
n_dim	int	The number of dimensions.	2

Returns:

Name	Type	Description
system	System	The system, with the positions walkers.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import fractional_brownian_motion
>>> fractional_brownian_motion(10, H=0.7).elements[0].values['H']
0.7

Source code in src/skcomplex/systems/physics/soft_matter/_anomalous.py

@register_system(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'fractional noise', 'long memory', 'sde'),
    references=('Argun et al. 2021', 'Mandelbrot and Van Ness 1968'),
    defaults={'until': 10.0, 'dt': 0.05, 'solver': FractionalEulerMaruyamaSolver(dt=0.01)},
)
def fractional_brownian_motion(
    n_walkers: int = 500,
    *,
    H: float = 0.3,
    D: float = 1.0,
    n_dim: int = 2,
) -> System:
    """Return fractional Brownian motion, `dx = sqrt(2 D) dB_H`, with `<x^2> = 2 D t^(2 H)` along every axis.

    The increments are correlated over the whole run: they anticorrelate for `H < 1/2`, which makes the walkers
    subdiffusive, and correlate for `H > 1/2`, which makes them superdiffusive.

    Args:
        n_walkers:
            The number of independent walkers.

        H:
            The Hurst index.

        D:
            The diffusion coefficient.

        n_dim:
            The number of dimensions.

    Returns:
        system:
            The system, with the positions `walkers.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import fractional_brownian_motion
        >>> fractional_brownian_motion(10, H=0.7).elements[0].values['H']
        0.7
    """
    walkers = GroupElement(
        'walkers',
        element=FractionalWalkerElement,
        size=n_walkers,
        values={'x': np.zeros((n_walkers, n_dim)), 'D': D, 'H': H},
    )
    noise = RateInteraction(
        'noise',
        target='walkers',
        noise=lambda s: sqrt(2 * s.target.D),
        hurst=lambda s: s.target.H,
    )
    return System(elements=[walkers], interactions=[noise], space=ContinuousSpace(n_dim=n_dim))

hydrodynamically_coupled_traps(separation=3.0, *, tensor='rotne-prager-yamakawa', root='cholesky', radius=1.0, viscosity=0.001, kT=0.0041, k=0.2, on_time=0.3, off_time=0.1, dt=0.005)

Return two Brownian spheres in blinking optical traps, whose free motion is correlated through the fluid.

The traps hold the spheres at a distance separation along x, and switch off for off_time after every on_time. While they are off the spheres diffuse freely, and the force on one moves the other through the mobility, so their displacements correlate as M_12 / M_11. That is about 3 a / (2 r) along the line of centres and 3 a / (4 r) across it. The units are micrometres, seconds, and piconewtons, so the defaults are spheres of radius 1 um in water at room temperature, and each step of the system lasts dt seconds.

Parameters:

Name	Type	Description	Default
separation	float	The distance between the centres of the traps.	3.0
tensor	str	'free-draining', 'oseen', or 'rotne-prager-yamakawa'.	'rotne-prager-yamakawa'
root	str	'cholesky' or 'lanczos'.	'cholesky'
radius	float	The radius of the spheres.	1.0
viscosity	float	The viscosity of the fluid.	0.001
kT	float	The thermal energy.	0.0041
k	float	The stiffness of the traps.	0.2
on_time	float	The time the traps stay on in every period.	0.3
off_time	float	The time the traps stay off in every period.	0.1
dt	float	The length of a step.	0.005

Returns:

Name	Type	Description
system	System	The system, with the positions spheres.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import hydrodynamically_coupled_traps
>>> hydrodynamically_coupled_traps(4.0).elements[0].values['center']
[[-2.0, 0.0], [2.0, 0.0]]

Source code in src/skcomplex/systems/physics/soft_matter/_coupled_traps.py

@register_system(
    field='physics/soft_matter',
    tags=('brownian motion', 'hydrodynamics', 'optical tweezers', 'blinking traps', 'correlated noise'),
    references=('Argun et al. 2021', 'Batchelor 1976', 'Rotne and Prager 1969', 'Magazzù et al. 2019'),
    defaults={'until': 20000, 'dt': 1},
)
def hydrodynamically_coupled_traps(
    separation: float = 3.0,
    *,
    tensor: str = 'rotne-prager-yamakawa',
    root: str = 'cholesky',
    radius: float = 1.0,
    viscosity: float = 1e-3,
    kT: float = 4.1e-3,
    k: float = 0.2,
    on_time: float = 0.3,
    off_time: float = 0.1,
    dt: float = 5e-3,
) -> System:
    """Return two Brownian spheres in blinking optical traps, whose free motion is correlated through the fluid.

    The traps hold the spheres at a distance `separation` along `x`, and switch off for `off_time` after every
    `on_time`. While they are off the spheres diffuse freely, and the force on one moves the other through the
    mobility, so their displacements correlate as `M_12 / M_11`. That is about `3 a / (2 r)` along the line of centres
    and `3 a / (4 r)` across it. The units are micrometres, seconds, and piconewtons, so the defaults are spheres of
    radius 1 um in water at room temperature, and each step of the system lasts `dt` seconds.

    Args:
        separation:
            The distance between the centres of the traps.

        tensor:
            `'free-draining'`, `'oseen'`, or `'rotne-prager-yamakawa'`.

        root:
            `'cholesky'` or `'lanczos'`.

        radius:
            The radius of the spheres.

        viscosity:
            The viscosity of the fluid.

        kT:
            The thermal energy.

        k:
            The stiffness of the traps.

        on_time:
            The time the traps stay on in every period.

        off_time:
            The time the traps stay off in every period.

        dt:
            The length of a step.

    Returns:
        system:
            The system, with the positions `spheres.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import hydrodynamically_coupled_traps
        >>> hydrodynamically_coupled_traps(4.0).elements[0].values['center']
        [[-2.0, 0.0], [2.0, 0.0]]
    """
    centers = [[-separation / 2, 0.0], [separation / 2, 0.0]]
    spheres = GroupElement(
        'spheres',
        element=TrappedSphereElement,
        size=2,
        values={'x': centers, 'center': centers},
    )
    traps = BlinkingTrapsInteraction(
        'tweezers',
        target='spheres',
        k=k,
        on_time=on_time,
        off_time=off_time,
        tensor=tensor,
        root=root,
        radius=radius,
        viscosity=viscosity,
        temperature=kT,
        dt=dt,
    )
    return System(elements=[spheres], interactions=[traps], space=ContinuousSpace(n_dim=2), time=DiscreteTime())

kramers_escape(n_particles=200, *, D=0.1, random_state=0)

Return particles that escape over the barrier of U = x^4 / 4 - x^2 / 2 by thermal noise of intensity D.

The particles start in the left well, and they hop between the wells at the Kramers rate sqrt(U''(-1) |U''(0)|) / (2 pi) exp(-1 / (4 D)), which is sqrt(2) / (2 pi) exp(-1 / (4 D)).

Parameters:

Name	Type	Description	Default
n_particles	int	The number of independent particles.	200
D	float	The noise intensity.	0.1
random_state	int | Generator | None	The seed or generator of the small initial spread in the left well.	0

Returns:

Name	Type	Description
system	System	The system, with the positions particles.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import kramers_escape
>>> kramers_escape(10, D=0.2).elements[0].values['D']
0.2

Source code in src/skcomplex/systems/physics/soft_matter/_diffusion.py

@register_system(
    field='physics/soft_matter',
    tags=('noise', 'sde', 'escape', 'double well', 'rates'),
    references=('Kramers 1940', 'Hänggi, Talkner, and Borkovec 1990', 'Gros 2015'),
    defaults={'until': 100.0, 'dt': 1.0, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def kramers_escape(
    n_particles: int = 200,
    *,
    D: float = 0.1,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return particles that escape over the barrier of `U = x^4 / 4 - x^2 / 2` by thermal noise of intensity `D`.

    The particles start in the left well, and they hop between the wells at the Kramers rate
    `sqrt(U''(-1) |U''(0)|) / (2 pi) exp(-1 / (4 D))`, which is `sqrt(2) / (2 pi) exp(-1 / (4 D))`.

    Args:
        n_particles:
            The number of independent particles.

        D:
            The noise intensity.

        random_state:
            The seed or generator of the small initial spread in the left well.

    Returns:
        system:
            The system, with the positions `particles.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import kramers_escape
        >>> kramers_escape(10, D=0.2).elements[0].values['D']
        0.2
    """
    positions = -1.0 + np.random.default_rng(random_state).normal(0.0, 0.1, n_particles)
    group = GroupElement(
        'particles',
        element=DoubleWellElement,
        size=n_particles,
        values={'x': positions, 'amplitude': 0.0, 'D': D},
    )
    return System(elements=[group], interactions=_get_double_well_interactions('particles'))

levy_flight(n_walkers=500, *, mu=1.2, l0=1.0, n_dim=2)

Return the Lévy flight, a walk whose jumps have lengths with the tail l^(-1 - mu) in random directions.

For mu < 2 the jump lengths have no variance, so the walk is dominated by its longest jumps and the typical distance travelled after n steps grows as n^(1 / mu), faster than the n^(1/2) of a random walk.

Parameters:

Name	Type	Description	Default
n_walkers	int	The number of independent walkers.	500
mu	float	The tail exponent of the jump lengths.	1.2
l0	float	The shortest jump length.	1.0
n_dim	int	The number of dimensions.	2

Returns:

Name	Type	Description
system	System	The system, with the positions walkers.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import levy_flight
>>> levy_flight(10, mu=1.5).elements[0].values['mu']
1.5

Source code in src/skcomplex/systems/physics/soft_matter/_anomalous.py

@register_system(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'superdiffusion', 'random walk', 'stochastic map'),
    references=('Bertin 2021', 'Shlesinger, Zaslavsky, and Frisch 1995'),
    defaults={'until': 200.0, 'dt': 1.0},
)
def levy_flight(
    n_walkers: int = 500,
    *,
    mu: float = 1.2,
    l0: float = 1.0,
    n_dim: int = 2,
) -> System:
    """Return the Lévy flight, a walk whose jumps have lengths with the tail `l^(-1 - mu)` in random directions.

    For `mu < 2` the jump lengths have no variance, so the walk is dominated by its longest jumps and the typical
    distance travelled after `n` steps grows as `n^(1 / mu)`, faster than the `n^(1/2)` of a random walk.

    Args:
        n_walkers:
            The number of independent walkers.

        mu:
            The tail exponent of the jump lengths.

        l0:
            The shortest jump length.

        n_dim:
            The number of dimensions.

    Returns:
        system:
            The system, with the positions `walkers.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import levy_flight
        >>> levy_flight(10, mu=1.5).elements[0].values['mu']
        1.5
    """
    walkers = GroupElement(
        'walkers',
        element=LevyFlierElement,
        size=n_walkers,
        values={'x': np.zeros((n_walkers, n_dim)), 'mu': mu, 'l0': l0},
    )
    flights = MapInteraction('flights', target='walkers', state='x', function=_get_flight)
    return System(elements=[walkers], interactions=[flights], space=ContinuousSpace(n_dim=n_dim))

levy_walk(n_walkers=500, *, sigma=1.5, tau0=1.0, speed=1.0, n_dim=2, random_state=0)

Return the Lévy walk, flights at constant speed in random directions for times with the tail tau^(-1 - sigma).

For 1 < sigma < 2 the flight times have a mean but no variance, so long flights dominate and the ensemble mean squared displacement grows as t^(3 - sigma), faster than Brownian motion and slower than ballistic motion.

Parameters:

Name	Type	Description	Default
n_walkers	int	The number of independent walkers.	500
sigma	float	The tail exponent of the flight times.	1.5
tau0	float	The shortest flight time.	1.0
speed	float	The speed of every flight.	1.0
n_dim	int	The number of dimensions.	2
random_state	int | Generator | None	The seed or generator of the directions of the first flights.	0

Returns:

Name	Type	Description
system	System	The system, with the positions walkers.x, the velocities walkers.v, and the turn count turns.renewals.

Examples:

>>> from skcomplex.systems.physics.soft_matter import levy_walk
>>> levy_walk(10, sigma=1.2).elements[0].values['sigma']
1.2

Source code in src/skcomplex/systems/physics/soft_matter/_anomalous.py

@register_system(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'renewal process', 'superdiffusion', 'random walk'),
    references=('Argun et al. 2021', 'Zaburdaev, Denisov, and Klafter 2015', 'Shlesinger, West, and Klafter 1987'),
    defaults={'until': 500.0, 'dt': 2.5, 'solver': RK4Solver(dt=0.5)},
)
def levy_walk(
    n_walkers: int = 500,
    *,
    sigma: float = 1.5,
    tau0: float = 1.0,
    speed: float = 1.0,
    n_dim: int = 2,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Lévy walk, flights at constant speed in random directions for times with the tail `tau^(-1 - sigma)`.

    For `1 < sigma < 2` the flight times have a mean but no variance, so long flights dominate and the ensemble mean
    squared displacement grows as `t^(3 - sigma)`, faster than Brownian motion and slower than ballistic motion.

    Args:
        n_walkers:
            The number of independent walkers.

        sigma:
            The tail exponent of the flight times.

        tau0:
            The shortest flight time.

        speed:
            The speed of every flight.

        n_dim:
            The number of dimensions.

        random_state:
            The seed or generator of the directions of the first flights.

    Returns:
        system:
            The system, with the positions `walkers.x`, the velocities `walkers.v`, and the turn count `turns.renewals`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import levy_walk
        >>> levy_walk(10, sigma=1.2).elements[0].values['sigma']
        1.2
    """
    directions = np.random.default_rng(random_state).standard_normal((n_walkers, n_dim))
    velocities = speed * directions / np.linalg.norm(directions, axis=1, keepdims=True)
    walkers = GroupElement(
        'walkers',
        element=LevyWalkerElement,
        size=n_walkers,
        values={'x': np.zeros((n_walkers, n_dim)), 'v': velocities, 'sigma': sigma, 'tau0': tau0, 'speed': speed},
    )
    flight = RateInteraction('flight', target='walkers', state='x', rate=lambda s: s.target.v)
    turns = JumpInteraction(
        'turns',
        target='walkers',
        waiting=lambda s: draw('flight', 'pareto', b=s.target.sigma, scale=s.target.tau0),
        reset={'v': lambda s: [s.target.speed * item for item in _direction('direction', len(list(s.target.v)))]},
    )
    return System(elements=[walkers], interactions=[flight, turns], space=ContinuousSpace(n_dim=n_dim))

multiplicative_noise(n_particles=500, *, s0=1.0, ds=1.8, L=100.0, convention='anti-ito')

Return free particles in a box with position-dependent noise, read in the Itô, Stratonovich, or anti-Itô sense.

Only the anti-Itô convention gives the flat stationary density of particles at thermal equilibrium without forces. The particles start at the middle of the box.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of independent particles.	500
s0	float	The noise amplitude at the middle of the box.	1.0
ds	float	The change of the noise amplitude across the box.	1.8
L	float	The length of the box.	100.0
convention	str	The convention of the noise, 'ito', 'stratonovich', or 'anti-ito'.	'anti-ito'

Returns:

Name	Type	Description
system	System	The system, with the positions particles.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import multiplicative_noise
>>> multiplicative_noise(10, convention='ito').interactions[0].convention
'ito'

Source code in src/skcomplex/systems/physics/soft_matter/_diffusion.py

@register_system(
    field='physics/soft_matter',
    tags=('noise', 'sde', 'multiplicative noise', 'spurious drift'),
    references=('Argun et al. 2021', 'Volpe and Wehr 2016'),
    defaults={'until': 100.0, 'dt': 0.05, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def multiplicative_noise(
    n_particles: int = 500,
    *,
    s0: float = 1.0,
    ds: float = 1.8,
    L: float = 100.0,
    convention: str = 'anti-ito',
) -> System:
    """Return free particles in a box with position-dependent noise, read in the Itô, Stratonovich, or anti-Itô sense.

    Only the anti-Itô convention gives the flat stationary density of particles at thermal equilibrium without
    forces. The particles start at the middle of the box.

    Args:
        n_particles:
            The number of independent particles.

        s0:
            The noise amplitude at the middle of the box.

        ds:
            The change of the noise amplitude across the box.

        L:
            The length of the box.

        convention:
            The convention of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.

    Returns:
        system:
            The system, with the positions `particles.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import multiplicative_noise
        >>> multiplicative_noise(10, convention='ito').interactions[0].convention
        'ito'
    """
    group = GroupElement(
        'particles',
        element=ModulatedParticleElement,
        size=n_particles,
        values={'x': L / 2, 's0': s0, 'ds': ds, 'L': L},
    )
    noise = RateInteraction(
        'noise',
        target='particles',
        noise=lambda s: s.target.s0 + s.target.ds * (s.target.x / s.target.L - 0.5),
        convention=convention,
    )
    return System(elements=[group], interactions=[noise])

optical_trap(n_particles=200, *, k=(1.0, 0.25), gamma=1.0, kT=1.0, amplitude=(0.0, 0.0), frequency=0.0, blink_period=None, duty=0.5)

Return overdamped particles in a harmonic optical trap, still, oscillating, or blinking.

In a still trap every coordinate is normal with variance kT / k and relaxes in the time gamma / k.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of independent particles.	200
k	Sequence[float]	The stiffness along each axis.	(1.0, 0.25)
gamma	float	The friction coefficient.	1.0
kT	float	The thermal energy.	1.0
amplitude	Sequence[float]	The amplitude of the oscillation of the center along each axis.	(0.0, 0.0)
frequency	float	The frequency of the oscillation of the center.	0.0
blink_period	float | None	The period of a blinking trap, or None for a trap that is always on.	None
duty	float	The fraction of every blinking period during which the trap is on.	0.5

Returns:

Name	Type	Description
system	System	The system, with the positions beads.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import optical_trap
>>> optical_trap(10, blink_period=2.0).interactions[0].blink_period
2.0

Source code in src/skcomplex/systems/physics/soft_matter/_diffusion.py

@register_system(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'optical tweezers', 'forcing'),
    references=('Argun et al. 2021', 'Ashkin et al. 1986', 'Jones, Maragò, and Volpe 2015'),
    defaults={'until': 20.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def optical_trap(
    n_particles: int = 200,
    *,
    k: Sequence[float] = (1.0, 0.25),
    gamma: float = 1.0,
    kT: float = 1.0,
    amplitude: Sequence[float] = (0.0, 0.0),
    frequency: float = 0.0,
    blink_period: float | None = None,
    duty: float = 0.5,
) -> System:
    """Return overdamped particles in a harmonic optical trap, still, oscillating, or blinking.

    In a still trap every coordinate is normal with variance `kT / k` and relaxes in the time `gamma / k`.

    Args:
        n_particles:
            The number of independent particles.

        k:
            The stiffness along each axis.

        gamma:
            The friction coefficient.

        kT:
            The thermal energy.

        amplitude:
            The amplitude of the oscillation of the center along each axis.

        frequency:
            The frequency of the oscillation of the center.

        blink_period:
            The period of a blinking trap, or `None` for a trap that is always on.

        duty:
            The fraction of every blinking period during which the trap is on.

    Returns:
        system:
            The system, with the positions `beads.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import optical_trap
        >>> optical_trap(10, blink_period=2.0).interactions[0].blink_period
        2.0
    """
    n_dim = len(k)
    beads = GroupElement(
        'beads',
        element=BrownianParticleElement,
        size=n_particles,
        values={'x': np.zeros((n_particles, n_dim)), 'gamma': gamma, 'temperature': kT},
    )
    trap = OpticalTrapInteraction(
        'tweezers',
        target='beads',
        k=list(k),
        amplitude=list(amplitude),
        frequency=frequency,
        blink_period=blink_period,
        duty=duty,
    )
    bath = ThermalNoiseInteraction('bath', target='beads')
    return System(elements=[beads], interactions=[trap, bath], space=ContinuousSpace(n_dim=n_dim))

rigid_dimers(n_dimers=100, *, L=1.0, v0=0.0, D_t=0.1, D_r=1.0, random_state=0)

Return colloidal molecules: pairs of Brownian or active spheres joined by rigid rods, free in the plane.

Each sphere diffuses with D_t, and an active sphere also swims at the speed v0 along its own orientation, which diffuses with D_r. The rod keeps the spheres at the distance L. Diffusing both spheres freely and then restoring their distance, the scheme of Problem 5.1 of Argun et al., gives in the limit of small steps a centre that diffuses with D_t / 2 and a rod angle that diffuses with 2 D_t / L², so the rotational diffusion of the molecule falls as L^-2 in this free-draining model. Active spheres, Problem 9.3, swim together along the mean of their orientations and turn the rod with the part of their difference across it, which raises the long-time diffusion of the centre to D_t / 2 + v0² / (4 D_r). The dimers start at the origin with random rod angles and orientations.

Parameters:

Name	Type	Description	Default
n_dimers	int	The number of dimers.	100
L	float	The length of the rod.	1.0
v0	float	The swimming speed of each sphere, 0 for Brownian spheres.	0.0
D_t	float	The translational diffusion coefficient of each sphere.	0.1
D_r	float	The rotational diffusion coefficient of the orientation of each sphere.	1.0
random_state	int | Generator | None	The seed or generator of the initial angles and orientations.	0

Returns:

Name	Type	Description
system	System	The system, with the centres dimers.c, the rod angles dimers.phi, and the orientations dimers.theta1 and dimers.theta2. The spheres sit at c ± (L / 2) (cos phi, sin phi).

Examples:

>>> from skcomplex.systems.physics.soft_matter import rigid_dimers
>>> rigid_dimers(5, L=2.0).elements[0].values['L']
2.0

Source code in src/skcomplex/systems/physics/soft_matter/_rigid_bodies.py

@register_system(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'colloidal molecule', 'rotational diffusion', 'active matter'),
    references=('Argun et al. 2021',),
    defaults={'until': 20.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.005)},
)
def rigid_dimers(
    n_dimers: int = 100,
    *,
    L: float = 1.0,
    v0: float = 0.0,
    D_t: float = 0.1,
    D_r: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return colloidal molecules: pairs of Brownian or active spheres joined by rigid rods, free in the plane.

    Each sphere diffuses with `D_t`, and an active sphere also swims at the speed `v0` along its own orientation,
    which diffuses with `D_r`. The rod keeps the spheres at the distance `L`. Diffusing both spheres freely and then
    restoring their distance, the scheme of Problem 5.1 of Argun et al., gives in the limit of small steps a centre
    that diffuses with `D_t / 2` and a rod angle that diffuses with `2 D_t / L²`, so the rotational diffusion of the
    molecule falls as `L^-2` in this free-draining model. Active spheres, Problem 9.3, swim together along the mean of
    their orientations and turn the rod with the part of their difference across it, which raises the long-time
    diffusion of the centre to `D_t / 2 + v0² / (4 D_r)`. The dimers start at the origin with random rod angles and
    orientations.

    Args:
        n_dimers:
            The number of dimers.

        L:
            The length of the rod.

        v0:
            The swimming speed of each sphere, 0 for Brownian spheres.

        D_t:
            The translational diffusion coefficient of each sphere.

        D_r:
            The rotational diffusion coefficient of the orientation of each sphere.

        random_state:
            The seed or generator of the initial angles and orientations.

    Returns:
        system:
            The system, with the centres `dimers.c`, the rod angles `dimers.phi`, and the orientations `dimers.theta1`
            and `dimers.theta2`. The spheres sit at `c ± (L / 2) (cos phi, sin phi)`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import rigid_dimers
        >>> rigid_dimers(5, L=2.0).elements[0].values['L']
        2.0
    """
    rng = np.random.default_rng(random_state)
    angles = rng.uniform(0.0, 2 * np.pi, (3, n_dimers))
    values = {
        'c': np.zeros((n_dimers, 2)),
        'phi': angles[0],
        'theta1': angles[1],
        'theta2': angles[2],
        'L': L,
        'v0': v0,
        'D_t': D_t,
        'D_r': D_r,
    }
    dimers = GroupElement('dimers', element=DimerElement, size=n_dimers, values=values)
    return System(elements=[dimers], space=ContinuousSpace(n_dim=2))

scaled_brownian_motion(n_particles=500, *, alpha=0.5, K=1.0, t0=0.01)

Return scaled Brownian motion, dx = sqrt(2 K alpha (t + t0)^(alpha - 1)) dW, with <x^2> ~ 2 K t^alpha.

The diffusion coefficient changes deterministically in time, which is the same as running Brownian motion on the rescaled clock t^alpha. The offset t0 keeps the coefficient finite at the start.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of independent walkers.	500
alpha	float	The anomalous diffusion exponent, below 1 for subdiffusion and above 1 for superdiffusion.	0.5
K	float	The generalized diffusion coefficient.	1.0
t0	float	The time offset of the diffusion coefficient.	0.01

Returns:

Name	Type	Description
system	System	The system, with the positions walkers.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import scaled_brownian_motion
>>> scaled_brownian_motion(10, alpha=1.5).elements[0].values['alpha']
1.5

Source code in src/skcomplex/systems/physics/soft_matter/_diffusion.py

@register_system(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'sde', 'non-ergodic'),
    references=('Argun et al. 2021', 'Lim and Muniandy 2002', 'Jeon et al. 2014'),
    defaults={'until': 10.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def scaled_brownian_motion(
    n_particles: int = 500,
    *,
    alpha: float = 0.5,
    K: float = 1.0,
    t0: float = 0.01,
) -> System:
    """Return scaled Brownian motion, `dx = sqrt(2 K alpha (t + t0)^(alpha - 1)) dW`, with `<x^2> ~ 2 K t^alpha`.

    The diffusion coefficient changes deterministically in time, which is the same as running Brownian motion on the
    rescaled clock `t^alpha`. The offset `t0` keeps the coefficient finite at the start.

    Args:
        n_particles:
            The number of independent walkers.

        alpha:
            The anomalous diffusion exponent, below 1 for subdiffusion and above 1 for superdiffusion.

        K:
            The generalized diffusion coefficient.

        t0:
            The time offset of the diffusion coefficient.

    Returns:
        system:
            The system, with the positions `walkers.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import scaled_brownian_motion
        >>> scaled_brownian_motion(10, alpha=1.5).elements[0].values['alpha']
        1.5
    """
    return System(
        elements=[
            GroupElement(
                'walkers',
                element=ScaledWalkerElement,
                size=n_particles,
                values={'K': K, 'alpha': alpha, 't0': t0},
            ),
        ],
        interactions=[
            RateInteraction(
                'noise',
                target='walkers',
                noise=lambda s: sqrt(2 * s.target.K * s.target.alpha * (s.t + s.target.t0) ** (s.target.alpha - 1)),
            ),
        ],
    )

shaken_spring_block(n_blocks=50, *, k=1.0, anchor=0.0, mu_s=1.0, mu_d=0.8, amplitude=2.0, shake_time=2.0, relax_time=8.0, v_stick=0.02)

Return a chain of blocks joined by springs on a frictional substrate, shaken and relaxed in turn.

Every block of unit weight feels the springs to its neighbours, k (x_{i+1} + x_{i-1} - 2 x_i), an optional spring of stiffness anchor to its rest position, and Coulomb friction. A sliding block feels the dynamic friction -mu_d sign(v), and a block at rest stays at rest while the spring force on it is below the static friction mu_s. In every cycle a random force of intensity amplitude^2 shakes the blocks for shake_time, and then the chain relaxes for relax_time to a blocked configuration, where the elongations xi_i = x_{i+1} - x_i satisfy k |xi_{i+1} - xi_i| < mu_s. Repeated cycles sample blocked configurations, whose energy per spring grows with the shaking, and whose elongations are correlated over a length that grows with that energy. Without springs between blocks and with an anchor, the blocks are independent, and their blocked positions satisfy |x_i| < mu_s / anchor. A block counts as at rest below the speed v_stick, which must exceed the change of speed (mu_s + mu_d) h of a solver step h. The ends of the chain are free, and the blocks start at rest.

Parameters:

Name	Type	Description	Default
n_blocks	int	The number of blocks.	50
k	float	The stiffness of the springs between neighbouring blocks.	1.0
anchor	float	The stiffness of the springs that hold every block to its rest position.	0.0
mu_s	float	The static friction coefficient.	1.0
mu_d	float	The dynamic friction coefficient.	0.8
amplitude	float	The amplitude of the random shaking force.	2.0
shake_time	float	The duration of the shaking in every cycle.	2.0
relax_time	float	The duration of the relaxation in every cycle.	8.0
v_stick	float	The speed below which a block counts as at rest.	0.02

Returns:

Name	Type	Description
system	System	The system, with the positions chain.x and velocities chain.v of the blocks.

Examples:

>>> from skcomplex.systems.physics.soft_matter import shaken_spring_block
>>> shaken_spring_block(5).elements[0].prepare().size_
5

Source code in src/skcomplex/systems/physics/soft_matter/_friction.py

@register_system(
    field='physics/soft_matter',
    tags=('field', 'sde', 'friction', 'granular matter', 'blocked states'),
    references=('Bertin 2021', 'Edwards and Oakeshott 1989', 'Bertin, Dauchot, and Droz 2004'),
    defaults={'until': 200.0, 'dt': 10.0, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def shaken_spring_block(
    n_blocks: int = 50,
    *,
    k: float = 1.0,
    anchor: float = 0.0,
    mu_s: float = 1.0,
    mu_d: float = 0.8,
    amplitude: float = 2.0,
    shake_time: float = 2.0,
    relax_time: float = 8.0,
    v_stick: float = 0.02,
) -> System:
    """Return a chain of blocks joined by springs on a frictional substrate, shaken and relaxed in turn.

    Every block of unit weight feels the springs to its neighbours, `k (x_{i+1} + x_{i-1} - 2 x_i)`, an optional
    spring of stiffness `anchor` to its rest position, and Coulomb friction. A sliding block feels the dynamic friction
    `-mu_d sign(v)`, and a block at rest stays at rest while the spring force on it is below the static friction
    `mu_s`. In every cycle a random force of intensity `amplitude^2` shakes the blocks for `shake_time`, and then the
    chain relaxes for `relax_time` to a blocked configuration, where the elongations `xi_i = x_{i+1} - x_i` satisfy
    `k |xi_{i+1} - xi_i| < mu_s`. Repeated cycles sample blocked configurations, whose energy per spring grows with
    the shaking, and whose elongations are correlated over a length that grows with that energy. Without springs
    between blocks and with an anchor, the blocks are independent, and their blocked positions satisfy
    `|x_i| < mu_s / anchor`. A block counts as at rest below the speed `v_stick`, which must exceed the change of
    speed `(mu_s + mu_d) h` of a solver step `h`. The ends of the chain are free, and the blocks start at rest.

    Args:
        n_blocks:
            The number of blocks.

        k:
            The stiffness of the springs between neighbouring blocks.

        anchor:
            The stiffness of the springs that hold every block to its rest position.

        mu_s:
            The static friction coefficient.

        mu_d:
            The dynamic friction coefficient.

        amplitude:
            The amplitude of the random shaking force.

        shake_time:
            The duration of the shaking in every cycle.

        relax_time:
            The duration of the relaxation in every cycle.

        v_stick:
            The speed below which a block counts as at rest.

    Returns:
        system:
            The system, with the positions `chain.x` and velocities `chain.v` of the blocks.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import shaken_spring_block
        >>> shaken_spring_block(5).elements[0].prepare().size_
        5
    """
    chain = FieldElement(
        'chain',
        element=FrictionalBlockElement,
        space=LatticeSpace(n_blocks, boundary='reflecting'),
        boundary={'x': 'zero-flux', 'v': 'zero-flux'},
        values={
            'x': np.zeros(n_blocks),
            'k': k,
            'anchor': anchor,
            'mu_s': mu_s,
            'mu_d': mu_d,
            'amplitude': amplitude,
            'shake_time': shake_time,
            'relax_time': relax_time,
            'v_stick': v_stick,
        },
    )
    return System(elements=[chain])

stochastic_resonance(x=-1.0, *, amplitude=0.1, frequency=0.01, D=0.12)

Return a particle in a double well, dx = (x - x^3 + A cos(2 pi f t)) dt + sqrt(2 D) dW.

The forcing alone is too weak to move the particle between the wells, but noise of the right intensity makes the hops follow the forcing.

Parameters:

Name	Type	Description	Default
x	float	The initial position.	-1.0
amplitude	float	The amplitude of the periodic forcing.	0.1
frequency	float	The frequency of the forcing.	0.01
D	float	The noise intensity.	0.12

Returns:

Name	Type	Description
system	System	The system, with the position particle.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import stochastic_resonance
>>> stochastic_resonance(D=0.2).elements[0].D
0.2

Source code in src/skcomplex/systems/physics/soft_matter/_brownian.py

@register_system(
    field='physics/soft_matter',
    tags=('noise', 'sde', 'forcing', 'resonance'),
    references=('Benzi, Sutera, and Vulpiani 1981', 'Gammaitoni et al. 1998'),
    defaults={'until': 500.0, 'dt': 0.5},
)
def stochastic_resonance(
    x: float = -1.0,
    *,
    amplitude: float = 0.1,
    frequency: float = 0.01,
    D: float = 0.12,
) -> System:
    """Return a particle in a double well, `dx = (x - x^3 + A cos(2 pi f t)) dt + sqrt(2 D) dW`.

    The forcing alone is too weak to move the particle between the wells, but noise of the right intensity makes the
    hops follow the forcing.

    Args:
        x:
            The initial position.

        amplitude:
            The amplitude of the periodic forcing.

        frequency:
            The frequency of the forcing.

        D:
            The noise intensity.

    Returns:
        system:
            The system, with the position `particle.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import stochastic_resonance
        >>> stochastic_resonance(D=0.2).elements[0].D
        0.2
    """
    return System(
        elements=[DoubleWellElement('particle', x=x, amplitude=amplitude, frequency=frequency, D=D)],
        interactions=_get_double_well_interactions('particle'),
    )

tilted_periodic_potential(n_particles=500, *, f=1.0, U0=0.5, L=1.0, gamma=1.0, kT=0.5)

Return overdamped particles in the tilted periodic potential V(x) = U0 cos(2 pi x / L) - f x.

The force f breaks the left and right symmetry, so the particles reach a nonequilibrium steady state with a uniform probability current J. Their mean velocity is J L = L (1 - exp(-f L / kT)) / int_0^L dx int_0^L dy exp((V(x) - V(x - y)) / kT) kT / gamma, which is f / gamma without the potential and exponentially small in the barrier for a weak force. The particles start at the origin.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of independent particles.	500
f	float	The constant force that tilts the potential.	1.0
U0	float	The amplitude of the periodic potential.	0.5
L	float	The period of the potential.	1.0
gamma	float	The friction coefficient.	1.0
kT	float	The thermal energy.	0.5

Returns:

Name	Type	Description
system	System	The system, with the positions particles.x.

Examples:

>>> from skcomplex.systems.physics.soft_matter import tilted_periodic_potential
>>> tilted_periodic_potential(10, f=2.0).interactions[0].f
2.0

Source code in src/skcomplex/systems/physics/soft_matter/_ratchets.py

@register_system(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'nonequilibrium steady state', 'current'),
    references=('Bertin 2021', 'Risken 1996', 'Reimann 2002'),
    defaults={'until': 50.0, 'dt': 0.25, 'solver': EulerMaruyamaSolver(dt=0.005)},
)
def tilted_periodic_potential(
    n_particles: int = 500,
    *,
    f: float = 1.0,
    U0: float = 0.5,
    L: float = 1.0,
    gamma: float = 1.0,
    kT: float = 0.5,
) -> System:
    """Return overdamped particles in the tilted periodic potential `V(x) = U0 cos(2 pi x / L) - f x`.

    The force `f` breaks the left and right symmetry, so the particles reach a nonequilibrium steady state with a
    uniform probability current `J`. Their mean velocity is
    `J L = L (1 - exp(-f L / kT)) / int_0^L dx int_0^L dy exp((V(x) - V(x - y)) / kT) kT / gamma`, which is
    `f / gamma` without the potential and exponentially small in the barrier for a weak force. The particles start
    at the origin.

    Args:
        n_particles:
            The number of independent particles.

        f:
            The constant force that tilts the potential.

        U0:
            The amplitude of the periodic potential.

        L:
            The period of the potential.

        gamma:
            The friction coefficient.

        kT:
            The thermal energy.

    Returns:
        system:
            The system, with the positions `particles.x`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import tilted_periodic_potential
        >>> tilted_periodic_potential(10, f=2.0).interactions[0].f
        2.0
    """
    particles = GroupElement(
        'particles',
        element=BrownianParticleElement,
        size=n_particles,
        values={'x': np.zeros((n_particles, 1)), 'gamma': gamma, 'temperature': kT},
    )
    washboard = TiltedWashboardInteraction('washboard', target='particles', U0=U0, L=L, f=f)
    bath = ThermalNoiseInteraction('bath', target='particles')
    return System(elements=[particles], interactions=[washboard, bath], space=ContinuousSpace(n_dim=1))