TNT Loop-Bridge (Two-Species Kinetics) — Handbook

1. Quick Reference

Use When:

Telechelic polymers with two distinct chain populations:

  • Bridges: Load-bearing chains connecting two different junctions

  • Loops: Non-load-bearing chains with both ends on the same junction

Material Examples:

  • HEUR (Hydrophobically-modified Ethoxylated URethane) thickeners

  • PEG-hydrophobe associating polymers

  • Telechelic ionomers (e.g., Surlyn)

  • Flower micelle networks

  • Triblock copolymer solutions (ABA with associating end-blocks)

Key Characteristics:

  • Non-monotonic viscosity with shear rate

  • Shear thickening at moderate rates (loop-to-bridge conversion)

  • Shear thinning at high rates (bridge breakage dominates)

  • Stress overshoot in startup flows with possible initial thickening

Parameters (6):

Parameter

Default

Units

Description

\(G\)

1000

Pa

Network modulus (all bridges active)

\(\tau_b\)

1.0

s

Bridge lifetime

\(\tau_a\)

0.1

s

Loop-to-bridge association time

\(\nu\)

1.0

dimensionless

Force sensitivity (Bell exponent)

\(f_{B,eq}\)

0.5

dimensionless

Equilibrium bridge fraction

\(\eta_s\)

0.0

Pa·s

Solvent viscosity

Test Modes:

All six protocols supported:

  • FLOW_CURVE: Non-monotonic viscosity with thickening then thinning

  • STARTUP: Stress overshoot with transient thickening

  • CREEP: Population evolution under constant stress

  • RELAXATION: Bridge-to-loop conversion during stress decay

  • OSCILLATION: Linear viscoelasticity (SAOS)

  • LAOS: Nonlinear oscillatory response (population cycling)

Key Equation:

\[\frac{df_B}{dt} = \frac{1 - f_B}{\tau_a} - f_B \cdot k_{off}(\mathbf{S})\]

where \(k_{off}(\mathbf{S}) = \frac{1}{\tau_b} \exp\left[\nu \left(\text{tr}(\mathbf{S}) - 3\right)\right]\).

Only bridges contribute to stress:

\[\boldsymbol{\sigma} = G \cdot f_B \cdot (\mathbf{S} - \mathbf{I}) + 2\eta_s \mathbf{D}\]

2. Notation Guide

State Variables:

\(f_B\)

Bridge fraction (dimensionless, 0 to 1)

\(f_L\)

Loop fraction, \(f_L = 1 - f_B\)

\(\mathbf{S}\)

Bridge conformation tensor (dimensionless, 3x3 symmetric)

\(\text{tr}(\mathbf{S})\)

Trace of conformation tensor, \(S_{xx} + S_{yy} + S_{zz}\)

Material Parameters:

\(G\)

Network modulus (Pa), sets stress magnitude

\(\tau_b\)

Bridge lifetime (s), reciprocal detachment rate

\(\tau_a\)

Loop-to-bridge association time (s)

\(\nu\)

Force sensitivity (dimensionless), Bell activation exponent

\(f_{B,eq}\)

Equilibrium bridge fraction, \(\tau_b/(\tau_a + \tau_b)\)

\(\eta_s\)

Solvent viscosity (Pa·s)

Kinetic Rates:

\(k_{on}\)

Loop-to-bridge association rate, \(1/\tau_a\)

\(k_{off}(\mathbf{S})\)

Force-dependent bridge dissociation rate

\(k_{off,0}\)

Equilibrium detachment rate, \(1/\tau_b\)

Tensors:

\(\boldsymbol{\kappa}\)

Velocity gradient tensor, \(\nabla \mathbf{v}\)

\(\mathbf{D}\)

Rate-of-deformation tensor, symmetric part of \(\boldsymbol{\kappa}\)

\(\mathbf{I}\)

Identity tensor

Dimensionless Numbers:

\(\text{Wi}_b\)

Weissenberg number for bridges, \(\dot{\gamma} \tau_b\)

\(\text{Wi}_a\)

Weissenberg number for association, \(\dot{\gamma} \tau_a\)

3. Overview

Physical Picture

The TNT Loop-Bridge model describes telechelic polymers — linear chains with associating groups (stickers) at both ends. In solution, these stickers aggregate into micelles or junction zones. Each chain can exist in one of two states:

  1. Bridge: Both ends attached to different junctions, chain is stretched and carries network stress

  2. Loop: Both ends attached to the same junction, chain is relaxed and contributes no stress

At equilibrium, the bridge fraction \(f_{B,eq}\) is determined by the balance between:

  • Loop-to-bridge conversion (rate \(k_{on} = 1/\tau_a\))

  • Bridge-to-loop relaxation (rate \(k_{off,0} = 1/\tau_b\))

Under flow, this balance is perturbed:

  • Flow convection stretches bridges, increasing \(k_{off}\) (force-dependent detachment)

  • Strain can promote loop-to-bridge conversion (geometric effect)

  • Non-monotonic viscosity results: thickening at moderate rates (more bridges), thinning at high rates (bridge breakage dominates)

Historical Context

The model was introduced by:

  • Annable, Buscall, Ettelaie, and Whittlestone (1993) in Journal of Rheology as a two-species extension of the Tanaka-Edwards transient network framework

  • Built on earlier work by Tanaka & Edwards (1992) on transient networks

  • Inspired by Marrucci, Bhatt, and Ball (1993) on telechelic dynamics

  • Extended by Vaccaro & Marrucci (2000) with detailed stress-population coupling

  • Experimentally validated by Tripathi, Tam, and McKinley (2006) for HEUR solutions

The model is particularly important for understanding:

  • Shear thickening in associating polymers (e.g., paint thickeners)

  • Flow-induced structuring in telechelic networks

  • Stress overshoot and transient thickening

  • Nonlinear rheology of flower micelle networks

Key Features

Two-Species Kinetics:

Unlike single-species transient networks (e.g., basic TNT), this model tracks the bridge fraction \(f_B(t)\) as a dynamic variable. Only bridges contribute to stress, so the effective modulus is \(G_{\text{eff}} = G \cdot f_B\).

Force-Dependent Dissociation:

Bridge detachment follows Bell’s theory:

\[k_{off}(\mathbf{S}) = k_{off,0} \exp\left[\nu \left(\text{tr}(\mathbf{S}) - 3\right)\right]\]

where \(\text{tr}(\mathbf{S}) - 3\) measures chain extension (since \(\text{tr}(\mathbf{I}) = 3\)). High extension accelerates breakage.

Stress-Population Coupling:

Population kinetics and conformation evolution are fully coupled:

  • \(f_B\) controls the effective modulus

  • \(\mathbf{S}\) controls \(k_{off}\) (through force dependence)

  • Both evolve simultaneously under flow

Non-Monotonic Viscosity:

At low rates: \(\dot{\gamma} \tau_a \sim 1\) promotes loop-to-bridge conversion → viscosity increases.

At high rates: \(\dot{\gamma} \tau_b \sim 1\) and \(\text{tr}(\mathbf{S}) \gg 3\) → bridge breakage dominates → viscosity decreases.

The viscosity maximum occurs at \(\dot{\gamma}^* \sim \frac{1}{\nu \tau_b}\).

4. Physical Foundations

Telechelic Architecture

Telechelic polymers are linear chains with two functional end-groups (telechelics = “distant attractors” in Greek). Common examples:

  • HEUR: PEG backbone with hydrophobic end-caps (C16-C18 alkanes)

  • Ionomers: Polyolefin chains with ionic end-groups (e.g., Surlyn)

  • ABA triblocks: Polystyrene-PEO-Polystyrene in selective solvent

In aqueous solution, hydrophobic end-groups aggregate into micelles (typically 20-100 stickers per micelle). Each chain can bridge two micelles or form a loop within one micelle.

Bridge-Loop Equilibrium

At rest, the bridge fraction is determined by entropic and enthalpic balance:

  • Loop formation: High entropy (chain is relaxed), but only one micelle occupied

  • Bridge formation: Low entropy (chain is stretched), but connects network

The equilibrium condition \(\mu_B = \mu_L\) (equal chemical potentials) gives:

\[f_{B,eq} = \frac{\tau_b}{\tau_a + \tau_b} = \frac{k_{off,0}}{k_{on} + k_{off,0}}\]

This is independent of \(G\) (which only sets stress magnitude) and \(\nu\) (which requires extension for activation).

Typical values:

  • \(\tau_a < \tau_b\) → \(f_{B,eq} > 0.5\) (most chains are bridges)

  • \(\tau_a > \tau_b\) → \(f_{B,eq} < 0.5\) (most chains are loops)

Three-State Chain Population

The full loop-bridge framework distinguishes three chain states:

  • Bridges: Both ends attached to different junction points — these are the only chains that carry stress (elastically active)

  • Loops: Both ends attached to the same junction — carry no stress (act as dangling ends contributing viscous drag)

  • Danglers: One end free (detached) — carry no stress

The conservation constraint is:

\[f_B + f_L + f_D = 1\]

where \(f_B\), \(f_L\), \(f_D\) are the fractions of bridges, loops, and danglers respectively. In the simplified two-state model implemented in RheoJAX, danglers are absorbed into the loop population: \(f_B + f_L = 1\).

Stress-Bearing Mechanism

Only bridges contribute to stress. This is the key distinction from single-species networks.

Loops are relaxed: they exert no force on junctions. Bridges are stretched: they exert force \(\sim k_B T \cdot (\mathbf{R} - \mathbf{R}_0)\) on their end junctions.

In mean-field approximation, all bridges have the same average conformation \(\mathbf{S}\), giving stress:

\[\boldsymbol{\sigma}_{\text{bridge}} = G \cdot f_B \cdot (\mathbf{S} - \mathbf{I})\]

where:

  • \(G = n_{chain} k_B T\) (total chain density × thermal energy)

  • \(f_B\) is the fraction carrying stress

  • \((\mathbf{S} - \mathbf{I})\) is the extension beyond equilibrium

Solvent contributes Newtonian viscosity \(2\eta_s \mathbf{D}\).

Shear-Induced Association

Flow stretches bridges, which has two effects:

  1. Force-dependent detachment: High \(\text{tr}(\mathbf{S})\) increases \(k_{off}\) exponentially (Bell’s theory)

  2. Geometric promotion of bridging: Shear separates micelles, increasing the probability that a detached chain reattaches to a different micelle (not yet captured in mean-field model)

At moderate Weissenberg numbers (\(\dot{\gamma} \tau_a \sim 1\)), the second effect dominates → \(f_B\) increases above \(f_{B,eq}\) → shear thickening.

At high Weissenberg numbers (\(\dot{\gamma} \tau_b \sim 1\) and \(\text{tr}(\mathbf{S}) \gg 3\)), force-dependent detachment dominates → \(f_B\) decreases → shear thinning.

Force-Dependent Dissociation (Bell Theory)

Bond breakage under force follows Kramers theory for thermally activated escape over a barrier:

\[k_{off}(F) = k_{off,0} \exp\left[\frac{F \cdot \delta}{k_B T}\right]\]

where \(\delta\) is the activation distance. For Gaussian chains, force is proportional to extension:

\[F \sim k_B T \cdot \left(\sqrt{\text{tr}(\mathbf{S})} - \sqrt{3}\right)\]

giving:

