Hébraud–Lequeux (HL) Model — Handbook

Quick Reference

• Use when: Mean-field modeling of soft glassy materials, yield-stress fluids, foams, emulsions, pastes

• Parameters: 3 (\(\alpha\), \(\sigma_c\), \(\tau\))

• Key equation: \(\partial_t P(\sigma, t) = -\dot{\gamma}(t) \partial_\sigma P + D(t) \partial^2_\sigma P - \frac{1}{\tau} \Theta(|\sigma|-\sigma_c) P + \Gamma(t) \delta(\sigma)\)

• Test modes: Flow curve, creep, relaxation, startup, oscillation (SAOS/LAOS)

• Material examples: Foams, emulsions, pastes, concentrated colloidal suspensions, soft glassy materials

Notation Guide

Symbol	Units	Description
\(P(\sigma, t)\)	1/Pa	Probability density function of local stresses
\(\sigma\)	Pa	Local shear stress on mesoscopic block
\(\sigma_c\)	Pa	Local yield stress threshold
\(\dot{\gamma}\)	1/s	Macroscopic applied shear rate
\(\tau\)	s	Plastic relaxation time (microscopic yield time)
\(D(t)\)	Pa^2/s	Mechanical noise (stress diffusivity)
\(\Gamma(t)\)	1/s	Plastic activity rate (total yielding rate)
\(\alpha\)	—	Noise coupling parameter (control parameter)

Overview

The Hébraud–Lequeux (HL) model is a seminal mean-field elastoplastic model for soft glassy materials (SGMs), introduced by Hébraud and Lequeux in 1998 [1]. It describes the rheology of yield-stress fluids, foams, emulsions, and pastes by considering the statistical evolution of local stresses.

The model occupies a central place in the theory of amorphous materials as a minimal framework that captures the interplay between elasticity, plasticity, and mechanical noise. Unlike thermal glasses where rearrangements are driven by temperature (\(k_B T\)), in soft glasses rearrangements are often driven by stress fluctuations arising from plastic events elsewhere in the material.

Note

This implementation uses high-performance JAX kernels with a Finite Volume Method (FVM) solver, enabling efficient fitting to flow curves, creep, relaxation, and LAOS data.

Historical Context

The HL model emerged in the late 1990s alongside the Soft Glassy Rheology (SGR) model as physicists sought to understand the “jamming” transition and rheology of disordered materials.

Origins

The model was motivated by the limitations of purely phenomenological models (like Bingham or Herschel-Bulkley) which describe what happens but not why. Hébraud and Lequeux sought a microscopic justification for the flow of concentrated suspensions.

They built upon the earlier work of kinetic elastoplastic models but introduced a crucial self-consistency condition: the “mechanical noise” that facilitates local yielding is itself generated by the yielding events. This feedback loop leads to a dynamical phase transition between a solid (glassy) state and a fluid state.

Key Publications Timeline

• 1990: Kinetic models for flow in disordered media (early mean-field approaches)

• 1997: SGR model introduced by Sollich et al. (energy landscape approach)

• 1998: Hébraud & Lequeux publish “Mode-coupling theory for the pasty rheology of soft glassy materials” [1] (stress landscape approach)

• 2004: Picard et al. extend ideas to spatial models (Elastoplastic Models, EPM) [7]

• 2009: Bocquet, Colin & Ajdari derive similar equations from kinetic theory [5]

• 2018: Nicolas et al. review establishes HL as the mean-field limit of modern EPMs [10]

Relationship to SGR and STZ

The HL model sits in a landscape of microscopic theories:

• SGR (Soft Glassy Rheology): Disorder is in the energy barrier heights (trap depths). Noise is a thermal-like effective temperature \(x\). Focuses on aging and broad relaxation spectra.

• STZ (Shear Transformation Zone): Dynamics are controlled by creation/annihilation of specific structural defects (STZs). Focuses on molecular mechanisms in metallic glasses.

• HL (Hébraud–Lequeux): Disorder is in the local stress state. Noise is self-generated mechanical diffusion. Focuses on the yield stress transition and flow curves.

Mean-Field vs. Spatial Models

HL is a mean-field theory: it assumes the mechanical noise is uniform (“grey noise”). It corresponds to the infinite-range limit of spatial Elastoplastic Models (EPMs). While it misses spatial correlations (like shear bands or avalanches), it correctly predicts the macroscopic rheology and phase diagram of many soft glasses.

Physical Foundations

Mesoscopic Block Picture

The HL model discretizes the material into mesoscopic blocks of size \(\xi\) (correlation length of plastic events, typically 10-100 particle diameters). Each block is characterized by:

1. Local stress \(\sigma\): The shear stress carried by the block.

2. Yield threshold \(\sigma_c\): A critical stress above which the block becomes unstable.

3. Elastic response: Below yield, stress increases linearly with strain: \(\dot{\sigma} = G_0 \dot{\gamma}\).

Unlike EPMs which track the spatial position of every block, the HL model tracks the probability distribution \(P(\sigma, t)\) of finding a block with local stress \(\sigma\).

Mechanical Noise and Self-Consistency

The key innovation of HL is mechanical noise coupling: yielding events at one location create stress fluctuations that affect neighbors. In a mean-field approximation, these complex elastic interactions are modeled as a stochastic noise term.

The stress of a block evolves due to: 1. Macroscopic Loading: Deterministic drift \(\dot{\gamma} G_0\). 2. Mechanical Noise: Stochastic fluctuations due to neighbors yielding, modeled as a diffusive process with diffusivity \(D(t)\).

The diffusivity \(D(t)\) is self-consistently determined by the plastic activity:

\[D(t) = \alpha \Gamma(t)\]

where: - \(\Gamma(t)\) = plastic activity rate (fraction of blocks yielding per unit time) - \(\alpha\) = noise coupling parameter (material-dependent constant)

Physical interpretation: - Each yielding event redistributes stress to neighbors (via an Eshelby-like propagator). - In mean-field, the sum of many such “kicks” looks like Gaussian white noise (Central Limit Theorem). - The stronger the plastic activity \(\Gamma\), the stronger the noise \(D\).

