Hebraud-Lequeux Model Overview

• Hebraud-Lequeux model is a mean-field elastoplastic framework describing yielding transitions and flow behavior in amorphous solids via self-consistent mechanical noise.
• It employs a Fokker-Planck formulation to capture stress diffusion, plastic rearrangements, and avalanche statistics under varying shear protocols.
• The model’s extensions to thermal activation, aging, and active matter reveal its sensitivity to local yielding kinetics and finite-dimensional effects.

The Hebraud-Lequeux model, also written Hébraud-Lequeux, is a mean-field elastoplastic model for the rheology and yielding transition of amorphous solids, jammed materials, and related soft disordered media. It represents a material as an ensemble of mesoscopic regions carrying local stresses, with plastic rearrangements occurring once a local threshold is exceeded and with the effect of distant rearrangements encoded as a self-consistent mechanical noise. Its defining closure is that the stress diffusion coefficient is proportional to the global plastic activity, so that noise is generated by the system’s own yielding dynamics rather than imposed externally. In this way the model provides a minimal mean-field description of athermal flow, yield-stress behavior, avalanche statistics, and several later extensions to aging, thermal activation, stochastic resetting, and active rheology (Agoritsas et al., 2015).

1. Core stochastic formulation

In its standard form, the model evolves the probability density $p(\sigma,t)$ of local shear stresses $\sigma$ according to

$$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$

with

$$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$

Here $\dot\gamma$ is the imposed shear rate, $D(t)$ is the stress diffusion coefficient, $\Gamma(t)$ is the global yielding rate, and $\delta(\sigma)$ implements stress reset after yielding. The macroscopic stress is the first moment,

$$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$

In unit-threshold conventions, yielding occurs for $|\sigma|>1$, and $\sigma$0 is the glass transition point separating fluid-like and amorphous-solid behavior (Sollich et al., 2016).

The physical interpretation of the terms is standard. The drift term $\sigma$1 describes elastic loading under applied shear, the diffusion term $\sigma$2 represents random stress kicks from remote plastic events, and the sink-source pair $\sigma$3 removes unstable sites and reinserts them at zero stress. This self-consistent diffusive closure is the characteristic HL assumption: mechanical noise is not thermal, but is generated by plastic activity elsewhere in the medium (Agoritsas et al., 2015).

An equivalent formulation uses the local distance to instability $\sigma$4, defined as the extra stress needed locally for a plastic event. In that representation, $\sigma$5 is the density of local regions about to yield. At finite temperature one writes

$\sigma$6

where the second loss term describes thermal activation for locally stable sites, and the symbol $\sigma$7 in this equation characterizes the disorder potential rather than the mechanical-noise coupling (Popović et al., 2020).

2. Stationary rheology and the yielding transition

Under steady shear, the HL closure makes the stationary problem nonlinear: one first solves the Fokker-Planck equation at fixed $\sigma$8, computes the plastic activity $\sigma$9, and then imposes $$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$0. In the standard athermal model, the low-shear behavior exhibits three regimes controlled by the coupling $$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$1. With a single normalized threshold, the critical value is $$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$2; with threshold disorder $$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$3, it generalizes to $$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$4 (Agoritsas et al., 2015).

Regime	Diffusion scaling	Macroscopic stress scaling
$$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$5	$$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$6	$$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$7
$$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$8	$$\partial_t p(\sigma,t) = -\dot\gamma\, \partial_\sigma p + D(t)\, \partial_\sigma^2 p - h(\sigma) p + \Gamma(t) \delta(\sigma),$$9	$$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$0
$$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$1	$$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$2	$$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$3

The regime $$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$4 is the yield-stress regime. It gives a Herschel-Bulkley form with exponent $$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$5, or equivalently a flow exponent $$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$6 when the flow curve is written as $$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$7. In the standard HL model this $$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$8 value is the mean-field result associated with Gaussian mechanical noise (Ferrero et al., 2019).

A notable result is that structural disorder in the local thresholds does not qualitatively modify these stationary low-shear scalings within the HL framework. It changes prefactors and replaces $$h(\sigma)=\Theta(|\sigma|-1), \qquad \Gamma(t)=\int d\sigma\, h(\sigma)\, p(\sigma,t), \qquad D(t)=\alpha \Gamma(t).$$9 by the corresponding disorder average, but it does not alter the critical structure of the three regimes or the existence of a finite macroscopic yield stress. This sharply contrasts with Soft-Glassy-Rheology-type mean-field descriptions, in which disorder is a central kinematic ingredient for generating nontrivial rheology (Agoritsas et al., 2015).

Later work generalized the local yielding rule itself. If overstressed sites yield not at a constant rate but at a progressive rate

$\dot\gamma$0

then the HL flow exponent becomes

$\dot\gamma$1