\[k_{off}(\mathbf{S}) = k_{off,0} \exp\left[\nu \left(\text{tr}(\mathbf{S}) - 3\right)\right]\]

where \(\nu \sim \delta / a\) (activation distance over monomer size).

Typical values: \(\nu \sim 0.5\) to \(5\) (weak to strong force sensitivity).

Mean-Field Approximation

The model assumes all bridges have the same conformation \(\mathbf{S}\). This neglects:

  • Distribution of chain extensions (polydispersity effects)

  • Spatial heterogeneity (micelle clustering)

  • Correlations between bridge orientation and stress

These approximations are reasonable for:

  • Moderate extensions (\(\text{tr}(\mathbf{S}) < 10\))

  • Homogeneous flows (no spatial gradients)

  • Fast equilibration within each population

For highly stretched states or shear banding, more detailed models are needed.

5. Governing Equations

Population Kinetics

The bridge fraction evolves according to:

\[\frac{df_B}{dt} = k_{on} (1 - f_B) - k_{off}(\mathbf{S}) f_B\]

where:

  • \(k_{on} = 1/\tau_a\): Loop-to-bridge association rate (constant)

  • \(k_{off}(\mathbf{S}) = \frac{1}{\tau_b} \exp\left[\nu \left(\text{tr}(\mathbf{S}) - 3\right)\right]\): Force-dependent bridge dissociation

At equilibrium (\(df_B/dt = 0\), \(\mathbf{S} = \mathbf{I}\)):

\[f_{B,eq} = \frac{k_{on}}{k_{on} + k_{off,0}} = \frac{\tau_b}{\tau_a + \tau_b}\]

This provides a consistency check: \(f_{B,eq}\) should match the ratio of timescales.

Analytical Bridge Fraction Solution

For constant shear rate \(\dot{\gamma}_0\), the bridge fraction evolves as:

\[f_B(t) = f_B^{\text{eq}}(\dot{\gamma}_0) + \left[f_B(0) - f_B^{\text{eq}}(\dot{\gamma}_0)\right] \exp\!\left[-(k_{LB} + k_{BL})\,t\right]\]

where \(f_B^{\text{eq}}(\dot{\gamma}_0) = k_{LB} / (k_{LB} + k_{BL})\) is the equilibrium bridge fraction at shear rate \(\dot{\gamma}_0\). The approach to equilibrium occurs on the kinetic timescale \(\tau_{\text{kin}} = 1/(k_{LB} + k_{BL})\), which is independent of the chain relaxation time \(\tau_b\).

Two Characteristic Timescales

The loop-bridge model has two distinct timescales:

  1. Chain relaxation \(\tau_b\): Time for individual chain conformation to relax (sets the viscoelastic response of bridges)

  2. Bridge kinetics \(\tau_{\text{kin}} = 1/(k_{LB} + k_{BL})\): Time for the bridge fraction to equilibrate (sets the population dynamics)

When \(\tau_b \ll \tau_{\text{kin}}\), chains relax fast but the bridge fraction changes slowly — producing a two-step relaxation visible in G(t). When \(\tau_b \gg \tau_{\text{kin}}\), population dynamics are fast but chain relaxation is slow — the effective modulus adjusts quickly to \(G_{\text{eff}} = G \cdot f_B^{\text{eq}}\).

Bridge Conformation Evolution

The bridge conformation tensor evolves via upper-convected derivative with dissociation:

\[\frac{d\mathbf{S}}{dt} = \boldsymbol{\kappa} \cdot \mathbf{S} + \mathbf{S} \cdot \boldsymbol{\kappa}^T - k_{off}(\mathbf{S}) (\mathbf{S} - \mathbf{I})\]

Physical interpretation:

  • \(\boldsymbol{\kappa} \cdot \mathbf{S} + \mathbf{S} \cdot \boldsymbol{\kappa}^T\): Convective stretching (upper-convected)

  • \(-k_{off}(\mathbf{S}) (\mathbf{S} - \mathbf{I})\): Dissociation returns \(\mathbf{S} \to \mathbf{I}\) (newly attached bridges are unstretched)

Key assumption: When a bridge detaches and reattaches, it reforms with equilibrium conformation \(\mathbf{I}\) (instantaneous equilibration of loops).

Stress Tensor

Total stress is the sum of bridge elastic stress and solvent viscous stress:

\[\boldsymbol{\sigma} = G f_B (\mathbf{S} - \mathbf{I}) + 2\eta_s \mathbf{D}\]

Critical point: Only the bridge fraction contributes to elastic stress. The effective modulus is \(G_{\text{eff}}(t) = G \cdot f_B(t)\), which varies with flow history.

For simple shear (\(\kappa_{xy} = \dot{\gamma}\), \(D_{xy} = \dot{\gamma}/2\)):

\[\sigma_{xy} = G f_B S_{xy} + \eta_s \dot{\gamma}\]

Five-State ODE System (Rate-Controlled)

For steady or transient shear at imposed \(\dot{\gamma}\), the system is:

\[\begin{split}\frac{df_B}{dt} &= \frac{1 - f_B}{\tau_a} - \frac{f_B}{\tau_b} \exp\left[\nu \left(S_{xx} + S_{yy} + S_{zz} - 3\right)\right] \\ \frac{dS_{xx}}{dt} &= 2\dot{\gamma} S_{xy} - k_{off} (S_{xx} - 1) \\ \frac{dS_{yy}}{dt} &= -k_{off} (S_{yy} - 1) \\ \frac{dS_{zz}}{dt} &= -k_{off} (S_{zz} - 1) \\ \frac{dS_{xy}}{dt} &= \dot{\gamma} S_{yy} - k_{off} S_{xy}\end{split}\]

where \(k_{off} = \frac{1}{\tau_b} \exp\left[\nu \left(S_{xx} + S_{yy} + S_{zz} - 3\right)\right]\).

Initial conditions: \(f_B(0) = f_{B,eq}\), \(\mathbf{S}(0) = \mathbf{I}\).

Steady state: Solve \(df_B/dt = 0\) and \(d\mathbf{S}/dt = \mathbf{0}\) simultaneously (typically requires root-finding).

Six-State ODE System (Creep)

For imposed stress \(\sigma_0\), add strain as a state variable:

\[\begin{split}\frac{d\gamma}{dt} &= \frac{\sigma_0 - G f_B S_{xy}}{\eta_s} \\ \dot{\gamma} &= \frac{d\gamma}{dt}\end{split}\]

with the same five equations above, now coupled through \(\dot{\gamma}(t)\).

This requires implicit solution at each timestep (creep compliance).

Linearized Equations (SAOS)

For small-amplitude oscillatory shear \(\gamma(t) = \gamma_0 e^{i\omega t}\), linearize around \(f_B = f_{B,eq}\), \(\mathbf{S} = \mathbf{I}\):

\[\begin{split}f_B(t) &= f_{B,eq} + \delta f_B e^{i\omega t} \\ \mathbf{S}(t) &= \mathbf{I} + \delta \mathbf{S} e^{i\omega t}\end{split}\]

Substituting and keeping first-order terms:

\[\begin{split}i\omega \delta f_B &= -\frac{\delta f_B}{\tau_a} - \frac{\delta f_B}{\tau_b} - \frac{\nu f_{B,eq}}{\tau_b} \text{tr}(\delta \mathbf{S}) \\ i\omega \delta S_{xy} &= \dot{\gamma}_0 - \frac{\delta S_{xy}}{\tau_b}\end{split}\]

Solving gives complex modulus:

\[G^*(\omega) = G f_{B,eq} \frac{i\omega \tau_b}{1 + i\omega \tau_b} \left(1 + \frac{\nu f_{B,eq}}{1 + i\omega \tau_{\text{pop}}}\right)\]

where \(\tau_{\text{pop}} = \frac{\tau_a \tau_b}{\tau_a + \tau_b}\) is the population relaxation time.

Prediction: Two relaxation processes:

  1. Bridge conformation relaxation (time \(\tau_b\))

  2. Population redistribution (time \(\tau_{\text{pop}}\))