This self-consistency creates a feedback loop: More stress \(\to\) More yielding (\(\Gamma \uparrow\)) \(\to\) More noise (\(D \uparrow\)) \(\to\) Stress spreads out (diffuses) \(\to\) More blocks reach \(\sigma_c\) \(\to\) More yielding.

Glass Transition

The parameter \(\alpha\) controls the stability of this feedback loop, leading to a phase transition:

\(\alpha\) < 0.5 (Glass phase): - Noise is insufficient to maintain a fluid state. - System “freezes” into a distribution restricted to \(|\sigma| < \sigma_c\). - Yield stress \(\sigma_y > 0\) emerges. - Aging: Relaxation timescales grow with waiting time (simple aging, \(\mu=1\)).

\(\alpha \geq 0.5\) (Fluid phase): - Noise is sufficient to maintain a steady distribution extending beyond \(\sigma_c\). - System flows at any non-zero stress. - No yield stress (\(\sigma_y = 0\)). - Finite viscosity: \(\eta \sim 1/(\alpha - 0.5)\) diverges as \(\alpha \to 0.5^+\).

Governing Equations

The evolution of the probability density function \(P(\sigma, t)\) is governed by a Fokker-Planck equation with sink and source terms.

The Master Equation

\[\partial_t P(\sigma, t) = \underbrace{-G_0 \dot{\gamma}(t) \partial_\sigma P}_{\text{Advection}} + \underbrace{D(t) \partial^2_\sigma P}_{\text{Diffusion}} - \underbrace{\frac{1}{\tau} \Theta(|\sigma|-\sigma_c) P}_{\text{Yielding}} + \underbrace{\Gamma(t) \delta(\sigma)}_{\text{Reinjection}}\]

Term-by-term breakdown:

1. Advection: \(-G_0 \dot{\gamma}(t) \partial_\sigma P\) Describes the deterministic elastic loading. All blocks move through stress space with velocity \(v = G_0 \dot{\gamma}\).

2. Diffusion: \(D(t) \partial^2_\sigma P\) Describes the mechanical noise. Stress fluctuations cause blocks to perform a random walk in stress space.

3. Yielding (Sink): \(-\frac{1}{\tau} \Theta(|\sigma|-\sigma_c) P\) Blocks with stress exceeding the threshold \(\sigma_c\) yield (relax) at a rate \(1/\tau\). The Heaviside function \(\Theta\) ensures only unstable blocks yield.

4. Reinjection (Source): \(\Gamma(t) \delta(\sigma)\) Yielded blocks dissipate their stress and are “reborn” at zero stress state \(\sigma=0\). This conserves the total number of blocks (probability).

Self-Consistency Closure

The system is closed by defining the activity rate \(\Gamma(t)\) and diffusivity \(D(t)\).

Total yielding rate (Activity):

\[\Gamma(t) = \frac{1}{\tau} \int_{|\sigma| > \sigma_c} P(\sigma, t) \, d\sigma\]

Mechanical noise (Diffusion):

\[D(t) = \alpha \Gamma(t)\]

Constraints

1. Probability Conservation: The equation preserves normalization \(\int P(\sigma, t) \, d\sigma = 1\) because the sink term (integral of yielding) exactly matches the source term (reinjection).

2. Macroscopic Stress: The measurable stress is the first moment of the distribution: \(\Sigma(t) = \int \sigma P(\sigma, t) \, d\sigma\)

Protocol-Specific Equations

The general Master Equation simplifies under specific rheological protocols.

Flow Curve (Steady Shear)

In steady state (\(\partial_t P = 0\)) under constant shear rate \(\dot{\gamma}\), the equation becomes:

\[0 = -G_0 \dot{\gamma} \partial_\sigma P + D \partial^2_\sigma P - \frac{1}{\tau}\Theta(|\sigma|-\sigma_c)P + \Gamma \delta(\sigma)\]

Analytical solutions exist for the fluid phase. For the glass phase (\(\alpha < 0.5\)), the model predicts a Herschel-Bulkley flow curve near yield:

\[\Sigma(\dot{\gamma}) \approx \Sigma_y + A \dot{\gamma}^{1/2}\]

The exponent \(n = 0.5\) is a universal prediction of the HL model (and mean-field EPMs generally) for the yielding transition.

Start-up Shear

Starting from a relaxed (aged) state, applying constant \(\dot{\gamma}\) leads to:

1. Elastic Regime: Initially, the entire distribution shifts: \(P(\sigma, t) \approx P_0(\sigma - G_0 \dot{\gamma} t)\). Stress grows linearly: \(\Sigma \sim G_0 \gamma\).

2. Overshoot: As the distribution hits \(\sigma_c\), avalanches of yielding occur. This causes a peak in stress (overshoot) before settling to steady state.

3. Steady State: Eventually, a balance between loading, yielding, and diffusion is reached.

The height and position of the overshoot depend on the initial “age” (sharpness) of the distribution.

Stress Relaxation (Step Strain)

After a rapid strain step \(\gamma_0\), the shear rate becomes zero (\(\dot{\gamma}=0\)). The equation simplifies to pure diffusion and decay:

\[\partial_t P = D(t) \partial^2_\sigma P - \frac{1}{\tau}\Theta(|\sigma|-\sigma_c)P + \Gamma(t) \delta(\sigma)\]

Crucially: - Blocks with \(|\sigma| > \sigma_c\) yield quickly (time \(\tau\)). - Blocks with \(|\sigma| < \sigma_c\) only evolve via diffusion \(D(t)\). - But \(D(t) \propto \Gamma(t)\), so as yielding stops, diffusion stops. - For \(\alpha < 0.5\), the system freezes with finite residual stress. Relaxation is incomplete.

Creep (Step Stress)

Under constant stress \(\Sigma_{app}\), the shear rate \(\dot{\gamma}(t)\) evolves to satisfy:

\[\Sigma_{app} = \int \sigma P(\sigma, t) \, d\sigma\]

This is an “inverse problem” solved numerically via a servo-controller.

• Below yield (\(\Sigma_{app} < \Sigma_y\)): \(\dot{\gamma} \to 0\). The distribution widens but remains confined.

• Above yield (\(\Sigma_{app} > \Sigma_y\)): \(\dot{\gamma} \to \text{const}\). The system flows.

• Near yield: The model exhibits delayed yielding—a long period of slow creep followed by sudden fluidization.

Oscillation (SAOS/LAOS)

For \(\gamma(t) = \gamma_0 \sin(\omega t)\):

• SAOS (Small Amplitude): The distribution oscillates slightly around 0. Linear viscoelasticity.

• LAOS (Large Amplitude): The distribution sweeps back and forth past \(\sigma_c\). - At max strain rate (zero strain), the distribution is widest. - At max strain (zero rate), the distribution is shifted. - Yielding occurs periodically, creating higher harmonics in the stress response.

Aging Analysis

The HL model captures physical aging, a hallmark of soft glassy materials.

Simple vs. Complex Aging

The HL model exhibits simple aging with an aging exponent \(\mu = 1\).

• After a quench (or cessation of flow), the system evolves towards a “frozen” state.

• The diffusion coefficient decays as \(D(t) \sim 1/t\).

• Relaxation times grow linearly with the age of the system \(t_w\).

Lifetime Distribution

The probability of a block persisting without yielding for time \(t\) scales as the age of the system. The “traps” in HL are stress states far from \(\sigma_c\). As the system ages, the distribution \(P(\sigma)\) becomes narrower and more peaked around 0, meaning blocks are “safer” from yielding.

Step Strain Response

The relaxation modulus \(G(t, t_w)\) depends on waiting time \(t_w\):

\[G(t, t_w) \approx G_{plateau} + f(t/t_w)\]

This \(t/t_w\) scaling is the signature of simple aging.

Comparison with SGR

Both HL and SGR predict \(\mu=1\) aging. - SGR: Aging is due to elements finding deeper energy traps. - HL: Aging is due to the stress distribution becoming narrower (fewer elements near yield threshold). - Difference: In HL, aging stops completely if noise vanishes (\(\alpha \to 0\)). In SGR, aging is intrinsic to the glass phase.

What You Can Learn

From fitting the HL model to experimental data, you can extract insights about the microscopic yielding processes and phase behavior.

Parameter Interpretation

\(\alpha\) (Noise Coupling Parameter)

• Definition: Ratio of mechanical noise strength to plastic activity (\(\alpha = D/\Gamma\)).

• For Graduate Students: \(\alpha\) is the control parameter for the dynamical phase transition. It measures the “fragility” of the structure. High \(\alpha\) means a single yield event triggers many others (avalanches), fluidizing the material.

• For Practitioners: Determines if the material has a yield stress. - \(\alpha < 0.5\): Paste/Gel (Yield Stress exists). - \(\alpha \ge 0.5\): Fluid (No Yield Stress). - Typical values: 0.2-0.4 for strong gels, 0.45-0.55 for soft jamming systems.

\(\sigma_c\) (Local Yield Threshold)

• Definition: The stress at which a mesoscopic block becomes unstable.

• For Graduate Students: Corresponds to the elementary energy barrier per unit volume.

• For Practitioners: Closely related to the macroscopic yield stress \(\Sigma_y\). For \(\alpha \approx 0.3\), \(\Sigma_y \approx 0.5 \sigma_c\).

\(\tau\) (Plastic Relaxation Time)

• Definition: The time it takes for a yielded block to dissipate its stress.

• For Graduate Students: The microscopic timescale of plastic rearrangement.

• For Practitioners: Sets the timescale of transient responses. If stress overshoot peaks at 1s, \(\tau\) is likely order 0.1-1s.

Material Classification Table

Parameter Range	Material Behavior	Typical Materials	Rheology Signature
\(\alpha\) < 0.3	Strong Glass	Dense granular, hard gels	High yield stress, brittle fracture
0.3 < \(\alpha\) < 0.5	Soft Glass	Foams, emulsions, pastes	Moderate yield stress, shear thinning \(n = 0.5\)
\(\alpha \approx 0.5\)	Critical Gel	Jamming transition	Power-law fluid, divergence of viscosity
\(\alpha\) > 0.5	Viscous Fluid	Dilute suspensions	Newtonian-like, no yield stress

Mathematical Summary

Measurement Protocols

Steady Shear: \(\dot{\gamma}(t) = \text{const}\) Step Strain: \(\gamma(t) = \gamma_0 \Theta(t)\) Creep: \(\Sigma(t) = \Sigma_0 \Theta(t)\) Oscillation: \(\gamma(t) = \gamma_0 \sin(\omega t)\)

Scaling Predictions

HL Model Scaling Predictions

Measurement	Regime	HL Scaling Prediction
Flow Curve	Glass (\(\alpha < 0.5\))	\(\Sigma = \Sigma_y + A \dot{\gamma}^{0.5}\) (Herschel-Bulkley \(n = 0.5\))
	Fluid (\(\alpha \to 0.5\))	\(\Sigma \sim \dot{\gamma}^{1/2}\) (Critical scaling)
	Fluid (\(\alpha > 0.5\))	\(\Sigma \sim \dot{\gamma}\) (Newtonian at low rates)
Relaxation	Glass	\(\Sigma(t) \to \Sigma_{res} > 0\) (Incomplete relaxation)
Creep	Glass (\(\Sigma < \Sigma_y\))	\(\dot{\gamma} \to 0\) (Arrest)
Viscosity	Fluid limit	\(\eta \sim (\alpha - 0.5)^{-1}\) (Divergence at transition)

Extended Features

Shear Banding

The standard HL model assumes a homogeneous shear rate. However, for certain parameter ranges (or extensions of the model), the flow curve can become non-monotonic or the material unstable, leading to shear banding—spatial coexistence of flowing and arrested regions.