The original HL rule corresponds to $\dot\gamma$2, hence $\dot\gamma$3, while a square-root progressive rate $\dot\gamma$4, motivated by smooth effective potentials, gives $\dot\gamma$5. This establishes that the flow exponent is a dynamical quantity sensitive to microscopic yielding kinetics even within mean field (Ferrero et al., 2019).

3. Avalanche statistics and protocol dependence

The HL model is also a mean-field theory of plastic avalanches. If $\dot\gamma$6 denotes avalanche size, the distribution has the generic form

$\dot\gamma$7

A central result is that, unlike mean-field depinning, the exponent $\dot\gamma$8 depends on the dynamic protocol used to trigger avalanches (Jagla, 2015).

For random triggering, one destabilizes a randomly chosen site. The avalanche process then maps to the first-passage time distribution of a symmetric random walk with a static absorbing boundary, and one recovers

$\dot\gamma$9

the usual mean-field value also found in depinning. For quasistatic uniform loading, by contrast, all stresses are increased until the weakest site reaches instability. In this case the absorbing boundary in the random-walk mapping retracts as $D(t)$0, which enhances the survival probability of long avalanches and produces a shallower distribution with

$D(t)$1

The physical origin is the near-threshold form $D(t)$2, which differs from the depinning case where $D(t)$3 as $D(t)$4 (Jagla, 2015).

This protocol sensitivity is one of the sharpest distinctions between plastic yielding and depinning in mean field. In HL plasticity with noise, the vanishing density of near-threshold sites means that uniform loading creates correlated collective triggering, whereas in depinning uniform and random driving lead to the same $D(t)$5 exponent (Jagla, 2015).

A later mean-field treatment of ductile and brittle yielding within an HL framework further characterized the quasistatic avalanche regime. There the avalanche size distribution was written as

$D(t)$6

with $D(t)$7 and cutoff

$D(t)$8

This is consistent with the earlier conclusion that quasistatic HL avalanches are substantially shallower than the $D(t)$9 depinning benchmark and remain system-size sensitive because avalanches diverge with size throughout quasistatic loading (Parley et al., 2023).

4. Aging, memory, and linear response

In the absence of shear, the HL model displays a distinctive aging regime for $\Gamma(t)$0. The unperturbed dynamics obey

$\Gamma(t)$1

and the stress diffusion coefficient decays asymptotically as

$\Gamma(t)$2

Because $\Gamma(t)$3 is finite, only a finite amount of memory can be erased. The result is “initial condition-dependent freezing”: the interior part of the stress distribution remains predominantly unchanged, while only a shrinking boundary layer near the yielding threshold keeps evolving (Sollich et al., 2016).

The mechanical consequence is incomplete stress relaxation. For a small step strain applied at system age $\Gamma(t)$4, the shear relaxation function

$\Gamma(t)$5

does not decay to zero but approaches a plateau. Its long-time asymptotic form is

$\Gamma(t)$6

so the system becomes progressively elastic as it ages. In scaled time-difference form,

$\Gamma(t)$7

The characteristic relaxation time is proportional to age, $\Gamma(t)$8, which is the scaling expected in simple aging (Sollich et al., 2016).

The same freezing appears in frequency space. The complex modulus satisfies

$\Gamma(t)$9

so $\delta(\sigma)$0 while the dissipative part decays as

$\delta(\sigma)$1

A plausible implication is that the HL model captures a strongly arrested mean-field aging scenario in which dissipation becomes progressively suppressed, but this same feature also marks a limitation: the cited analysis notes that real soft glasses do not appear to display such extreme arrest (Sollich et al., 2016).

This limitation becomes clearer when HL is compared with mean-field models incorporating power-law mechanical noise instead of Gaussian diffusion. In that broader setting, the aging decay can be algebraic, stretched exponential, or exponential depending on the noise exponent $\delta(\sigma)$2, whereas the Gaussian HL case corresponds to the specific $\delta(\sigma)$3 decay law. This suggests that the HL aging phenomenology is tightly linked to its diffusive-noise assumption (Parley et al., 2020).

5. Thermal activation and rounding of the athermal transition

At finite temperature, the sharp athermal yielding transition at $\delta(\sigma)$4 is rounded by thermally activated flow below threshold. In the HL model, an analytical low-temperature solution can be obtained in the $\delta(\sigma)$5-representation, where $\delta(\sigma)$6 is the extra local stress needed for a plastic event and thermal activation occurs at rate $\delta(\sigma)$7 for $\delta(\sigma)$8 (Popović et al., 2020).

In the athermal limit $\delta(\sigma)$9, the steady-state distribution develops a gap $$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$0 below threshold, and the critical stress is

$$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$1

Near the transition,

$$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$2