At high frequencies (\(\omega \tau_b \gg 1\)), \(G' \to G f_{B,eq}\) (plateau modulus).

Flow Curve (Steady Shear)

Steady-state bridge fraction and conformation satisfy:

\[\begin{split}0 &= \frac{1 - f_B^{ss}}{\tau_a} - \frac{f_B^{ss}}{\tau_b} \exp\left[\nu \left(\text{tr}(\mathbf{S}^{ss}) - 3\right)\right] \\ 0 &= \boldsymbol{\kappa} \cdot \mathbf{S}^{ss} + \mathbf{S}^{ss} \cdot \boldsymbol{\kappa}^T - k_{off}^{ss} (\mathbf{S}^{ss} - \mathbf{I})\end{split}\]

For simple shear, these are five coupled nonlinear equations in \((f_B^{ss}, S_{xx}^{ss}, S_{yy}^{ss}, S_{zz}^{ss}, S_{xy}^{ss})\).

Viscosity is:

\[\eta(\dot{\gamma}) = \frac{\sigma_{xy}^{ss}}{\dot{\gamma}} = \frac{G f_B^{ss} S_{xy}^{ss}}{\dot{\gamma}} + \eta_s\]

Non-monotonic behavior:

  • Low \(\dot{\gamma}\): \(f_B^{ss} > f_{B,eq}\) (shear thickening)

  • High \(\dot{\gamma}\): \(f_B^{ss} < f_{B,eq}\) (shear thinning from force-dependent detachment)

  • Maximum at \(\dot{\gamma}^* \sim \frac{1}{\nu \tau_b}\)

Startup Transient

For step strain rate \(\dot{\gamma}\) applied at \(t = 0\):

\[\sigma_{xy}(t) = G f_B(t) S_{xy}(t) + \eta_s \dot{\gamma}\]

Stress overshoot: Occurs when \(d\sigma_{xy}/dt = 0\) before steady state.

Transient thickening: If \(f_B(t)\) increases faster than \(S_{xy}(t)\) initially, stress can exceed the final steady-state value before the overshoot.

This is a signature of loop-to-bridge conversion under flow.

Bridge Recovery During Relaxation

After cessation of flow, the bridge fraction \(f_B\) recovers toward its quiescent equilibrium value. If flow had depleted bridges (i.e., \(f_B < f_B^{\text{eq,0}}\)), the effective modulus may increase during relaxation as bridges reform — a counter-intuitive feature where the material appears to stiffen while stress is relaxing.

This produces a characteristic non-monotonic stress relaxation: initial rapid decay (chain relaxation) followed by partial stress recovery (bridge reformation), then final decay to zero.

6. Parameter Table

Parameter

Symbol

Default

Bounds

Units

Physical Meaning

Network modulus

\(G\)

1000

(1, 1e8)

Pa

Elastic stress per unit bridge fraction

Bridge lifetime

\(\tau_b\)

1.0

(1e-6, 1e4)

s

Inverse detachment rate at equilibrium

Association time

\(\tau_a\)

0.1

(1e-6, 1e4)

s

Time for loop to become bridge

Force sensitivity

\(\nu\)

1.0

(0.01, 20)

dimensionless

Bell exponent, controls thinning onset

Equilibrium bridge fraction

\(f_{B,eq}\)

0.5

(0.01, 0.99)

dimensionless

Bridge fraction at rest

Solvent viscosity

\(\eta_s\)

0.0

(0.0, 1e4)

Pa·s

Background viscosity

Consistency Constraint:

The parameter \(f_{B,eq}\) is not independent: it should satisfy

\[f_{B,eq} = \frac{\tau_b}{\tau_a + \tau_b}\]

within experimental uncertainty. Fitting can either:

  1. Fix \(f_{B,eq}\) to this relation (5-parameter fit)

  2. Allow independent variation to test the model (6-parameter fit with consistency check)

Parameter Scaling:

For numerical stability, timescales should satisfy:

\[10^{-6} \leq \frac{\tau_a}{\tau_b} \leq 10^6\]

Extreme ratios can cause stiffness in the ODE system.

7. Parameter Interpretation

Network Modulus (\(G\))

Physical meaning: Elastic modulus if all chains were bridges (\(f_B = 1\)).

Related to chain density and temperature:

\[G = n_{\text{chain}} k_B T\]

where \(n_{\text{chain}}\) is the number density of polymer chains.

Fitting: Determined from the plateau modulus in SAOS:

\[G_0' = G f_{B,eq}\]

If \(f_{B,eq}\) is known independently, \(G\) is directly obtained.

Typical values:

  • Dilute HEUR (1 wt%): \(G \sim 10\) to \(100\) Pa

  • Semi-dilute (5 wt%): \(G \sim 100\) to \(1000\) Pa

  • Concentrated (>10 wt%): \(G \sim 1000\) to \(10000\) Pa

Bridge Lifetime (\(\tau_b\))

Physical meaning: Average time a bridge remains attached before detaching (at equilibrium, \(\mathbf{S} = \mathbf{I}\)).

Controlled by:

  • Sticker-micelle interaction energy \(\Delta G_{\text{bind}}\) (higher → longer \(\tau_b\))

  • Temperature (higher → shorter \(\tau_b\), Arrhenius)

  • Micelle aggregation number (larger → more stable)

Fitting: Terminal relaxation time from SAOS:

\[\tau_b \approx \frac{1}{\omega_c}\]

where \(\omega_c\) is the crossover frequency (\(G' = G''\)).

Typical values:

  • Fast dynamics (low molecular weight): \(\tau_b \sim 0.01\) to \(0.1\) s

  • Moderate: \(\tau_b \sim 0.1\) to \(10\) s

  • Slow (strong association): \(\tau_b \sim 10\) to \(1000\) s

Association Time (\(\tau_a\))

Physical meaning: Time for a loop to convert into a bridge.

Controlled by:

  • Sticker diffusion rate (how fast a free end finds a new micelle)

  • Micelle density (higher → shorter \(\tau_a\))

  • Chain flexibility (stiffer → longer \(\tau_a\))

Fitting: Determines the rate of shear thickening. Compare transient startup at different \(\dot{\gamma}\):

  • If thickening occurs at \(t \sim \tau_a\), association is being probed

  • Faster thickening → shorter \(\tau_a\)

Also affects the width of the thickening regime in the flow curve.

Typical values:

  • Fast reassociation: \(\tau_a \sim 0.001\) to \(0.01\) s

  • Moderate: \(\tau_a \sim 0.01\) to \(1\) s

  • Slow (diffusion-limited): \(\tau_a \sim 1\) to \(100\) s

Ratio interpretation:

  • \(\tau_a / \tau_b < 1\) → Most chains are bridges at rest (\(f_{B,eq} > 0.5\))

  • \(\tau_a / \tau_b > 1\) → Most chains are loops at rest (\(f_{B,eq} < 0.5\)), large thickening potential

Force Sensitivity (\(\nu\))

Physical meaning: Strength of force-dependent dissociation (Bell exponent).

Related to the activation distance \(\delta\) over which the sticker must move to escape the micelle:

\[\nu \sim \frac{\delta}{a}\]

where \(a\) is the monomer size.

Fitting: Controls the onset of shear thinning:

  • High \(\nu\) → steep increase in \(k_{off}\) with extension → early thinning

  • Low \(\nu\) → weak force dependence → extended thickening regime

Experimental signature: Slope of \(\log(\eta)\) vs \(\log(\dot{\gamma})\) in the thinning regime:

\[\frac{d \log \eta}{d \log \dot{\gamma}} \propto -\nu\]

Typical values:

  • Weak force sensitivity: \(\nu \sim 0.1\) to \(0.5\)

  • Moderate: \(\nu \sim 0.5\) to \(2\)

  • Strong: \(\nu \sim 2\) to \(10\)

Equilibrium Bridge Fraction (\(f_{B,eq}\))

Physical meaning: Fraction of chains in bridge state at rest.

Determined by free energy balance:

\[f_{B,eq} = \frac{\tau_b}{\tau_a + \tau_b}\]

Independent measurement: Estimate from plateau modulus if \(G\) is known:

\[f_{B,eq} = \frac{G_0'}{G}\]

Fitting strategy:

  1. Fit \(\tau_a\), \(\tau_b\) from rheology

  2. Compute \(f_{B,eq}\) from ratio

  3. Check consistency with \(G_0' / G\)

Discrepancies indicate:

  • Model inadequacy (e.g., dangling ends, super-bridges)

  • Parameter correlation (need more experimental constraints)

Typical values:

  • Loop-dominated: \(f_{B,eq} \sim 0.1\) to \(0.4\)

  • Balanced: \(f_{B,eq} \sim 0.4\) to \(0.6\)

  • Bridge-dominated: \(f_{B,eq} \sim 0.6\) to \(0.9\)

Solvent Viscosity (\(\eta_s\))

Physical meaning: Background viscosity from solvent and unassociated chains.

For aqueous HEUR: \(\eta_s \approx \eta_{\text{water}} \sim 0.001\) Pa·s (often negligible).

For concentrated solutions: \(\eta_s\) includes polymer overlap contributions.

Fitting: High-shear limiting viscosity:

\[\eta_{\infty} = \eta_s\]

as \(f_B \to 0\) and \(\mathbf{S} \to \mathbf{I}\) at \(\dot{\gamma} \to \infty\).

Typical values:

  • Dilute aqueous: \(\eta_s \sim 0.001\) to \(0.01\) Pa·s

  • Semi-dilute: \(\eta_s \sim 0.01\) to \(1\) Pa·s

  • Polymer melt: \(\eta_s \sim 1\) to \(1000\) Pa·s

8. Validity and Assumptions

The TNT Loop-Bridge model makes the following simplifications:

1. Two-Species Approximation

Assumes only two states: bridges and loops. Neglects:

  • Dangling ends (one sticker attached, one free)

  • Free chains (both stickers detached)

  • Super-bridges (one chain bridging multiple micelles)

  • Multi-loop configurations

Valid when: Sticker-micelle binding energy is high (\(\Delta G_{\text{bind}} \gg k_B T\)), so free ends are rare.

2. Mean-Field Conformation

All bridges have the same average conformation \(\mathbf{S}\). Neglects:

  • Distribution of chain extensions

  • Spatial heterogeneity

  • Correlation between orientation and local stress

Valid when: Moderate extensions (\(\text{tr}(\mathbf{S}) < 10\)), homogeneous flows, fast equilibration.

3. Instantaneous Loop Reformation

When a bridge detaches, the loop immediately equilibrates to \(\mathbf{S} = \mathbf{I}\). Neglects:

  • Memory of previous conformation

  • Transient stretching of loops

Valid when: Loop relaxation time \(\ll \tau_b\) (loops are much faster than bridges).

4. Affine Deformation

Chain conformation follows the macroscopic velocity gradient (upper-convected derivative). Neglects:

  • Chain slip

  • Non-affine motion near junctions

  • Brownian fluctuations

Valid when: \(\text{Wi} < 10\) to \(100\) (moderate Weissenberg number).

5. No Hydrodynamic Interactions

Micelles are treated as fixed junctions; solvent flow around micelles is ignored. Neglects:

  • Micelle mobility

  • Cooperative rearrangements

Valid when: Micelle size \(\ll\) chain contour length, dilute to semi-dilute regime.

6. Homogeneous Flow

Assumes no spatial gradients (0D model). Cannot describe:

  • Shear banding

  • Flow instabilities

  • Spatial heterogeneity

Valid when: Gap width \(\ll\) characteristic length scale of structure evolution.

7. Bell-Type Dissociation

Force-dependent detachment follows exponential form. Neglects:

  • Non-Gaussian chain statistics (finite extensibility)

  • Stick-slip dynamics

  • Multiple activation pathways

Valid when: Extensions are moderate (\(\text{tr}(\mathbf{S}) < 10\) to \(20\)).

When to Use This Model:

  • HEUR solutions in dilute to semi-dilute regime

  • Telechelic ionomers with well-defined end-groups

  • Flower micelle networks with two-state kinetics

  • Shear thickening followed by thinning in flow curves

  • Stress overshoot in startup flows

When NOT to Use:

  • High extensions (\(\text{tr}(\mathbf{S}) > 20\)) → need finite extensibility (FENE-type)

  • Shear banding → need spatial (1D) model

  • Multiple relaxation processes → need more species or modes

  • Chain entanglements dominate → need reptation-based model

9. Regimes and Behavior

Linear Viscoelastic Regime (\(\text{Wi} \ll 1\))

Behavior: Single Maxwell-like relaxation with effective modulus.

Modulus:

\[G_0' = G f_{B,eq}\]

Relaxation time: \(\tau_b\) (bridge detachment controls terminal relaxation).

Population: \(f_B \approx f_{B,eq}\) (small perturbations).

Key prediction: Storage modulus plateau at high frequency reflects the equilibrium bridge fraction.

Shear Thickening Regime (\(\text{Wi}_a \sim 1\), \(\text{Wi}_b < 1\))

Conditions: \(\dot{\gamma} \tau_a \sim 1\), extensions are still moderate.

Mechanism: Flow promotes loop-to-bridge conversion faster than force-dependent detachment.

Bridge fraction: \(f_B > f_{B,eq}\) (excess bridges).

Viscosity: Increases above zero-shear value:

\[\eta(\dot{\gamma}) > \eta_0 = \frac{G f_{B,eq} \tau_b}{1 + \eta_s / (G f_{B,eq} \tau_b)}\]

Experimental signature: Positive slope in \(\log \eta\) vs \(\log \dot{\gamma}\) plot.

Material examples: HEUR at 5-10 wt%, moderate shear rates (10-100 1/s).

Critical Shear Rate (\(\dot{\gamma}^*\))

Definition: Shear rate at which viscosity is maximum.

Scaling estimate:

\[\dot{\gamma}^* \sim \frac{1}{\nu \tau_b}\]

Physical picture: Balance between:

  • Thickening (loop-to-bridge conversion, rate \(\sim 1/\tau_a\))

  • Thinning (force-dependent detachment, rate \(\sim (1/\tau_b) e^{\nu (\text{tr}(\mathbf{S}) - 3)}\))

At \(\dot{\gamma}^*\), extensions are \(\text{tr}(\mathbf{S}) \sim 5\) to \(10\).

Parameter dependence:

  • Large \(\nu\) → low \(\dot{\gamma}^*\) (early onset of thinning)

  • Large \(\tau_b\) → low \(\dot{\gamma}^*\) (longer-lived bridges are more sensitive to force)

Shear Thinning Regime (\(\text{Wi}_b \gg 1\))

Conditions: \(\dot{\gamma} \tau_b \gg 1\), high extensions \(\text{tr}(\mathbf{S}) \gg 3\).

Mechanism: Force-dependent detachment dominates, \(k_{off} \gg k_{on}\).

Bridge fraction: \(f_B < f_{B,eq}\) (depletion of bridges).

Viscosity: Power-law decrease:

\[\eta(\dot{\gamma}) \sim \dot{\gamma}^{-p}\]

where \(p \sim 0.5\) to \(1\) (depends on \(\nu\)).

Experimental signature: Negative slope in \(\log \eta\) vs \(\log \dot{\gamma}\), stress plateau or weak increase.

Material examples: HEUR at high shear (>1000 1/s), approaching \(\eta_s\).

Startup Transient Behavior

Stress overshoot: Occurs when \(d\sigma/dt = 0\) before steady state.

Overshoot time: \(t_{\text{peak}} \sim \tau_b\) (bridge relaxation).

Overshoot magnitude: Depends on \(\text{Wi}_a\):

  • Low \(\text{Wi}_a\) (\(< 0.1\)): Weak overshoot, \(\sigma_{\text{peak}} \approx 1.1 \sigma_{ss}\)

  • Moderate \(\text{Wi}_a\) (\(\sim 1\)): Strong overshoot with transient thickening, \(\sigma_{\text{peak}} > 1.5 \sigma_{ss}\)

  • High \(\text{Wi}_a\) (\(> 10\)): Overshoot returns to steady-state thinning

Transient thickening: If \(f_B(t)\) increases faster than \(S_{xy}(t)\) grows, stress can exceed the final value before relaxing to steady state.

This is a signature of loop-to-bridge conversion and distinguishes telechelic networks from simple Maxwell models.

Creep Compliance

Step stress: \(\sigma_0\) applied at \(t = 0\).

Short-time behavior: Elastic response \(\gamma \sim \sigma_0 / (G f_{B,eq})\).

Intermediate time: Population redistribution, \(f_B\) adjusts to new stress.

Long-time behavior: Viscous flow \(\gamma \sim t / \eta_s\) if \(\eta_s > 0\).

Compliance:

\[J(t) = \frac{\gamma(t)}{\sigma_0}\]

Shows double-exponential relaxation (bridge relaxation + population adjustment) followed by linear flow.

Creep Rupture Via Bridge Collapse

Under sustained applied stress, the bridge fraction can progressively decrease:

  1. Stress stretches bridge chains, increasing their breakage rate

  2. Reduced \(f_B\) decreases the effective modulus \(G_{\text{eff}} = G \cdot f_B\)

  3. Lower modulus means higher strain for the same stress → more chain stretch

  4. Positive feedback: above a critical stress, \(f_B \to 0\) (all chains become loops)

This bridge collapse mechanism produces creep rupture — delayed catastrophic failure under sustained load. The critical stress depends on the ratio \(k_{BL}/k_{LB}\) and the force sensitivity of the bridge-to-loop transition.

Oscillatory Shear (SAOS and LAOS)

SAOS (small amplitude): Two relaxation processes:

  1. Bridge conformation (time \(\tau_b\))

  2. Population redistribution (time \(\tau_{\text{pop}} = \tau_a \tau_b / (\tau_a + \tau_b)\))

Cole-Cole plot: \(G''\) vs \(G'\) shows deviation from single Maxwell (two arcs).

LAOS (large amplitude): Population cycling — \(f_B(t)\) oscillates due to periodic stretching/relaxation.

Higher harmonics appear when \(k_{off}(t)\) becomes nonlinear (high strain amplitudes).

10. What You Can Learn from Fitting This Model

Equilibrium Bridge Fraction

From: Plateau modulus in SAOS, \(G_0'\).

Relation:

\[f_{B,eq} = \frac{G_0'}{G}\]

Interpretation: Fraction of chains actively contributing to network stress at rest.

Application: Design of thickeners (maximize \(f_{B,eq}\) for high viscosity at rest).

Bridge Lifetime and Association Time

From: Terminal relaxation frequency \(\omega_c\) (where \(G' = G''\)) and flow curve thickening regime.

Relations:

\[\begin{split}\tau_b &\sim \frac{1}{\omega_c} \\ \tau_a &\sim \frac{1}{\dot{\gamma}_{\text{thick}}}\end{split}\]

where \(\dot{\gamma}_{\text{thick}}\) is the shear rate where thickening begins.

Interpretation:

  • \(\tau_b\): Sticker-micelle binding strength

  • \(\tau_a\): Chain diffusion and reassociation rate

Application: Tuning molecular architecture (e.g., end-group hydrophobicity) to control timescales.

Force Sensitivity (Bell Exponent)

From: Slope of thinning regime in flow curve, onset of thinning \(\dot{\gamma}^*\).

Relation:

\[\nu \sim \frac{\log(k_{off,\text{high}} / k_{off,0})}{\text{tr}(\mathbf{S}_{\text{high}}) - 3}\]

Interpretation: Sensitivity of detachment to chain extension (activation barrier shape).

Application: Predicting failure under high-rate deformation (e.g., extrusion, spraying).

Shear Thickening Amplitude

From: Maximum viscosity \(\eta_{\text{max}}\) compared to zero-shear \(\eta_0\).

Relation:

\[\frac{\eta_{\text{max}}}{\eta_0} \propto \frac{f_{B,\text{max}}}{f_{B,eq}}\]

where \(f_{B,\text{max}}\) is the peak bridge fraction.

Interpretation: Amount of excess bridging achievable under flow.

Dependence:

  • Large \(\tau_a / \tau_b\) ratio → large thickening (more loops available for conversion)

  • Small \(\nu\) → sustained thickening (weak force dependence)

Application: Optimizing thickening efficiency in coatings, adhesives.

Onset of Shear Thinning

From: Critical shear rate \(\dot{\gamma}^*\) where \(\eta(\dot{\gamma})\) is maximum.

Relation:

\[\dot{\gamma}^* \sim \frac{1}{\nu \tau_b}\]

Interpretation: Transition from association-dominated to detachment-dominated regime.

Application: Processing windows for molding, extrusion (avoid thinning).

Stress Overshoot and Transient Thickening

From: Startup experiments at various \(\dot{\gamma}\).

Metrics:

  • Overshoot magnitude: \(\sigma_{\text{peak}} / \sigma_{ss}\)

  • Overshoot time: \(t_{\text{peak}}\)

  • Transient thickening: Does \(\sigma(t) > \sigma_{ss}\) occur?

Interpretation:

  • Large overshoot → significant population redistribution

  • Transient thickening → loop-to-bridge conversion is fast relative to conformation evolution

Application: Predicting behavior in rapid deformations (printing, dispensing).

Network Connectivity

From: Comparison of \(f_{B,eq}\) (from \(G_0'\)) with \(\tau_a / (\tau_a + \tau_b)\) (from timescales).

Consistency check: Should match within 10-20%.

Discrepancies indicate:

  • Presence of dangling ends (reduces \(G_0'\))

  • Multi-species populations (more than loops/bridges)

  • Spatial heterogeneity (micelle clustering)

Application: Validating model assumptions, guiding molecular design.

11. Experimental Design and Data Requirements

Small-Amplitude Oscillatory Shear (SAOS)

Objective: Measure \(G'(\omega)\) and \(G''(\omega)\) to extract:

  • Plateau modulus \(G_0'\) → \(f_{B,eq}\) (if \(G\) known)

  • Terminal relaxation time \(\tau_b\)

  • Population relaxation time \(\tau_{\text{pop}}\)

Protocol:

  1. Strain sweep to identify linear regime (typically \(\gamma_0 < 0.1\))

  2. Frequency sweep from \(\omega = 10^{-2}\) to \(10^2\) rad/s

  3. Check for low-frequency plateau (terminal behavior) and high-frequency plateau (\(G_0'\))

Expected features:

  • Two relaxation processes (double-Maxwellian)

  • \(G'\) plateau at high \(\omega\) → \(G_0' = G f_{B,eq}\)

  • Crossover at \(\omega_c \sim 1/\tau_b\)

Data quality:

  • At least 10 points per decade in \(\omega\)

  • Signal-to-noise ratio >100:1 for \(G'\), \(G''\)

  • Temperature control \(\pm 0.1\) K (association is thermally activated)

Steady Shear Flow Curve

Objective: Measure \(\sigma(\dot{\gamma})\) or \(\eta(\dot{\gamma})\) to observe:

  • Shear thickening at moderate \(\dot{\gamma}\)

  • Viscosity maximum at \(\dot{\gamma}^*\)

  • Shear thinning at high \(\dot{\gamma}\)

Protocol:

  1. Shear rate sweep from \(\dot{\gamma} = 10^{-3}\) to \(10^3\) 1/s (log-spaced)

  2. Allow steady state at each point (wait \(> 5 \tau_b\))

  3. Check reversibility (up-down sweep) to detect thixotropy

Expected features:

  • Zero-shear plateau \(\eta_0\) at low \(\dot{\gamma}\)

  • Viscosity peak at \(\dot{\gamma}^*\)

  • Power-law thinning at high \(\dot{\gamma}\)

  • High-shear plateau \(\eta_\infty = \eta_s\)

Data quality:

  • At least 15-20 points covering the full curve

  • Steady-state criterion: \(|\dot{\sigma}/\sigma| < 0.01\) over time \(\sim \tau_b\)

  • Check for edge fracture, slip at high \(\dot{\gamma}\)

Startup Shear Flow

Objective: Measure transient \(\sigma(t)\) after step \(\dot{\gamma}\) to extract:

  • Stress overshoot magnitude and time

  • Transient thickening (if \(\sigma(t) > \sigma_{ss}\) before overshoot)

  • Approach to steady state

Protocol:

  1. Pre-shear at low \(\dot{\gamma}\) to equilibrate

  2. Rest for time \(> 10 \tau_b\) (return to \(f_B = f_{B,eq}\))

  3. Step to target \(\dot{\gamma}\), measure \(\sigma(t)\) for \(t = 0\) to \(10 \tau_b\)

  4. Repeat for multiple \(\dot{\gamma}\) spanning thickening and thinning regimes

Expected features:

  • Overshoot at \(t \sim \tau_b\)

  • Transient thickening for \(\dot{\gamma} \sim 1/\tau_a\)

  • Steady-state approach as \(\sigma(t) \to \sigma_{ss}\)

Data quality:

  • Sampling rate \(> 10 / \tau_b\) (resolve overshoot peak)

  • Repeat 3-5 times to check reproducibility

  • Monitor strain: \(\gamma(t) = \int \dot{\gamma} dt\) should be linear

Creep Compliance

Objective: Measure \(\gamma(t)\) after step stress \(\sigma_0\) to extract:

  • Initial compliance \(J_0 = 1/(G f_{B,eq})\)

  • Retardation times \(\tau_b\), \(\tau_{\text{pop}}\)

  • Long-time viscous flow \(\gamma \sim t/\eta_s\)

Protocol:

  1. Apply constant stress \(\sigma_0\) (within linear regime if possible)

  2. Measure \(\gamma(t)\) for \(t = 0\) to \(100 \tau_b\)

  3. Repeat for multiple \(\sigma_0\) to check linearity

Expected features:

  • Instantaneous compliance jump \(J_0\)

  • Retardation (slow relaxation) at intermediate times

  • Linear flow \(J(t) \sim t/\eta_s\) if \(\eta_s > 0\)

Data quality:

  • High-resolution strain measurement (resolution \(< 0.001\))

  • Check for instrument compliance (subtract out)

  • Temperature drift \(< 0.1\) K over measurement

Large-Amplitude Oscillatory Shear (LAOS)

Objective: Measure higher harmonics in \(\sigma(t)\) during oscillation to probe:

  • Population cycling \(f_B(t)\)

  • Nonlinear force-dependent detachment

  • Intracycle shear thickening-thinning

Protocol:

  1. Apply sinusoidal strain \(\gamma(t) = \gamma_0 \sin(\omega t)\)

  2. Vary \(\gamma_0\) from 0.1 (linear) to 5 (highly nonlinear)

  3. Measure stress waveform \(\sigma(t)\) at high resolution (>100 points per cycle)

  4. Fourier decompose to extract harmonics \(G_1'(\omega), G_1''(\omega), G_3'(\omega), G_3''(\omega), \ldots\)

Expected features:

  • Odd harmonics dominate (material symmetry)

  • Third harmonic \(G_3'\) increases with \(\gamma_0\) (nonlinear elastic response)

  • Intracycle thickening-thinning visible in Lissajous plots

Data quality:

  • Waveform sampling >200 points/cycle

  • Multiple cycles (10+) to check periodicity

  • Fourier truncation error <1%

Temperature Control

Association kinetics are thermally activated:

\[\tau_b(T) \sim \exp\left[\frac{\Delta G_{\text{bind}}}{k_B T}\right]\]

Requirement: Temperature control \(\pm 0.1\) K to avoid drift in \(\tau_a\), \(\tau_b\).

Protocol: Equilibrate sample for 10-15 minutes at target temperature before testing.

12. Computational Implementation

State Variables

For rate-controlled flows (imposed \(\dot{\gamma}\)):

\[\mathbf{y} = [f_B, S_{xx}, S_{yy}, S_{zz}, S_{xy}]^T\]

For stress-controlled flows (creep):

\[\mathbf{y} = [f_B, S_{xx}, S_{yy}, S_{zz}, S_{xy}, \gamma]^T\]

ODE Right-Hand Side (Rate-Controlled)

Implemented in rheojax.models.tnt._kernels.tnt_loop_bridge_ode_rhs:

def tnt_loop_bridge_ode_rhs(y, t, gamma_dot, params):
    """
    y = [f_B, S_xx, S_yy, S_zz, S_xy]
    params = [G, tau_b, tau_a, nu, f_B_eq, eta_s]
    """
    f_B, S_xx, S_yy, S_zz, S_xy = y
    G, tau_b, tau_a, nu, f_B_eq, eta_s = params

    # Trace of conformation tensor
    tr_S = S_xx + S_yy + S_zz

    # Force-dependent detachment rate
    k_off = (1.0 / tau_b) * jnp.exp(nu * (tr_S - 3.0))

    # Population kinetics
    df_B_dt = (1.0 - f_B) / tau_a - f_B * k_off

    # Conformation evolution (upper-convected)
    dS_xx_dt = 2.0 * gamma_dot * S_xy - k_off * (S_xx - 1.0)
    dS_yy_dt = -k_off * (S_yy - 1.0)
    dS_zz_dt = -k_off * (S_zz - 1.0)
    dS_xy_dt = gamma_dot * S_yy - k_off * S_xy

    return jnp.array([df_B_dt, dS_xx_dt, dS_yy_dt, dS_zz_dt, dS_xy_dt])

JAX compilation: Use jax.jit for 10-50x speedup:

ode_rhs_jit = jax.jit(tnt_loop_bridge_ode_rhs, static_argnums=(2,))

ODE Right-Hand Side (Creep)

For imposed stress \(\sigma_0\), add implicit equation for \(\dot{\gamma}(t)\):

def tnt_loop_bridge_creep_rhs(y, t, sigma_0, params):
    """
    y = [f_B, S_xx, S_yy, S_zz, S_xy, gamma]
    """
    f_B, S_xx, S_yy, S_zz, S_xy, gamma = y
    G, tau_b, tau_a, nu, f_B_eq, eta_s = params

    # Current stress (implicit equation)
    sigma_xy = G * f_B * S_xy + eta_s * gamma_dot

    # Solve for gamma_dot: sigma_xy = sigma_0
    gamma_dot = (sigma_0 - G * f_B * S_xy) / eta_s

    # Rest is same as rate-controlled
    tr_S = S_xx + S_yy + S_zz
    k_off = (1.0 / tau_b) * jnp.exp(nu * (tr_S - 3.0))

    df_B_dt = (1.0 - f_B) / tau_a - f_B * k_off
    dS_xx_dt = 2.0 * gamma_dot * S_xy - k_off * (S_xx - 1.0)
    dS_yy_dt = -k_off * (S_yy - 1.0)
    dS_zz_dt = -k_off * (S_zz - 1.0)
    dS_xy_dt = gamma_dot * S_yy - k_off * S_xy
    d_gamma_dt = gamma_dot

    return jnp.array([df_B_dt, dS_xx_dt, dS_yy_dt, dS_zz_dt, dS_xy_dt, d_gamma_dt])

Initial Conditions

For all transient protocols, start from equilibrium:

y0 = jnp.array([f_B_eq, 1.0, 1.0, 1.0, 0.0])  # [f_B, S_xx, S_yy, S_zz, S_xy]

For creep, add \(\gamma(0) = 0\):

y0_creep = jnp.array([f_B_eq, 1.0, 1.0, 1.0, 0.0, 0.0])

Steady-State Root-Finding (Flow Curve)

For flow curve, solve \(dy/dt = 0\) using Newton-Raphson:

from jax.scipy.optimize import root_scalar

def residual(y_ss, gamma_dot, params):
    dy_dt = tnt_loop_bridge_ode_rhs(y_ss, 0.0, gamma_dot, params)
    return dy_dt

# Initial guess: equilibrium
y_guess = jnp.array([f_B_eq, 1.0, 1.0, 1.0, 0.0])

# Solve for each gamma_dot
y_ss = root_scalar(residual, y_guess, args=(gamma_dot, params))

Bracketing: For non-monotonic viscosity, root-finding may fail. Use continuation from low to high \(\dot{\gamma}\), using previous solution as initial guess.

SAOS Linearization

For small-amplitude oscillatory shear, linearize around equilibrium:

\[G^*(\omega) = G f_{B,eq} \frac{i\omega \tau_b}{1 + i\omega \tau_b} \left(1 + \frac{\nu f_{B,eq}}{1 + i\omega \tau_{\text{pop}}}\right)\]

where \(\tau_{\text{pop}} = \tau_a \tau_b / (\tau_a + \tau_b)\).

Implemented analytically (no ODE solve needed):

def saos_modulus(omega, params):
    G, tau_b, tau_a, nu, f_B_eq, eta_s = params
    tau_pop = (tau_a * tau_b) / (tau_a + tau_b)

    # Bridge relaxation
    G_bridge = G * f_B_eq * (1j * omega * tau_b) / (1 + 1j * omega * tau_b)

    # Population correction
    G_pop_corr = (nu * f_B_eq) / (1 + 1j * omega * tau_pop)

    G_star = G_bridge * (1 + G_pop_corr)
    return G_star.real, G_star.imag  # (G', G'')

LAOS Simulation

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

# Time-dependent strain rate
def gamma_dot_laos(t, gamma_0, omega):
    return gamma_0 * omega * jnp.cos(omega * t)

# Solve ODE with time-varying gamma_dot
def ode_rhs_laos(y, t, gamma_0, omega, params):
    gamma_dot = gamma_dot_laos(t, gamma_0, omega)
    return tnt_loop_bridge_ode_rhs(y, t, gamma_dot, params)

# Integrate for 10 cycles
t_end = 10 * (2 * jnp.pi / omega)
sol = odeint(ode_rhs_laos, y0, t_span, args=(gamma_0, omega, params))

Extract harmonics via FFT after discarding transient cycles (first 5).

Numerical Stability

Stiffness: When \(\tau_a \ll \tau_b\) or \(\nu\) is large, the system becomes stiff (fast and slow timescales).

Solution: Use implicit solver (e.g., jax-scipy.integrate.solve_ivp with method='BDF').

Adaptive timestepping: Essential for capturing overshoot peak and steady-state approach.

Tolerances: Set rtol=1e-6, atol=1e-8 for accurate stress predictions.

Vectorization: Use jax.vmap to solve for multiple \(\dot{\gamma}\) in parallel:

solve_startup_vmap = jax.vmap(solve_startup, in_axes=(None, 0, None))
# Vectorize over gamma_dot array
sigma_array = solve_startup_vmap(y0, gamma_dot_array, params)

13. Fitting Guidance and Best Practices

Parameter Initialization

Step 1: Estimate \(G\) and \(f_{B,eq}\) from SAOS:

Plateau modulus at high frequency:

\[G_0' \approx G f_{B,eq}\]

If \(f_{B,eq}\) is unknown, assume \(f_{B,eq} \approx 0.5\) (balanced) as initial guess.

Step 2: Estimate \(\tau_b\) from terminal relaxation:

Crossover frequency \(\omega_c\) where \(G' = G''\):

\[\tau_b \approx \frac{1}{\omega_c}\]

Step 3: Estimate \(\tau_a\) from consistency:

If \(f_{B,eq}\) is known:

\[\tau_a = \tau_b \left(\frac{1 - f_{B,eq}}{f_{B,eq}}\right)\]

Step 4: Estimate \(\nu\) from flow curve slope:

Onset of thinning \(\dot{\gamma}^*\):

\[\nu \sim \frac{1}{\tau_b \dot{\gamma}^*}\]

Step 5: Set \(\eta_s = 0\) initially (negligible solvent viscosity).

Fitting Strategy (Hierarchical)

Stage 1: SAOS only (fix \(\nu\), \(\eta_s = 0\))

Fit \(G\), \(\tau_b\), \(\tau_a\) (or \(f_{B,eq}\)) to \(G'(\omega)\), \(G''(\omega)\).

Objective: Minimize

\[\chi^2 = \sum_i \left[\frac{G'_{\text{pred}}(\omega_i) - G'_{\text{data}}(\omega_i)}{\sigma_{G'}}\right]^2 + \sum_i \left[\frac{G''_{\text{pred}}(\omega_i) - G''_{\text{data}}(\omega_i)}{\sigma_{G''}}\right]^2\]

Stage 2: Add flow curve (release \(\nu\))

Fix \(G\), \(\tau_b\), \(\tau_a\) from Stage 1. Fit \(\nu\) to \(\eta(\dot{\gamma})\).

Objective: Minimize

\[\chi^2 = \sum_j \left[\frac{\eta_{\text{pred}}(\dot{\gamma}_j) - \eta_{\text{data}}(\dot{\gamma}_j)}{\sigma_{\eta}}\right]^2\]

Stage 3: Global refinement (all parameters)

Use Stage 1-2 results as initial guess. Fit all parameters simultaneously to SAOS + flow curve.

Stage 4: Add transient data (validate)

Fix all parameters from Stage 3. Predict startup \(\sigma(t)\) and compare to data (no fitting).

If discrepancies exist, refine \(\tau_a\) or \(\nu\) slightly.

Parameter Bounds

Always enforce physical bounds:

Parameter

Lower Bound

Upper Bound

\(G\)

1 Pa

\(10^8\) Pa

\(\tau_b\)

\(10^{-6}\) s

\(10^4\) s

\(\tau_a\)

\(10^{-6}\) s

\(10^4\) s

\(\nu\)

0.01

20

\(f_{B,eq}\)

0.01

0.99

\(\eta_s\)

0 Pa·s

\(10^4\) Pa·s

Weighting of Data

Log-space weighting: For quantities spanning decades (e.g., \(\eta\), \(\omega\)), use log-space residuals:

\[\chi^2 = \sum_i \left[\log(\eta_{\text{pred}}) - \log(\eta_{\text{data}})\right]^2\]

This prevents high-viscosity points from dominating the fit.

Uncertainty weighting: If experimental uncertainties \(\sigma_i\) are known, use inverse-variance weighting:

\[\chi^2 = \sum_i \frac{(y_{\text{pred},i} - y_{\text{data},i})^2}{\sigma_i^2}\]

Consistency Checks

1. Equilibrium bridge fraction:

Compare \(f_{B,eq}\) from fit with \(\tau_b / (\tau_a + \tau_b)\). Should agree within 10%.

2. Plateau modulus:

Check \(G f_{B,eq} \approx G_0'\) from SAOS.

3. Timescale separation:

Verify \(\tau_{\text{pop}} = \tau_a \tau_b / (\tau_a + \tau_b) < \tau_b\).

4. Critical shear rate:

Check \(\dot{\gamma}^* \sim 1 / (\nu \tau_b)\) matches observed viscosity maximum.

Common Pitfalls

1. Parameter correlation:

\(G\) and \(f_{B,eq}\) are highly correlated in SAOS (only product \(G f_{B,eq}\) is constrained). Need flow curve or independent \(f_{B,eq}\) measurement.

2. Overfitting transient data:

Startup curves can be fit well with incorrect parameters (e.g., too high \(\nu\) compensated by too low \(\tau_a\)). Always validate against multiple protocols.

3. Ignoring solvent viscosity:

If high-shear viscosity does not plateau at \(\eta_s \approx 0\), \(\eta_s > 0\) is required. Neglecting it biases \(G\), \(\tau_b\) estimates.

4. Non-unique solutions:

Non-monotonic viscosity can have multiple parameter sets giving similar fits. Use physical constraints (e.g., \(f_{B,eq} < 0.9\) for loop-dominated systems).

15. Usage Examples

Example 1: Basic SAOS Prediction

from rheojax.models.tnt import TNTLoopBridge
import jax.numpy as jnp
import matplotlib.pyplot as plt

# Create model with default parameters
model = TNTLoopBridge()

# SAOS prediction
omega = jnp.logspace(-2, 2, 100)
G_prime, G_double_prime = model.predict(omega, test_mode='oscillation')

# Plot
plt.loglog(omega, G_prime, 'o-', label="G'")
plt.loglog(omega, G_double_prime, 's-', label="G''")
plt.xlabel('Frequency (rad/s)')
plt.ylabel('Modulus (Pa)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

Example 2: Flow Curve with Non-Monotonic Viscosity

# Predict flow curve
gamma_dot = jnp.logspace(-3, 3, 50)
viscosity = model.predict(gamma_dot, test_mode='flow_curve')

# Plot
plt.loglog(gamma_dot, viscosity, 'o-')
plt.xlabel('Shear rate (1/s)')
plt.ylabel('Viscosity (Pa·s)')
plt.title('Non-Monotonic Viscosity')
plt.grid(True, alpha=0.3)
plt.show()

Example 3: Startup with Transient Thickening

# Startup at moderate shear rate
t = jnp.linspace(0, 20, 500)
gamma_dot = 1.0  # Wi_a ~ 0.1

stress = model.predict(t, test_mode='startup', gamma_dot=gamma_dot)

# Plot
plt.plot(t, stress, '-')
plt.xlabel('Time (s)')
plt.ylabel('Stress (Pa)')
plt.title(f'Startup at gamma_dot = {gamma_dot} 1/s')
plt.grid(True, alpha=0.3)
plt.show()

Example 4: Fitting to Experimental Data

from rheojax.core.data import RheoData

# Load experimental SAOS data
omega_exp = ...  # rad/s
G_star_exp = ...  # Pa (complex)

# Fit model
rheo_data = RheoData(x=omega_exp, y=G_star_exp, test_mode='oscillation')
result = model.fit(rheo_data)

# Print fitted parameters
params = model.get_parameter_values()
print(f"G = {params['G']:.1f} Pa")
print(f"tau_b = {params['tau_b']:.3f} s")
print(f"tau_a = {params['tau_a']:.3f} s")
print(f"nu = {params['nu']:.2f}")
print(f"f_B_eq = {params['f_B_eq']:.2f}")

Example 5: Bayesian Inference

# Bayesian inference with NumPyro
result_bayes = model.fit_bayesian(
    rheo_data,
    num_warmup=1000,
    num_samples=2000,
    num_chains=4,
    seed=42
)

# Check convergence
print(f"R-hat: {result_bayes.diagnostics['r_hat']}")
print(f"ESS: {result_bayes.diagnostics['ess']}")

# Credible intervals
intervals = model.get_credible_intervals(
    result_bayes.posterior_samples,
    credibility=0.95
)
print(f"tau_a: [{intervals['tau_a'][0]:.3f}, {intervals['tau_a'][1]:.3f}] s")

Example 6: Parameter Sensitivity

# Vary nu to see effect on flow curve
nu_values = [0.5, 1.0, 2.0, 5.0]
gamma_dot = jnp.logspace(-2, 3, 50)

for nu in nu_values:
    model.set_parameter('nu', nu)
    eta = model.predict(gamma_dot, test_mode='flow_curve')
    plt.loglog(gamma_dot, eta, label=f'nu = {nu}')

plt.xlabel('Shear rate (1/s)')
plt.ylabel('Viscosity (Pa·s)')
plt.legend()
plt.title('Effect of Force Sensitivity')
plt.grid(True, alpha=0.3)
plt.show()

14. Failure Mode: Loop Saturation

Under extreme flow conditions, all bridge chains convert to loops:

\[f_B \to 0, \quad f_L \to 1\]

In this limit, the network loses all elasticity (\(G_{\text{eff}} \to 0\)) and behaves as a viscous fluid. This loop saturation represents the complete destruction of the stress-bearing network.

Physical signatures:

  • Flow curve shows a dramatic viscosity drop at high rates

  • Startup stress overshoot followed by near-complete stress collapse

  • Recovery after flow cessation governed by \(\tau_{\text{kin}}\) (bridge reformation time)

16. See Also

TNT Shared Reference:

TNT Base Model:

Generalizations:

Alternative Models for Similar Materials:

17. API Reference

class rheojax.models.tnt.TNTLoopBridge[source]

Bases: TNTBase

Loop-bridge kinetics model for telechelic polymer networks.

Implements reversible loop-bridge kinetics for telechelic polymers where chains can exist as load-bearing bridges (both ends on different junctions) or non-load-bearing loops (both ends on same junction).

The bridge fraction f_B evolves dynamically via attachment (loops → bridges) and force-activated dissociation (bridges → loops via Bell kinetics).

State Variables

  • f_B: Bridge fraction (0 to 1)

  • S: Conformation tensor of bridges [S_xx, S_yy, S_zz, S_xy]

Key Equations

Bridge fraction kinetics:

df_B/dt = (1 - f_B)/τ_a - f_B·β(stretch)
β(stretch) = (1/τ_b)·exp(ν·(stretch - 1))

Conformation tensor (bridges only):

dS/dt = L·S + S·L^T + g_0·I - β(stretch)·S

Stress (only bridges carry load):

σ = f_B·G·S_xy + η_s·γ̇
param G:

Network modulus (fully bridged, Pa), bounds (1e0, 1e8)

type G:

float, default 1e3

param tau_b:

Bridge detachment time (s), bounds (1e-6, 1e4)

type tau_b:

float, default 1.0

param tau_a:

Loop attachment time (s), bounds (1e-6, 1e4)

type tau_a:

float, default 5.0

param nu:

Bell force sensitivity (dimensionless), bounds (0.01, 20)

type nu:

float, default 1.0

param f_B_eq:

Equilibrium bridge fraction, bounds (0.01, 0.99)

type f_B_eq:

float, default 0.5

param eta_s:

Solvent viscosity (Pa·s), bounds (0.0, 1e4)

type eta_s:

float, default 0.0

parameters

Model parameters

Type:

ParameterSet

fitted_

Whether the model has been fitted

Type:

bool

Examples

Basic usage:

>>> model = TNTLoopBridge()
>>> gamma_dot = np.logspace(-2, 2, 50)
>>> sigma = model.predict(gamma_dot, test_mode='flow_curve')

Startup with bridge fraction tracking:

>>> t = np.linspace(0, 100, 500)
>>> stress, f_B = model.simulate_startup(
...     t, gamma_dot=10.0, return_bridge_fraction=True
... )

See also

TNTSingleMode

Basic single-mode TNT (constant breakage)

TNTCates

Living polymer (wormlike micelle) model

__init__()[source]

Initialize TNT loop-bridge model.

initialize_from_creep(t, strain, sigma_applied=1.0)[source]

Initialize parameters from creep data for numerical stability.

Key insight: For creep simulation, eta_s must be non-zero to prevent infinite initial shear rate when S_xy starts at 0.

Parameters:

  • t (ndarray) – Time array (s)

  • strain (ndarray) – Strain array

  • sigma_applied (float) – Applied constant stress (Pa)

Return type:

None

initialize_from_relaxation(t, modulus)[source]

Initialize parameters from stress relaxation data.

Uses conservative tau_b estimate to ensure numerical stability with typical pre-shear rates (Wi = gamma_dot * tau_b should be < ~100).

Parameters:

  • t (ndarray) – Time array (s)

  • modulus (ndarray) – Relaxation modulus G(t) (Pa)

Return type:

None

property G: float

Get network modulus G (Pa).

property tau_b: float

Get bridge detachment time τ_b (s).

property tau_a: float

Get loop attachment time τ_a (s).

property nu: float

Get Bell force sensitivity ν (dimensionless).

property f_B_eq: float

Get equilibrium bridge fraction f_B_eq (dimensionless).

property eta_s: float

Get solvent viscosity η_s (Pa·s).

property G_eff: float

Get effective modulus G_eff = f_B_eq·G (Pa).

This is the linearized modulus at equilibrium.

property eta_0: float

Get zero-shear viscosity η₀ = f_B_eq·G·τ_b + η_s (Pa·s).

get_equilibrium_state()[source]

Return equilibrium state [f_B_eq, 1, 1, 1, 0].

At rest: f_B = f_B_eq, S = I (unstretched, isotropic)

Returns:

Equilibrium state [f_B, S_xx, S_yy, S_zz, S_xy]

Return type:

Array

model_function(X, params, test_mode=None, **kwargs)[source]

NumPyro/BayesianMixin model function.

Routes to appropriate prediction based on test_mode. This is the stateless function used for both NLSQ optimization and Bayesian inference.

Parameters:

  • X (array-like) – Independent variable

  • params (array-like) – Parameter values: [G, tau_b, tau_a, nu, f_B_eq, eta_s]

  • test_mode (str, optional) – Override stored test mode

Returns:

Predicted response

Return type:

jnp.ndarray

predict_flow_curve(gamma_dot, return_components=False)[source]

Predict steady shear stress and viscosity.

Parameters:

  • gamma_dot (ndarray) – Shear rate array (1/s)

  • return_components (bool) – If True, return (sigma, eta, N1)

Returns:

Shear stress σ (Pa), or (σ, η, N₁) if return_components=True

Return type:

ndarray | tuple[ndarray, ndarray, ndarray]

predict_saos(omega, return_components=True)[source]

Predict SAOS storage and loss moduli.

In the linear regime, loop-bridge model reduces to effective Maxwell: G_eff = f_B_eq·G, τ_eff = τ_b

Parameters:

  • omega (ndarray) – Angular frequency array (rad/s)

  • return_components (bool) – If True, return (G’, G’’)

Returns:

(G’, G’’) if return_components=True, else |G*|

Return type:

tuple[ndarray, ndarray] | ndarray

predict_normal_stresses(gamma_dot)[source]

Predict first and second normal stress differences.

N₁ = f_B·G·(S_xx - S_yy) N₂ = 0 (upper-convected derivative)

Parameters:

gamma_dot (ndarray) – Shear rate array (1/s)

Returns:

(N₁, N₂) in Pa

Return type:

tuple[ndarray, ndarray]

simulate_startup(t, gamma_dot, return_bridge_fraction=False)[source]

Simulate startup flow at constant shear rate.

Parameters:

  • t (ndarray) – Time array (s)

  • gamma_dot (float) – Applied shear rate (1/s)

  • return_bridge_fraction (bool) – If True, also return f_B(t)

Returns:

Shear stress σ(t), or (σ, f_B) if return_bridge_fraction=True

Return type:

ndarray | tuple[ndarray, ndarray]

simulate_relaxation(t, gamma_dot_preshear, return_bridge_fraction=False)[source]

Simulate stress relaxation after cessation of steady shear.

Parameters:

  • t (ndarray) – Time array (s), starting from t=0 (cessation)

  • gamma_dot_preshear (float) – Shear rate before cessation (1/s)

  • return_bridge_fraction (bool) – If True, also return f_B(t)

Returns:

Relaxing stress σ(t), or (σ, f_B) if return_bridge_fraction=True

Return type:

ndarray | tuple[ndarray, ndarray]

simulate_creep(t, sigma_applied, return_rate=False)[source]

Simulate creep deformation under constant stress.

Parameters:

  • t (ndarray) – Time array (s)

  • sigma_applied (float) – Applied constant stress (Pa)

  • return_rate (bool) – If True, also return shear rate γ̇(t)

Returns:

Strain γ(t), or (γ, γ̇) if return_rate=True

Return type:

ndarray | tuple[ndarray, ndarray]

simulate_laos(t, gamma_0, omega, n_cycles=None)[source]

Simulate Large-Amplitude Oscillatory Shear (LAOS).

Parameters:

  • t (ndarray) – Time array (s), or None to auto-generate

  • gamma_0 (float) – Strain amplitude (dimensionless)

  • omega (float) – Angular frequency (rad/s)

  • n_cycles (int | None) – Number of oscillation cycles (overrides t)

Returns:

Dictionary with keys: ‘t’, ‘strain’, ‘stress’, ‘strain_rate’, ‘f_B’

Return type:

dict[str, ndarray]

get_bridge_fraction_vs_rate(gamma_dot)[source]

Compute steady-state bridge fraction vs shear rate.

Parameters:

gamma_dot (ndarray) – Shear rate array (1/s)

Returns:

(gamma_dot, f_B_steady)

Return type:

tuple[ndarray, ndarray]

extract_laos_harmonics(laos_result, n_harmonics=5)[source]

Extract Fourier harmonics from LAOS stress response.

Parameters:

  • laos_result (dict[str, ndarray]) – Result from simulate_laos()

  • n_harmonics (int) – Number of harmonics to extract

Returns:

Dictionary with ‘n’, ‘sigma_prime’, ‘sigma_double_prime’, ‘intensity’, ‘I3_I1’

Return type:

dict[str, ndarray]

__repr__()

Return string representation.

Return type:

str

static compute_stretch(S_xx, S_yy, S_zz)

Compute stretch ratio from conformation tensor.

stretch = sqrt(tr(S)/3)

At equilibrium (S=I): stretch = 1. stretch > 1 indicates chains are extended beyond equilibrium.

Parameters:

  • S_xx (float) – Diagonal conformation tensor components

  • S_yy (float) – Diagonal conformation tensor components

  • S_zz (float) – Diagonal conformation tensor components

Returns:

Stretch ratio (dimensionless, ≥ 0)

Return type:

float

deborah_number(omega)

Compute Deborah number De = τ_b·ω.

Parameters:

omega (float) – Angular frequency (rad/s)

Returns:

Deborah number (dimensionless)

Return type:

float

fit(X, y, method='nlsq', check_compatibility=False, use_log_residuals=None, use_multi_start=None, n_starts=5, perturb_factor=0.3, deformation_mode=None, poisson_ratio=0.5, auto_init=False, return_result=False, check_physics=False, uncertainty=None, **kwargs)

Fit the model to data using NLSQ optimization.

This method uses NLSQ (GPU-accelerated nonlinear least squares) by default for fast point estimation. The optimization result is stored for potential warm-starting of Bayesian inference.

For very wide frequency ranges (>10 decades), multi-start optimization is automatically enabled to escape local minima.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – Target values

  • method (str) – Optimization method (‘nlsq’ by default for compatibility)

  • check_compatibility (bool) – Whether to check model-data compatibility before fitting. If True, warns when model may not be appropriate for data. Default is False for backward compatibility.

  • use_log_residuals (bool | None) – Whether to use log-space residuals for fitting. Recommended for wide frequency ranges (>8 decades) to prevent optimizer bias. If None (default), automatically detected based on data range. Explicit True/False overrides auto-detection.

  • use_multi_start (bool | None) – Whether to use multi-start optimization to escape local minima. Recommended for very wide ranges (>10 decades). If None (default), automatically enabled for >10 decades.

  • n_starts (int) – Number of random starts for multi-start optimization (default: 5)

  • perturb_factor (float) – Perturbation magnitude for multi-start random starts (default: 0.3). Parameters are perturbed by ± perturb_factor * (value or range). Larger values (0.7-0.9) explore wider parameter space.

  • auto_init (bool) – If True, calls auto_p0() to estimate initial parameters from data before running the optimizer (default: False).

  • return_result (bool) – If True, returns a FitResult instead of self. This intentionally breaks method chaining for workflows that need structured result objects (default: False).

  • check_physics (bool) – If True, runs post-fit physics validation and emits RheoJaxPhysicsWarning for any violations (default: False).

  • uncertainty (str | None) – Post-fit uncertainty method. "hessian" for fast Cramér-Rao bounds, "bootstrap" for residual bootstrap CIs, or None to skip (default: None).

  • **kwargs – Additional fitting options passed to _fit()

Return type:

BaseModel | Any

Returns:

self for method chaining (default), or FitResult if return_result=True.

Example

>>> model = Maxwell()
>>> model.fit(t, G_data)  # Uses NLSQ by default
>>> model.fit(t, G_data, method='nlsq', max_iter=1000)
>>> model.fit(t, G_data, check_compatibility=True)  # Check compatibility
>>> model.fit(omega, G_star, use_log_residuals=True)  # Force log-residuals
>>> model.fit(mastercurve, None, use_multi_start=True, n_starts=10)  # Multi-start
>>> result = model.fit(t, G_data, return_result=True)  # Structured result
>>> result = model.fit(t, G_data, auto_init=True, check_physics=True,
...                    return_result=True)  # Full pipeline
fit_bayesian(X, y=None, num_warmup=1000, num_samples=2000, num_chains=4, initial_values=None, test_mode=None, seed=None, deformation_mode=None, poisson_ratio=0.5, **nuts_kwargs)

Perform Bayesian inference using NumPyro NUTS sampler.

This method delegates to BayesianMixin.fit_bayesian() to run NUTS sampling for Bayesian parameter estimation. If initial_values is not provided and the model has been previously fitted with fit(), the NLSQ point estimates are automatically used for warm-starting.

Multi-chain sampling is enabled by default (num_chains=4) to provide reliable convergence diagnostics (R-hat, ESS) and parallel execution on multi-GPU systems.

Parameters:

  • X (numpy.typing.ArrayLike) – Independent variable data (input features) or RheoData object

  • y (Optional[numpy.typing.ArrayLike]) – Dependent variable data (observations to fit). If X is RheoData, y is ignored and extracted from X.

  • num_warmup (int) – Number of warmup/burn-in iterations (default: 1000)

  • num_samples (int) – Number of posterior samples per chain (default: 2000)

  • num_chains (int) – Number of MCMC chains (default: 4). Multiple chains enable proper R-hat computation and parallel execution. Chain method is auto-selected: ‘parallel’ on multi-GPU, ‘vectorized’ on single GPU/CPU.

  • initial_values (dict[str, float] | None) – Optional dict of initial parameter values for warm-start. If None and model is fitted, uses NLSQ estimates.

  • test_mode (str | None) – Explicit test mode (e.g., ‘relaxation’, ‘creep’, ‘oscillation’). If None, inferred from RheoData.metadata[‘test_mode’] or defaults to ‘relaxation’. Overrides RheoData metadata if provided.

  • seed (int | None) – Random seed for reproducibility. If None, uses seed=0 for deterministic results. Set to different values for independent runs.

  • **nuts_kwargs – Additional arguments passed to NUTS sampler (e.g., target_accept_prob, chain_method)

Return type:

BayesianResult

Returns:

BayesianResult containing posterior samples, summary statistics, and convergence diagnostics (R-hat, ESS, divergences)

Example

>>> model = Maxwell()
>>> # Warm-start from NLSQ with explicit mode
>>> model.fit(t, G_data, test_mode='relaxation')  # NLSQ optimization
>>> result = model.fit_bayesian(t, G_data, test_mode='relaxation')
>>>
>>> # RheoData with embedded mode (recommended)
>>> rheo_data = RheoData(x=omega, y=G_star, metadata={'test_mode': 'oscillation'})
>>> result = model.fit_bayesian(rheo_data)
>>>
>>> # Or provide explicit initial values
>>> result = model.fit_bayesian(
...     t, G_data,
...     initial_values={'G0': 1e5, 'eta': 1e3},
...     test_mode='creep'
... )
fit_predict(X, y, **kwargs)

Fit model and return predictions.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – Target values

  • **kwargs – Additional fitting options

Return type:

numpy.typing.ArrayLike

Returns:

Model predictions on training data

classmethod from_dict(data)

Create model from dictionary.

Parameters:

data (dict[str, Any]) – Dictionary representation

Return type:

BaseModel

Returns:

Model instance

get_bayesian_result()

Get stored Bayesian inference result.

Return type:

BayesianResult | None

Returns:

BayesianResult from fit_bayesian(), or None if not run

Example

>>> model.fit_bayesian(t, G_data)
>>> result = model.get_bayesian_result()
>>> print(result.diagnostics['r_hat'])
get_credible_intervals(posterior_samples, credibility=0.95)

Compute highest density intervals (HDI) for posterior samples.

Computes the highest posterior density interval for each parameter, which is the shortest interval containing the specified probability mass. This is preferred over equal-tailed intervals for skewed distributions.

Parameters:

  • posterior_samples (dict[str, ndarray]) – Dictionary mapping parameter names to posterior sample arrays (from BayesianResult.posterior_samples)

  • credibility (float) – Probability mass to include in interval (default: 0.95) Common values: 0.68 (1 sigma), 0.95 (2 sigma), 0.997 (3 sigma)

Return type:

dict[str, tuple[float, float]]

Returns:

Dictionary mapping parameter names to (lower, upper) credible interval tuples. All values are float64.

Example

>>> result = model.fit_bayesian(X, y)
>>> intervals_95 = model.get_credible_intervals(
...     result.posterior_samples, credibility=0.95
... )
>>> print(f"95% CI for a: {intervals_95['a']}")
static get_equilibrium_conformation()

Return equilibrium conformation tensor S_eq = I.

In the absence of flow, chains adopt their equilibrium configuration: S = I (isotropic, unstretched).

Returns:

Equilibrium state [S_xx, S_yy, S_zz, S_xy] = [1, 1, 1, 0]

Return type:

Array

get_nlsq_result()

Get stored NLSQ optimization result.

Returns:

OptimizationResult from NLSQ fit, or None if not fitted

Example

>>> model.fit(t, G_data)
>>> result = model.get_nlsq_result()
>>> if result:
...     print(f"Converged: {result.success}")
get_parameter_dict()

Get all parameters as a dictionary.

Returns:

Dictionary of parameter name -> value

Return type:

dict[str, float]

get_parameter_uncertainties()

Get standard errors for fitted parameters from NLSQ covariance.

Returns:

std_error}, or None if covariance unavailable

Return type:

dict of {param_name

get_params(deep=True)

Get model parameters.

Parameters:

deep (bool) – If True, return parameters of sub-objects

Return type:

dict[str, Any]

Returns:

Dictionary of parameter names and values

initialize_from_flow_curve(gamma_dot, sigma)

Initialize parameters from flow curve data.

Parameters:

  • gamma_dot (ndarray) – Shear rate array (1/s)

  • sigma (ndarray) – Stress array (Pa)

Return type:

None

initialize_from_saos(omega, G_prime, G_double_prime)

Initialize parameters from SAOS data.

Uses the crossover frequency and modulus to estimate τ_b and G. In the linear regime, TNT reduces to Maxwell: crossover at ω_c·τ_b = 1.

Parameters:

  • omega (ndarray) – Angular frequency array (rad/s)

  • G_prime (ndarray) – Storage modulus G’ (Pa)

  • G_double_prime (ndarray) – Loss modulus G’’ (Pa)

Return type:

None

property pcov_

Parameter covariance matrix from NLSQ fit.

Returns:

ndarray of shape (n_params, n_params), or None if not fitted

property popt_

Optimal parameter values from NLSQ fit.

Returns:

ndarray of shape (n_params,), or None if not fitted

precompile(test_mode='relaxation', X=None, y=None)

Precompile NLSQ residual functions to eliminate JIT cold-start.

Triggers JIT compilation by running a minimal fit (max_iter=1) with dummy data. The model parameters are reset afterwards so the model is left in its original state.

This is useful for interactive sessions or benchmarks where the ~870ms first-fit JIT overhead should be excluded.

Parameters:

  • test_mode (str) – Test mode to precompile for (default: ‘relaxation’).

  • X (Optional[numpy.typing.ArrayLike]) – Optional input data for shape inference. If None, uses a 10-point logspace array.

  • y (Optional[numpy.typing.ArrayLike]) – Optional output data. If None, generates ones matching X.

Return type:

float

Returns:

Compilation time in seconds.

Example

>>> model = Maxwell()
>>> t = model.precompile(test_mode='relaxation')
>>> print(f"Compiled in {t:.2f}s")
>>> model.fit(X, y)  # No JIT overhead
precompile_bayesian(X=None, y=None, test_mode=None, num_chains=4)

Precompile NUTS kernel to eliminate JIT overhead in subsequent calls.

Triggers JIT compilation of the NumPyro model by running a minimal sampling (1 warmup, 1 sample). This caches the compiled kernel so that subsequent fit_bayesian() calls are 2-5x faster.

Parameters:

  • X (ndarray | RheoData | None) – Sample input data for determining array shapes. If None, uses a default 10-point linspace [0.01, 100].

  • y (ndarray | None) – Sample output data. If None, generates dummy data.

  • test_mode (str | TestModeEnum | None) – Test mode to precompile for. If None, defaults to ‘relaxation’.

Returns:

Compilation time in seconds.

Return type:

float

Example

>>> model = Maxwell()
>>> compile_time = model.precompile_bayesian(test_mode='relaxation')
>>> print(f"Compiled in {compile_time:.1f}s")
>>> # Now fit_bayesian() will be faster
>>> result = model.fit_bayesian(X, y)  # No compilation overhead
predict(X, test_mode=None, deformation_mode=None, poisson_ratio=None, **kwargs)

Make predictions.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • test_mode (str | None) – Optional test mode (‘oscillation’, ‘relaxation’, ‘creep’, ‘flow’). If provided, sets model’s test_mode before prediction. Useful for data generation without fitting.

  • deformation_mode (str | DeformationMode | None) – Optional deformation mode for output conversion. If None, uses the mode stored from fit(). If tensile, converts G* predictions to E* space.

  • poisson_ratio (float | None) – Poisson’s ratio for conversion. If None, uses value stored from fit() (default 0.5).

  • **kwargs – Additional arguments passed to the internal _predict method.

Return type:

numpy.typing.ArrayLike

Returns:

Model predictions (in E* space if deformation_mode is tensile)

sample_prior(num_samples=1000, seed=None)

Sample from prior distributions over parameter bounds.

Samples from uniform prior distributions defined by parameter bounds. This is useful for prior predictive checks and understanding the prior’s influence on the posterior.

Parameters:

  • num_samples (int) – Number of samples to draw from prior (default: 1000)

  • seed (int | None) – Random seed for reproducibility (default: None)

Return type:

dict[str, ndarray]

Returns:

Dictionary mapping parameter names to arrays of prior samples. Each array has shape (num_samples,) and dtype float64.

Raises:

Example

>>> model = MyModel()
>>> prior_samples = model.sample_prior(num_samples=500, seed=42)
>>> print(prior_samples["a"].shape)  # (500,)
score(X, y)

Compute model score (R² by default).

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – True target values

Return type:

float

Returns:

Model score (R² coefficient)

set_parameter_dict(params)

Set parameters from a dictionary.

Parameters:

params (dict[str, float]) – Dictionary of parameter name -> value

Return type:

None

set_params(**params)

Set model parameters.

Parameters:

**params – Parameter names and values

Return type:

BaseModel

Returns:

self for method chaining

to_dict()

Serialize model to dictionary.

Return type:

dict[str, Any]

Returns:

Dictionary representation of model

weissenberg_number(gamma_dot)

Compute Weissenberg number Wi = τ_b·γ̇.

Parameters:

gamma_dot (float) – Shear rate (1/s)

Returns:

Weissenberg number (dimensionless)

Return type:

float

Key Methods:

  • fit(rheo_data): Fit model to SAOS, flow curve, or transient data using NLSQ

  • predict(x, test_mode, **kwargs): Predict stress, viscosity, or modulus

  • fit_bayesian(rheo_data, num_warmup, num_samples): Bayesian inference with NumPyro NUTS

  • get_parameter_values(): Return dict of fitted parameters

  • set_parameter(name, value): Update parameter value

Test Mode Support:

Test Mode

Description

FLOW_CURVE

Steady shear: \(\eta(\dot{\gamma})\) with non-monotonic behavior

STARTUP

Transient shear: \(\sigma(t)\) with overshoot and possible thickening

CREEP

Constant stress: \(\gamma(t)\) with retardation and flow

RELAXATION

Stress relaxation: \(\sigma(t)\) after step strain (bridge-to-loop conversion)

OSCILLATION

SAOS: \(G'(\omega)\), \(G''(\omega)\) with two relaxations

LAOS

Large amplitude: \(\sigma(t)\) with higher harmonics

18. References

Foundational Papers

  1. Annable, T., Buscall, R., Ettelaie, R., & Whittlestone, D. (1993) “The rheology of solutions of associating polymers: Comparison of experimental behavior with transient network theory” Journal of Rheology, 37(4), 695-726. DOI: 10.1122/1.550391

    Original TNT loop-bridge model for HEUR solutions.

  2. Tanaka, F., & Edwards, S. F. (1992) “Viscoelastic properties of physically cross-linked networks: Transient network theory” Macromolecules, 25(5), 1516-1523. DOI: 10.1021/ma00031a024

    Foundation of transient network theory.

  3. Marrucci, G., Bhargava, S., & Cooper, S. L. (1993) “Models of shear-thickening behavior in physically crosslinked networks” Macromolecules, 26(24), 6483-6488. DOI: 10.1021/ma00076a027

    Theory of flow-induced association in telechelic polymers.

Force-Dependent Kinetics

  1. Bell, G. I. (1978) “Models for the specific adhesion of cells to cells” Science, 200(4342), 618-627. DOI: 10.1126/science.347575

    Bell’s theory of force-dependent bond dissociation.

  2. Vaccaro, A., & Marrucci, G. (2000) “A model for the nonlinear rheology of associating polymers” Journal of Non-Newtonian Fluid Mechanics, 92(2-3), 261-273. DOI: 10.1016/S0377-0257(00)00095-1

    Stress-population coupling in telechelic networks.

Experimental Validation

  1. Tripathi, A., Tam, K. C., & McKinley, G. H. (2006) “Rheology and dynamics of associative polymers in shear and extension: Theory and experiments” Macromolecules, 39(5), 1981-1999. DOI: 10.1021/ma051614x

    Comprehensive rheological study of HEUR solutions, validation of loop-bridge kinetics.

  2. Pellens, L., Corrales, R. G., & Mewis, J. (2004) “General nonlinear rheological behavior of associative polymers” Journal of Rheology, 48(2), 379-393. DOI: 10.1122/1.1645516

    Shear thickening and thinning in associating polymers.

Reviews and Extensions

  1. Rubinstein, M., & Semenov, A. N. (2001) “Dynamics of entangled solutions of associating polymers” Macromolecules, 34(4), 1058-1068. DOI: 10.1021/ma0013049

    Theory of dynamics in telechelic networks with entanglements.

  2. Berret, J. F. (2006) “Rheology of wormlike micelles: Equilibrium properties and shear banding transitions” Molecular Gels, 667-720. DOI: 10.1007/1-4020-3689-2_20

    Related systems: wormlike micelles with loop-bridge-like kinetics.

  3. Green, M. S., & Tobolsky, A. V. (1946) “A new approach to the theory of relaxing polymeric media” Journal of Chemical Physics, 14(2), 80-92. DOI: 10.1063/1.1724109

    Classic transient network theory (Green-Tobolsky model).

Numerical Methods

  1. Fang, J., Kröger, M., & Öttinger, H. C. (2000) “A thermodynamically admissible reptation model for fast flows of entangled polymers. II. Model predictions for shear and extensional flows” Journal of Rheology, 44(6), 1293-1317. DOI: 10.1122/1.1308522

    Numerical methods for stiff differential equations in rheology.