Spatial Extensions (Non-local HL)

The mean-field approximation can be relaxed by introducing spatial coupling. Instead of a global \(\Gamma(t)\), one defines a local \(\Gamma(x, t)\) and allows stress diffusion to occur spatially:

\[D(x, t) = \alpha \Gamma(x, t) + l_{corr}^2 \nabla^2 \Gamma\]

This “Spatial HL” or “Kinetic Elasto-Plastic Model” connects HL to modern EPMs.

Thixotropy

While HL captures aging (slow evolution), it does not explicitly model “structure” parameters like flocculation. However, the probability distribution \(P(\sigma)\) acts as a fingerprint of the structural state. A narrow distribution corresponds to a “structured/aged” state; a wide distribution corresponds to a “rejuvenated/flowed” state.

Usage

Basic Usage

from rheojax.models import HebraudLequeux
import numpy as np

# Initialize model
model = HebraudLequeux()

# Fit to flow curve data
gamma_dot = np.logspace(-2, 1, 20)
sigma_data = np.array([...])  # Experimental stress data
model.fit(gamma_dot, sigma_data, test_mode='steady_shear')

# Predict at new shear rates
gamma_dot_pred = np.logspace(-3, 2, 50)
sigma_pred = model.predict(gamma_dot_pred, test_mode='steady_shear')

Glassy State Example

from rheojax.models import HebraudLequeux
import numpy as np

# Initialize model with glassy parameters
model = HebraudLequeux()
model.parameters.set_value("alpha", 0.3)     # Glassy phase (< 0.5)
model.parameters.set_value("sigma_c", 10.0)  # Pa
model.parameters.set_value("tau", 0.1)       # s

# Fit Flow Curve
gdot = np.logspace(-2, 1, 20)
# Load experimental data: stress_data = load_data(...)
stress_data = 10.0 + 5.0 * gdot**0.5  # Example: Herschel-Bulkley
model.fit(gdot, stress_data, test_mode='steady_shear')

# Predict Creep
time = np.linspace(0, 100, 1000)
gamma_creep = model.predict(time, test_mode='creep', sigma_applied=12.0)

Bayesian Inference

Bayesian inference is highly recommended for the HL model because the phase transition at \(\alpha = 0.5\) means small changes in \(\alpha\) can qualitatively change predictions. Quantifying uncertainty in \(\alpha\) is therefore critical.

from rheojax.models import HebraudLequeux
import numpy as np

model = HebraudLequeux()

# Generate or load data
gdot = np.logspace(-2, 1, 20)
stress_data = np.array([...])  # Experimental data

# Bayesian inference with NUTS
result = model.fit_bayesian(
    gdot, stress_data,
    test_mode='steady_shear',
    num_warmup=1000,
    num_samples=2000,
    num_chains=4
)

# Get credible intervals
intervals = model.get_credible_intervals(result.posterior_samples)
print(f"alpha: {intervals['alpha']}")
print(f"sigma_c: {intervals['sigma_c']}")

# Critical question: Is the material a glass?
alpha_samples = result.posterior_samples['alpha']
prob_glass = (alpha_samples < 0.5).mean()
print(f"Probability of glass phase: {prob_glass:.2%}")

Startup and Relaxation

from rheojax.models import HebraudLequeux
import numpy as np

model = HebraudLequeux()

# Fit to startup data (stress overshoot)
time = np.linspace(0, 10, 100)
gdot = 1.0  # Constant shear rate
stress_startup = np.array([...])  # Measured stress(t)
model.fit(time, stress_startup, test_mode='startup', gdot=gdot)

# Fit to relaxation data (step strain)
gamma0 = 0.1  # Step strain amplitude
modulus_data = np.array([...])  # G(t) = sigma(t)/gamma0
model.fit(time, modulus_data, test_mode='relaxation', gamma0=gamma0)

LAOS Analysis

Large Amplitude Oscillatory Shear reveals nonlinear material response:

from rheojax.models import HebraudLequeux
import numpy as np

model = HebraudLequeux()

# LAOS simulation
gamma0 = 1.0  # Strain amplitude
omega = 1.0   # Angular frequency (rad/s)
n_cycles = 5

time = np.linspace(0, n_cycles * 2 * np.pi / omega, 500)
stress_laos = np.array([...])  # Measured stress(t)

model.fit(time, stress_laos, test_mode='laos', gamma0=gamma0, omega=omega)

# Predict LAOS response
stress_pred = model.predict(time, test_mode='laos')

# Construct Lissajous curve
strain = gamma0 * np.sin(omega * time)
# Plot stress vs strain for Lissajous-Bowditch curve

Advanced Usage: Multiple Protocols

from rheojax.models import HebraudLequeux
import numpy as np
import matplotlib.pyplot as plt

# Initialize model
model = HebraudLequeux()

# Set to glassy state
model.parameters.set_value("alpha", 0.3)
model.parameters.set_value("sigma_c", 10.0)  # Pa
model.parameters.set_value("tau", 0.1)       # s

# 1. Fit Flow Curve
gdot = np.logspace(-2, 1, 20)
stress_data = 10.0 + 5.0 * gdot**0.5  # Example data
model.fit(gdot, stress_data, test_mode='steady_shear')

# 2. Predict multiple protocols
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Flow curve
gdot_pred = np.logspace(-3, 2, 50)
sigma_pred = model.predict(gdot_pred, test_mode='steady_shear')
axes[0, 0].loglog(gdot_pred, sigma_pred)
axes[0, 0].set_xlabel('Shear rate (1/s)')
axes[0, 0].set_ylabel('Stress (Pa)')
axes[0, 0].set_title('Flow Curve')

# Creep (above yield)
time = np.linspace(0, 100, 1000)
J_pred = model.predict(time, test_mode='creep', sigma_applied=12.0)
axes[0, 1].plot(time, J_pred)
axes[0, 1].set_xlabel('Time (s)')
axes[0, 1].set_ylabel('Compliance J(t) (1/Pa)')
axes[0, 1].set_title('Creep (σ > σ_y)')