This gap controls the activated regime: for $$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$3, flow arises from rare thermal events crossing the effective barrier set by $$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$4 (Popović et al., 2020).

The resulting low-temperature strain rate takes the scaling form

$$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$5

and, more generally,

$$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$6

For the HL model, the athermal flow exponent is $$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$7, so $$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$8. In the activated regime below threshold,

$$\langle \sigma \rangle(t)=\int d\sigma\, \sigma\, p(\sigma,t).$$9

These relations tie the rounding of the transition to both the athermal criticality and the disorder-potential exponent $|\sigma|>1$0 in the thermal activation rate (Popović et al., 2020).

The same study reported that this scaling law collapses not only HL data but also data from a 2D elasto-plastic model and previously published molecular-dynamics simulations of a 2D Lennard-Jones glass. This suggests that HL retains analytical value as a tractable reference theory for thermal rounding even when its original formulation is athermal (Popović et al., 2020).

6. Reformulations, generalizations, and limitations

Several later developments reinterpret the HL model as one member of a larger class of self-consistent stochastic processes. One reformulation casts it as diffusion with stochastic resetting, where the diffusivity depends on the average resetting rate. In the generalized model,

$|\sigma|>1$1

with the original HL case corresponding to $|\sigma|>1$2. In that resetting language, the standard HL model shows a continuous phase transition at $|\sigma|>1$3: for $|\sigma|>1$4, the small-field response is linear with finite $|\sigma|>1$5; for $|\sigma|>1$6, the diffusivity vanishes as $|\sigma|>1$7 and the response becomes nonlinear with a finite order parameter as $|\sigma|>1$8. The same extension predicts no transition for $|\sigma|>1$9 and a discontinuous phase transition for $\sigma$00 between finite-diffusivity and vanishing-diffusivity regimes (Bertin, 2022).

Another generalization adds two timescales absent from the original specification: a fluid lifetime $\sigma$01 during which a yielded element remains fluidized, and a delayed stress-propagation kernel with rate $\sigma$02. In the limit $\sigma$03 and $\sigma$04 one recovers the standard HL closure $\sigma$05, but for finite values the stationary solution can become linearly unstable, leading to spontaneous oscillations and stick-slip motion under imposed constant shear rate. The proposed mechanism is synchronization of yielding and recovery events (Bouchaud et al., 2015).

The model has also been extended to active matter. In a two-component HL variant containing passive blocks with $\sigma$06 and active blocks with $\sigma$07, both coupled through the same self-consistent diffusion rule $\sigma$08, the rheology can become non-monotonic and can display negative stresses at positive strain rates. Depending on the active fraction $\sigma$09 and the diffusion coupling $\sigma$10, one obtains passive, turnover, or fully “backwards” flow curves, with positive or negative yield stress or viscosity (Ioratim-Uba et al., 27 Aug 2025).

At the same time, several works delimit the range of validity of the original Gaussian-noise HL assumption. In a scalar yielding model coupled through an Eshelby kernel, replacing structured elastic interactions by renewed random coupling yields HL exactly; in that sense HL corresponds to a random-coupling limit with uncorrelated random-walk mechanical noise and Hurst exponent $\sigma$11 (Aguirre et al., 2018). This helps explain why some exponents are mean-field benchmarks while others are not. Subsequent work distinguished “static” exponents, such as those governing avalanche size distributions and the density of sites near threshold, from “dynamic” exponents, notably the flow exponent $\sigma$12 and dynamical exponent $\sigma$13, which change when the local yielding rate changes from constant to progressive. HL calculations display the same $\sigma$14 shift in $\sigma$15 seen in spatial elastoplastic simulations (Ferrero et al., 2019).

A stronger criticism concerns the statistics of mechanical noise itself. An improved mean-field treatment in which noise has fat tails rather than a Gaussian law predicts

$\sigma$16

with logarithmic corrections, rather than the HL value $\sigma$17. That analysis argues that the original HL assumption is not a suitable explanation for experimentally observed Herschel-Bulkley exponents, which are instead significantly affected by finite-dimensional avalanche geometry and propagation effects (Lin et al., 2017). Related generalized HL scenarios with anomalous, Lévy-type noise have been used to describe convex absorbing transitions in soft matter, where the order-parameter exponent varies continuously as $\sigma$18 for $\sigma$19 and saturates at $\sigma$20 in the Gaussian regime $\sigma$21 (Jocteur et al., 22 Jun 2026).

Taken together, these extensions position the Hebraud-Lequeux model as both a foundational mean-field theory and a controlled reference point. Its main enduring contribution is the self-consistent idea that activity generates the noise that destabilizes other regions. Its main limitation is equally clear: once the noise statistics, spatial structure, or local yielding kinetics depart from the Gaussian, memoryless assumptions of the original model, several dynamical exponents and even the character of the transition can change substantially.