# Relaxation
G_pred = model.predict(time, test_mode='relaxation', gamma0=0.1)
axes[1, 0].semilogy(time, G_pred)
axes[1, 0].set_xlabel('Time (s)')
axes[1, 0].set_ylabel('G(t) (Pa)')
axes[1, 0].set_title('Stress Relaxation')

# Startup
stress_startup = model.predict(time[:200], test_mode='startup', gdot=1.0)
axes[1, 1].plot(time[:200], stress_startup)
axes[1, 1].set_xlabel('Time (s)')
axes[1, 1].set_ylabel('Stress (Pa)')
axes[1, 1].set_title('Startup Flow')

plt.tight_layout()
plt.show()

Nonlinear Response

Large Amplitude Step Strain

For step strain of amplitude \(\gamma_0\) comparable to unity, the stress response depends nonlinearly on \(\gamma_0\). The initial stress distribution is shifted:

\[P(\sigma, 0^+) = P_0(\sigma - G_0 \gamma_0)\]

If \(G_0 \gamma_0 > \sigma_c\), a portion of the distribution immediately exceeds the yield threshold and yields rapidly. This causes:

1. Instantaneous stress drop as unstable blocks yield

2. Modified relaxation dynamics as the remaining distribution evolves

3. Strain softening for large amplitudes

The nonlinear relaxation modulus \(G(\gamma_0, t)\) decreases with increasing \(\gamma_0\), a phenomenon known as Type III nonlinearity.

LAOS Nonlinearity

Under large amplitude oscillatory shear \(\gamma(t) = \gamma_0 \sin(\omega t)\):

• For \(\gamma_0 \ll 1\): Linear viscoelastic response (elliptical Lissajous curves)

• For \(\gamma_0 \sim O(1)\): Significant yielding occurs each cycle

Lissajous-Bowditch Interpretation:

The parametric plot \(\sigma\) vs \(\gamma\) reveals material nonlinearity:

Lissajous Curve Shapes in HL Model

Shape	Physical Interpretation	HL Regime
Ellipse	Linear viscoelastic	\(\gamma_0 \ll 1\)
Rectangle with rounded corners	Elastic + plastic dissipation	\(\gamma_0 \sim 1\), glass phase
Parallelogram	Dominant plastic yielding	\(\gamma_0 \gg 1\)

Higher Harmonic Generation:

The stress response contains odd harmonics due to the symmetric yielding process:

\[\sigma(t) = \gamma_0 \sum_{n=1,3,5,...} \left[ G'_n \sin(n\omega t) + G''_n \cos(n\omega t) \right]\]

The third harmonic ratio \(I_{3/1} = |G^*_3|/|G^*_1|\) quantifies nonlinearity:

• \(I_{3/1} \propto \gamma_0^2\) for weak nonlinearity

• \(I_{3/1} \to O(1)\) for strong yielding

Validity and Assumptions

Model Assumptions

The HL model makes several simplifying assumptions:

1. Uniform yield threshold: All blocks have identical \(\sigma_c\)

2. Mean-field coupling: No spatial correlations between yielding events

3. Overdamped dynamics: Inertial effects are negligible

4. Athermal yielding: Thermal activation plays no role (\(k_B T \ll \sigma_c \xi^3\))

5. Instantaneous stress reset: Upon yielding, blocks immediately relax to \(\sigma=0\)

6. Affine deformation: Local strain follows macroscopic strain

Data Requirements

For robust fitting:

• Flow curve: At least 10 points spanning 2+ decades in \(\dot{\gamma}\)

• Transients: Time resolution \(\Delta t < \tau/10\)

• LAOS: At least 3-5 complete oscillation cycles

Limitations

Mean-field approximation:

Neglects spatial correlations. Real yielding events trigger neighbors (avalanches), which can lead to shear banding and strain localization not captured by the basic HL model.

Uniform threshold:

Real materials have a distribution of local yield stresses \(\rho(\sigma_c)\). Extensions exist but increase model complexity.

Sharp yield criterion:

The Heaviside \(\Theta(|\sigma| - \sigma_c)\) is a hard threshold. Softer yielding criteria can be implemented.

No microstructure tracking:

Unlike thixotropic models (DMT, MIKH), HL does not track explicit “structure” beyond the stress distribution \(P(\sigma)\).

When HL May Fail

Consider alternative models when:

• Shear banding observed: Use spatial EPM or nonlocal HL

• Pronounced thixotropy: Use DMT or MIKH models

• Temperature-sensitive: Use thermally-activated models (SGR)

• Polymer contribution: Use Maxwell + HL hybrid

• Non-monotonic flow curve: May indicate instabilities beyond HL scope

Fitting Guidance

Initialization Strategy

Step 1: Determine phase (\(\alpha < 0.5\) or \(\alpha \geq 0.5\)) - Yield stress observed: Start with \(\alpha\) = 0.3 (glass). - No yield stress: Start with \(\alpha\) = 0.7 (fluid).

Step 2: Fit flow curve - sigma_c: Constrain near observed yield stress \(\Sigma_y\). - alpha: Adjust to match curvature. Lower \(\alpha\) gives flatter (more plastic) curves.

Step 3: Fit transient data - tau: Controls relaxation speed.

Parameter Bounds

Parameter	Typical Range	Physical Constraint
alpha	0.1-1.5	\(\alpha < 0.5\) (glass), \(\alpha \geq 0.5\) (fluid)
sigma_c	\(0.5\text{--}2 \times \Sigma_y\)	Local yield threshold
tau	0.01-100 s	Relaxation timescale

Common Fitting Issues

Issue: \(\alpha\) converges to unphysical values (> 1)

Cause: Data doesn’t constrain the model, or material isn’t well-described by HL. Solution: Set tight bounds alpha.bounds = (0.1, 0.8). Check if data shows yield stress behavior.

Issue: Flow curve slope doesn’t match \(n = 0.5\)

Cause: Material may have different yielding exponent (power-law solid, polymer contribution). Solution: Consider Herschel-Bulkley with free exponent. HL universally predicts \(n = 0.5\).

Issue: Poor fit to transient data with good flow curve fit

Cause: \(\tau\) may be poorly constrained from steady-state data alone. Solution: Fit transient (startup, relaxation) data simultaneously or use Bayesian inference.

Issue: Relaxation doesn’t plateau (continues decaying)

Cause: Data may be in the fluid phase (\(\alpha\) > 0.5), or the decay timescale exceeds observation. Solution: Check phase state. Glass phase shows incomplete relaxation; fluid phase shows full decay.

Issue: Numerical oscillations in predictions

Cause: CFL violation in FVM solver (too large dt for given grid and shear rate). Solution: Reduce model.grid_n_bins or use adaptive time-stepping (automatic in RheoJAX).

Multi-Protocol Fitting

For reliable parameter estimation, fit multiple protocols simultaneously:

import numpy as np
from scipy.optimize import minimize
from rheojax.models import HebraudLequeux

def combined_objective(params, model, data_dict):
    """Objective combining multiple protocols."""
    model.parameters.set_values(params)

    residuals = []

    # Flow curve
    if 'flow_curve' in data_dict:
        gdot, stress = data_dict['flow_curve']
        stress_pred = model.predict(gdot, test_mode='steady_shear')
        residuals.append((stress - stress_pred) / stress.std())

    # Startup
    if 'startup' in data_dict:
        time, stress, gdot = data_dict['startup']
        stress_pred = model.predict(time, test_mode='startup', gdot=gdot)
        residuals.append((stress - stress_pred) / stress.std())

    return np.sum(np.concatenate(residuals)**2)

Implementation Details

The RheoJAX implementation uses an explicit Finite Volume Method (FVM) solver written in JAX.

• Advection: First-order upwind scheme (stable for hyperbolic transport).

• Diffusion: Central difference scheme.

• Time-stepping: Adaptive RK2 or Operator Splitting with Forward Euler.

• JIT Compilation: The entire solver is JIT-compiled for GPU acceleration (100x speedup).

Numerical Parameters: * sigma_max: Default 5.0 (normalized units). * n_bins: Default 501. * dt: Default 0.005.

Finite Volume Discretization

The stress axis \(\sigma \in [-\sigma_{max}, \sigma_{max}]\) is discretized into \(N\) bins:

\[\sigma_i = -\sigma_{max} + i \cdot \Delta\sigma, \quad i = 0, 1, \ldots, N-1\]

where \(\Delta\sigma = 2\sigma_{max}/(N-1)\).

The probability density is stored as cell-centered values \(P_i \approx P(\sigma_i, t)\).

Advection (Upwind Scheme):

\[\left.\frac{\partial P}{\partial t}\right|_{adv} = -v \frac{P_i - P_{i-1}}{\Delta\sigma} \quad \text{for } v > 0\]

where \(v = G_0 \dot{\gamma}\). This first-order upwind scheme is unconditionally stable for advection-dominated transport.

Diffusion (Central Difference):

\[\left.\frac{\partial P}{\partial t}\right|_{diff} = D \frac{P_{i+1} - 2P_i + P_{i-1}}{(\Delta\sigma)^2}\]

The explicit scheme requires the CFL condition:

\[\Delta t < \frac{(\Delta\sigma)^2}{2D}\]

Yielding (Sink Term):

\[\left.\frac{\partial P}{\partial t}\right|_{yield} = -\frac{1}{\tau} P_i \quad \text{for } |\sigma_i| > \sigma_c\]

Reinjection (Source Term):

The total mass lost to yielding is reinjected at \(\sigma = 0\):

\[\left.\frac{\partial P}{\partial t}\right|_{source} = \frac{\Gamma}{\Delta\sigma} \delta_{i, N/2}\]

where \(\Gamma = \tau^{-1} \sum_{|\sigma_j| > \sigma_c} P_j \Delta\sigma\).

Stability and Accuracy

The explicit Euler scheme has stability constraints:

1. Advection CFL: \(\Delta t < \Delta\sigma / |v|\)

2. Diffusion CFL: \(\Delta t < (\Delta\sigma)^2 / (2D)\)

The solver automatically sub-steps to satisfy these conditions.

For highly dynamic situations (rapid flow startup), the solver may take many sub-steps per user-specified \(\Delta t\). This is handled transparently via JAX’s lax.fori_loop.

Accuracy considerations:

• First-order upwind introduces numerical diffusion proportional to \(|v|\Delta\sigma\)

• Central difference for physical diffusion is second-order accurate

• Overall spatial accuracy: \(O(\Delta\sigma)\)

• Temporal accuracy: \(O(\Delta t)\) (forward Euler)

For high-precision work, increase n_bins (e.g., 1001) at the cost of computation time.

Protocol Kernel Details

Each rheological protocol is implemented as a specialized kernel:

Flow Curve Kernel (flow_curve_kernel):

• Runs simulation until steady state is reached

• Default: 20,000 steps at dt=0.005 (100 time units)

• Uses vmap to parallelize over multiple shear rates

Startup Kernel (startup_kernel):

• Initializes with narrow Gaussian at \(\sigma = 0\) (relaxed state)

• Applies constant shear rate and records stress history

• Returns (time, stress) arrays

Relaxation Kernel (relaxation_kernel):

• Initializes with Gaussian shifted by step strain

• Applies zero shear rate (\(\dot{\gamma} = 0\)) and records stress decay

• Returns (time, stress) arrays

Creep Kernel (creep_kernel):

• Uses servo-controller to maintain constant stress

• Adjusts shear rate dynamically: \(\dot{\gamma}_{n+1} = \dot{\gamma}_n + k_p \cdot (\Sigma_{target} - \Sigma_n) \cdot \Delta t\)

• Returns (time, strain) arrays

LAOS Kernel (laos_kernel):

• Applies oscillatory shear rate: \(\dot{\gamma}(t) = \gamma_0 \omega \cos(\omega t)\)

• Records full stress waveform

• Returns (time, stress) arrays

GPU Acceleration

The JAX implementation enables significant acceleration:

Performance Comparison

Operation	CPU Time	GPU Time (CUDA)
Single flow curve point	~50 ms	~0.5 ms
20-point flow curve (vmapped)	~1 s	~10 ms
Bayesian inference (1000 samples)	~30 min	~2 min

First call includes JIT compilation overhead (~5-30s). Subsequent calls are fast.

Theoretical Derivations

Yield Stress from Master Equation

The yield stress \(\Sigma_y\) can be derived by analyzing the steady-state master equation. In the glass phase (\(\alpha < 0.5\)), the distribution \(P(\sigma)\) is confined to \(|\sigma| < \sigma_c\) at zero shear rate.

As \(\dot{\gamma} \to 0^+\), the steady-state stress approaches:

\[\Sigma_y = \lim_{\dot{\gamma} \to 0^+} \int \sigma P(\sigma; \dot{\gamma}) \, d\sigma\]

This can be computed analytically for certain limits. The key result is that \(\Sigma_y\) exists and is positive only for \(\alpha < 0.5\).

Herschel-Bulkley Exponent Derivation

Near the yield point, the flow curve scales as:

\[\Sigma - \Sigma_y \sim \dot{\gamma}^n\]

The exponent \(n = 1/2\) arises from the balance of advection and diffusion in the wings of the distribution (near \(\sigma = \pm\sigma_c\)).

Scaling argument:

At small \(\dot{\gamma}\), the distribution has tails extending just beyond \(\sigma_c\). The width of these tails scales as:

\[\delta\sigma \sim \sqrt{D \tau} \sim \sqrt{\alpha \Gamma \tau}\]

The activity \(\Gamma\) is proportional to the tail population:

\[\Gamma \sim \frac{\delta\sigma}{\tau}\]

Self-consistently, this gives \(\delta\sigma \sim \sqrt{\alpha/\tau}\), independent of \(\dot{\gamma}\).

The excess stress \(\Sigma - \Sigma_y\) is proportional to the advective flux into the tail:

\[\Sigma - \Sigma_y \sim \dot{\gamma} \cdot (\text{time in tail}) \sim \dot{\gamma}^{1/2}\]

This universal exponent \(n = 1/2\) is a robust prediction of mean-field elastoplastic models.

Viscosity Divergence

In the fluid phase (\(\alpha > 0.5\)), the zero-shear viscosity is:

\[\eta_0 = \lim_{\dot{\gamma} \to 0} \frac{\Sigma}{\dot{\gamma}}\]

As \(\alpha \to 0.5^+\), the viscosity diverges:

\[\eta_0 \sim (\alpha - 0.5)^{-\beta}\]

with \(\beta = 1\) in the simple HL model. This divergence signals the approach to the glass transition.

References

[1] (1,2)

Hébraud, P. and Lequeux, F. “Mode-coupling theory for the pasty rheology of soft glassy materials.” Physical Review Letters, 81(14), 2934 (1998). https://doi.org/10.1103/PhysRevLett.81.2934

[2]

Fielding, S. M., Sollich, P., and Cates, M. E. “Aging and rheology in soft materials.” Journal of Rheology, 44(2), 323-369 (2000). https://doi.org/10.1122/1.551088

[3]

Sollich, P., Lequeux, F., Hébraud, P., and Cates, M. E. “Rheology of soft glassy materials.” Physical Review Letters, 78, 2020 (1997). https://doi.org/10.1103/PhysRevLett.78.2020

[4]

Coussot, P., Nguyen, Q. D., Huynh, H. T., and Bonn, D. “Viscosity bifurcation in thixotropic, yielding fluids.” Journal of Rheology, 46, 573-589 (2002). https://doi.org/10.1122/1.1459447

[5] (1,2)

Bocquet, L., Colin, A., and Ajdari, A. “Kinetic theory of plastic flow in soft glassy materials.” Physical Review Letters, 103, 036001 (2009). https://doi.org/10.1103/PhysRevLett.103.036001

[6]

Fielding, S. M., Sollich, P., & Cates, M. E. “Aging and rheology in soft materials.” Journal of Rheology, 44, 323-369 (2000). https://doi.org/10.1122/1.551088

[7]

Picard, G., Ajdari, A., Lequeux, F., & Bocquet, L. “Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids.” European Physical Journal E, 15, 371-381 (2004). https://doi.org/10.1140/epje/i2004-10054-8

[8]

Fielding, S. M. “Viscoelasticity and rheology near the soft glassy rheology transition.” Physical Review E, 76, 016311 (2007). https://doi.org/10.1039/b707980j

[9]

Bouchbinder, E. & Langer, J. S. “Nonequilibrium thermodynamics of driven amorphous materials. I. Internal degrees of freedom and volume deformation.” Physical Review E, 80, 031131 (2009). https://doi.org/10.1103/PhysRevE.80.031131

[10] (1,2)

Nicolas, A., Ferrero, E. E., Martens, K., & Barrat, J.-L. “Deformation and flow of amorphous solids: Insights from elastoplastic models.” Reviews of Modern Physics, 90, 045006 (2018). https://doi.org/10.1103/RevModPhys.90.045006

[11]

Benzi, R., Divoux, T., Barber, C., Grosso, G., & Sbragaglia, M. “Stress-dependent elastic modulus and the yieldstress fluid.” Physical Review E, 91, 023301 (2015).

[12]

Ferrero, E. E., Martens, K., & Barrat, J.-L. “Relaxation in yield stress systems through elastically interacting activated events.” Physical Review Letters, 113, 248301 (2014). https://doi.org/10.1103/PhysRevLett.113.248301

[13]

Moorcroft, R. L. & Fielding, S. M. “Criteria for shear banding in time-dependent flows of complex fluids.” Physical Review Letters, 110, 086001 (2013). https://doi.org/10.1103/PhysRevLett.110.086001

[14]

Divoux, T., Barentin, C., & Manneville, S. “From stress-induced fluidization processes to Herschel-Bulkley behaviour in simple yield stress fluids.” Soft Matter, 7, 8409-8418 (2011). https://doi.org/10.1039/c1sm05607g

[15]

Jop, P., Forterre, Y., & Pouliquen, O. “A constitutive law for dense granular flows.” Nature, 441, 727-730 (2006). https://doi.org/10.1038/nature04801

[16]

Falk, M. L. & Langer, J. S. “Dynamics of viscoplastic deformation in amorphous solids.” Physical Review E, 57, 7192 (1998). https://doi.org/10.1103/PhysRevE.57.7192

[17]

Martens, K., Bocquet, L., & Barrat, J.-L. “Connecting diffusion and dynamical heterogeneities in actively deformed amorphous systems.” Physical Review Letters, 106, 156001 (2011). https://doi.org/10.1103/PhysRevLett.106.156001

[18]

Lin, J., Lerner, E., Rosso, A., & Wyart, M. “Scaling description of the yielding transition in soft amorphous solids at zero temperature.” Proceedings of the National Academy of Sciences, 111, 14382-14387 (2014). https://doi.org/10.1073/pnas.1406391111

Further Reading

Textbooks and Reviews

• Larson, R. G. The Structure and Rheology of Complex Fluids. Oxford University Press (1999). Classic textbook covering the rheology of soft materials including yield stress fluids.

• Bonn, D., Denn, M. M., Berthier, L., Divoux, T., & Manneville, S. “Yield stress materials in soft condensed matter.” Reviews of Modern Physics, 89, 035005 (2017). Comprehensive review of yield stress materials from experimental and theoretical perspectives.

• Nicolas, A., et al. [10] The definitive review of elastoplastic models including HL. Essential reading for understanding the theoretical landscape.

Experimental Context

• Coussot, P. Rheophysics: Matter in All Its States. Springer (2014). Physical intuition for complex fluids including yield stress materials.

• Divoux, T., Barentin, C., & Manneville, S. [14] Connects HL predictions to experimental observations of creep and delayed yielding.

Theoretical Extensions

• Bocquet, Colin, & Ajdari [5] Derives HL-like equations from microscopic kinetic theory.

• Ferrero, Martens, & Barrat [12] Extends HL ideas to systems with elastically interacting activated events.

Related Models in RheoJAX

• SGR Conventional (Soft Glassy Rheology) — Handbook — Energy trap-based glass model

• SGR GENERIC (Thermodynamically Consistent) — Thermodynamically consistent SGR extension

• Shear Transformation Zone (STZ) — Defect-based plasticity model

• Elasto-Plastic Models (EPM) — Spatial elastoplastic model with avalanches

• DMT Thixotropic Models — Thixotropic models with structure kinetics

• Fluidity Models — Fluidity-based yield stress models

Glossary

Activity

The total rate of plastic yielding events, \(\Gamma(t) = \tau^{-1} \int_{|\sigma|>\sigma_c} P d\sigma\).

Coupling Parameter

The parameter \(\alpha\) controlling mechanical noise strength. Determines glass (\(< 0.5\)) vs fluid (\(\geq 0.5\)) behavior.

Elastic Loading

The deterministic drift of stress distribution due to applied shear: \(-G_0 \dot{\gamma} \partial_\sigma P\).

Glass Transition

The dynamical phase transition at \(\alpha = 0.5\) separating yield-stress solids from viscous fluids.

Mechanical Noise

Stress fluctuations from nearby yielding events, modeled as diffusion: \(D \partial^2_\sigma P\).

Mesoscopic Block

A fundamental unit of the model—a region of size \(\xi\) characterized by a single local stress \(\sigma\).

Reinjection

The source term \(\Gamma \delta(\sigma)\) returning yielded blocks to zero stress.

Self-Consistency

The closure relation \(D = \alpha \Gamma\) linking noise to activity.

Simple Aging

Aging where relaxation times scale linearly with waiting time (\(\mu = 1\)).

Stress Distribution

The probability density \(P(\sigma, t)\) of finding a block with local stress \(\sigma\).

Yield Threshold

The critical local stress \(\sigma_c\) above which blocks become unstable.

API Reference

class rheojax.models.hl.HebraudLequeux

Bases: BaseModel

Hébraud–Lequeux (HL) Model for Soft Glassy Materials.

The HL model (1998) is a mean-field description of yield-stress fluids where mesoscopic blocks of stress evolve via elastic loading, plastic yielding, and stress diffusion (mechanical noise) generated by yielding events elsewhere.

It predicts: - Finite yield stress for coupling parameter alpha < 0.5 - Herschel-Bulkley flow curves (stress ~ gdot^0.5) near yield - Creep and delayed yielding - Stress overshoots in startup flow - Non-linear LAOS response

Parameters:

• alpha – Coupling parameter (dimensionless). Controls phase state. alpha < 0.5: Glassy (yield stress) alpha >= 0.5: Fluid (no yield stress)

• tau – Microscopic yield timescale (s).

• sigma_c – Critical yield stress threshold (Pa).

parameters

ParameterSet containing alpha, tau, sigma_c.

__init__()

Initialize Hébraud–Lequeux Model.

model_function(X, params, test_mode=None, **kwargs)

Model function for Bayesian inference (NumPyro NUTS).

Parameters:

• X (ndarray) – Input array

• params (ndarray) – Parameter values [alpha, tau, sigma_c]

• test_mode (str | None) – Override test mode

• **kwargs – Protocol kwargs forwarded by BayesianMixin. Falls back to _last_fit_kwargs for values not provided.

get_phase_state()

Return the phase state based on alpha.

Return type:

str

See Also

• SGR Conventional (Soft Glassy Rheology) — Handbook — SGR model (energy-based glass transition)

• Shear Transformation Zone (STZ) — STZ model (defect-based plasticity)

• Elasto-Plastic Models (EPM) — EPM model (spatial avalanches)

• Section 3: Advanced Topics (Weeks 7-12) — Advanced elastoplastic modeling