Closed-Form Critical-State Caputo Flow for a Teardrop Bounding Surface: Operator Semantics and Factorial Consequences

1. Introduction

Fractional-order plasticity entered soil mechanics along two connected lines: the fractional viscoplasticity of Sumelka [1,2,3], which generalized the direction of plastic flow to a non-local fractional gradient and showed that non-normality and induced plastic anisotropy follow from the operator itself [4], and the critical-state fractional flow rules of Sun, Shen and co-workers [5,6,7], developed across applications from ballast to bounding-surface, numerical-scheme and true-triaxial formulations [8,9,10,11,12,13] and extended to viscoplastic soils and anisotropically overconsolidated clay [14,15], with the dynamics of the fractional viscoplastic operator examined numerically in [16]; a related non-orthogonal flow-rule line was developed by Lu and co-workers [17]. In the classification of the recent review [18], these are the past-stress and future-critical-state branches of stress-fractional plasticity. The critical-state branch anchors the operator so that the integration limits are the current and the critical stress states, which merge as the critical state is approached [9], and demonstrated — for the Modified Cam-Clay (MCC) surface [19] within the critical-state framework [20] — that state-dependent non-associativity then emerges without any plastic potential [9,11]. The analytical appeal of that line rests on a specific coincidence: MCC is polynomial in linear stress, so its Caputo fractional gradients reduce to incomplete Beta and Gamma functions in closed form [9,18].

The companion study [21] embedded a fractional flow operator in a 2×2×2 factorial with the teardrop bounding surface of Chatwong et al. [22] and a nonlinear AJOP compression law [23], and found the flow factor to be the single largest main effect (≈ −60% of the baseline plastic volumetric response at OCR = 5 for two independently calibrated clays). That study was, however, deliberately explicit that its operator was a fixed-window, past-stress Grünwald–Letnikov (GL) scheme — a loading-memory semantics in the past-stress branch of [18] — and not the critical-state Caputo definition [9], which it retained as motivation only. Two questions follow immediately. First, can the critical-state Caputo operator be constructed exactly for the teardrop surface at all, given that the surface is not polynomial in linear stress? Second, do the factorial conclusions of the companion study survive the change of window semantics?

This paper answers both. The first answer is a closed form: on the log-stress axis — formally, a Caputo derivative taken with respect to another function [24], here the logarithm, in the weighted/conjugated sense of [25] — the teardrop takes a power–exponential form, its critical-state terminal is exactly the reciprocal shape parameter, and its Caputo gradient is an incomplete-Beta–Kummer or Humbert-Φ₁ expression — the same class of special functions that made the MCC results attractive, now on a surface those results cannot reach. The second answer is a collapse: with the critical-state operator in place, the dominant flow main effect of the companion study falls below 1% at every tested state, for both soils and both loading paths, and is shown to be a transient that vanishes as the state reaches and dwells at the critical state. The contributions are: (i) the exact critical-state terminal t* = 1/Ψ and the closed-form Caputo gradients of the teardrop (Section 2), verified to near machine precision; (ii) a flow rule built on them that is exactly associated at the critical state (Section 2.3); (iii) a paired factorial re-evaluation showing the collapse of the flow main effect and its migration into the geometry×flow interaction (Section 4); and (iv) a clean discriminant between the two available routes to state-dependent non-associativity (Section 5).

Relation to the companion study, stated plainly: the calibrated parameters, path definitions, hardening semantics, factorial index definitions and the self-testing reference implementation are inherited from [21] to guarantee one-to-one comparability. Everything else is new to this paper: the operator belongs to the other branch of stress-fractional plasticity (critical-state window semantics rather than past-stress memory [18]); the closed-form ψ-Caputo gradients, the exact terminal t* = 1/Ψ and their 1,240-case verification are new analytical results; the principal numerical finding reverses the headline of [21]; and the dwell mechanism, the endpoint identity and the phase-transformation discriminant are new. The companion is the control experiment for this paper, not its template.

2. Closed-Form Critical-State Caputo Flow for the Teardrop

2.1. The Teardrop in Log-Stress Space and Its Exact Critical-State Terminal

In normalized triaxial space the teardrop yield curve of Chatwong et al. [22] is $\overline{q}$ = g($\overline{p}$) = $\overline{p}$·(Ω ln(1/$\overline{p}$))^1/Ψ, with $\overline{p}$ = p′/p′c, $\overline{q}$ = q/(Mp′c) and shape parameters Ψ, Ω, where M is the critical-state stress ratio of Table 2. All stresses and dilatancies are hereafter M-normalized in this sense — η̅ = η/M and D̅ = D/M, equivalently M ≡ 1 in the canonical implementation (ff_canonical.py in Supplementary S1) — with physical values recovered as η = Mη̅ and D = MD̅; bars mark the normalized stress variables, while the dilatancy symbols D, D_int, D_cap and D_GL denote M-normalized values throughout. Substituting the log-stress variable t = ln(1/$\overline{p}$) gives

g̃(t) = e^−t · (Ω t)^1/Ψ, t ≥ 0,

(1)

a power–exponential form: the algebraic component (Ωt)^1/Ψ is a pure power of t — the substitution $\overline{p}$ = e^−t converts the logarithmic composition into this power — while the factor e^−t is the exponential image of $\overline{p}$ itself; it is precisely this structure (power × exponential against the Caputo kernel) that renders every required integral closed-form — the log-stress axis being the natural axis of soil compression, on which the classical and AJOP compression laws are themselves defined [19,20,23]. The critical state of the surface — the apex, at which the associated flow direction is purely deviatoric — satisfies dg̃/dt = 0, which reduces to Ω t = Ω/Ψ, i.e.,

$$t^{*}=1/\Psi \mathit{exactly}, {\overline{p}}^{*}=e^{-1/\Psi }, \mathit{independent\; of} \Omega$$

(2)

For the four soils of the companion calibration this gives $\overline{p}$* = 0.490 (Ψ=1.4, Boston Blue and London Clay), 0.513 (Ψ=1.5, kaolin) and 0.386 (Ψ=1.05, Lower Cromer Till). We verified the closed form against direct numerical maximization to 10⁻⁸ for Ψ ∈ [1.05, 2.0]. This point is a property of the surface (the associated-flow critical state) and must not be confused with the zero of the embedded non-associated rule of [22], which lies at $\overline{p}$ ≈ 0.341 for the Boston Blue Clay shape (Section 5).

2.2. Closed-Form Caputo Integrals

With f(t) = $\overline{q}$₀ − g̃(t) at fixed $\overline{q}$₀ (the Caputo construction annihilates the constant), the fractional gradient — in the classical Caputo sense [26,27,28] — requires integrals of g̃′(s) = e^−s·Ω^1/Ψ·[(1/Ψ)s^1/Ψ−1 − s^1/Ψ] against the kernel |x−s|^−α, i.e., two instances of the master integral with β = 1/Ψ and β = 1/Ψ − 1. For the tip terminal (a = 0) the master integral is classical (entry 3.383.1 of [29]):

∫₀ᵗ s^β e^−s (t−s)^−α ds = B(β+1, 1−α) · t^β+1−α · ₁F₁(β+1; β+2−α; −t).

(3)

For the critical-state terminal a = t* > 0 we derived the shifted forms by the substitution s = a + (x−a)u and identification with the Euler-type integral representation of the Humbert confluent hypergeometric function Φ₁ [30,31]:

dry side (x > a): ∫ₐˣ s^β e^−s(x−s)^−α ds = a^β e^−a (x−a)^1−α/(1−α) · Φ₁(1, −β; 2−α; −(x−a)/a, −(x−a)),

(4)

wet side (x < a): ∫ₓᵃ s^β e^−s(s−x)^−α ds = x^β e^−x (a−x)^1−α/(1−α) · Φ₁(1−α, −β; 2−α; −(a−x)/x, −(a−x)).

(5)

Φ₁ was evaluated by its single-series expansion in the exponential variable with Gauss ₂F₁ factors, which remains valid under the analytic continuation of ₂F₁ for the parameter regimes encountered here (first argument down to −2.3). Verification against singularity-aware weighted quadrature [32] covered 400 cases for the tip-terminal form (3) over the teardrop exponents γ ∈ {1/Ψ − 1, 1/Ψ} (worst relative error 1.2×10⁻¹⁵), 200 cases at a second shifted terminal a = ln 2 in the moderate-|Z| regime with representative exponents γ ∈ {0.4, 0.7, 1.0, 1.5} (worst 7.4×10⁻¹⁶), and 640 cases for the critical-state terminal spanning both sides, both companion critical states, Ψ ∈ [1.05, 2.0] and α ∈ [0.3, 0.95] (worst 5.1×10⁻¹²); the a = ln 2 block stress-tests the master integrals at an additional shift value and is distinct from the dilatancy of the classical geometry, which is obtained from quadrature tables (Section 3.1). Figure 2 shows the full error distribution of a 320-case subset regenerated for this paper.

2.3. Flow Rule and Its Limits

The operator is defined as the terminal-directed, orientation-corrected Caputo operator in the log-stress coordinate t = −ln $\overline{p}$, taken toward the critical-state terminal t* = 1/Ψ (Section 2.2). It may equally be written as an orientation-corrected extension of the ψ-Caputo construction of Almeida [24] for the decreasing map ψ($\overline{p}$) = ln(1/$\overline{p}$): the kernel that keeps the integrand real and positive is the orientation-corrected form (ψ($\overline{p}$)−ψ(u))^−α on the dry side (t > t*) and (ψ(u)−ψ($\overline{p}$))^−α on the wet side (t < t*), rather than the increasing-ψ kernel of [24]. The two forms carry identical information: under u = e^−s the change of variable maps the Caputo integral in t onto this ψ-Caputo integral in $\overline{p}$, and as α → 1 the operator returns dF/dt = −$\overline{p}$ f′($\overline{p}$), i.e., associated flow in the log-stress (t–q) conjugate representation. The choice ψ = ln(1/$\overline{p}$) is the modeling postulate, argued from the fact that soil compression laws — classical and AJOP alike — are native to the logarithmic stress axis [19,20,23]. In the variable t = ψ($\overline{p}$) these are exactly the terminal-directed, orientation-corrected Caputo operators of Section 2.2, which is how the closed forms arise. The flow vector is assembled as follows. The yield function f($\overline{p}$,$\overline{q}$) = $\overline{q}$ − g($\overline{p}$) is linear in $\overline{q}$, so its deviatoric flow component is the integer derivative ∂f/∂$\overline{q}$ = 1, left untouched; the volumetric component replaces ∂f/∂ψ (the log-conjugate volumetric direction, work-paired with the increment d ln $\overline{p}$) by the terminal-directed, orientation-corrected Caputo operator of Proposition 1, left-sided on the dry side (t > t*) and right-sided on the wet side (t < t*). The Caputo construction annihilates the constant $\overline{q}$₀, so only integrals of g̃′ appear. The reported dilatancy is the component ratio D_cap = (volumetric)/(deviatoric) in this convention. The transformation to the $\overline{p}$-conjugate pairing is stated with its sign: dψ/d$\overline{p}$ = −1/$\overline{p}$ = −e^t, so the covector components obey ∂f/∂$\overline{p}$ = −e^t·∂f/∂ψ; because D_cap is built from g̃′ = −f′, the $\overline{p}$-conjugate dilatancy is D̅_p̅ = −e^t·D_cap: the factor −e^t < 0 is an orientation-reversing rescaling, converting the dilative-positive t-conjugate convention of Proposition 1 to the contractive-positive convention of the evolution equations while moving no zero; the physical dilatancy then follows as D = M·D̅_p̅, where M ≡ 1 in the normalized coordinates used throughout. Rescaling the flow magnitude does, however, alter the stress path and the hardening evolution en route, so convention-independence of the accumulated response is not a corollary of the zero analysis alone; it follows instead from an endpoint identity: the hardening law ties the increments exactly, dε_v^p = λ_eff(p′c)·d ln p′c, so the accumulated plastic volumetric strain of any converged configuration equals ∫ λ_eff d ln p′c between the initial and terminal states — a function of the endpoints only, independent of the dilatancy magnitude en route. Stated with its conditions: under the present hardening law and loading/convergence conditions, two configurations that reach the same absolute terminal state accumulate the same plastic volumetric response, irrespective of the path-dependent dilatancy magnitude and of the conjugacy convention. On the constant-p′ path used here, sharing the normalized terminal $\overline{p}$* does fix the absolute terminal, p′c* = p′/$\overline{p}$*, so the condition is met by construction; for other loading controls the coincidence of absolute terminals must be established rather than assumed. The GL window differs precisely because its effective terminal is displaced. The rerun confirms the identity in the cases checked: re-running the full Caputo factorial in the $\overline{p}$-conjugate convention leaves every flow-independent index unchanged and reduces the flow main effect to a numerically-zero contrast at OCR = 5 (|ΔF| < 0.1% of baseline for both soils; the exact last digit of such a near-zero contrast is floating-point-ordering dependent across platforms and is not a determinism target), while the terminal states of the integer and Caputo classical configurations agree to three decimals ($\overline{p}$ = 0.4996 versus 0.5000, accumulated strain −0.1375 versus −0.1377). Reported magnitudes of D_cap remain convention-specific and are labelled as such.

Proposition 1 (Closed-form critical-state Caputo flow). Let $t=-\mathrm{l}\mathrm{n}\overline{p}$, with critical-state terminal $t^{*}=1/\Psi$, exponent $\beta =1/\Psi -1$ and constant $c={\Omega }^{1/\Psi }$, so that $\tilde{g}\left(t\right)=c\hspace{0.17em}{t}^{\beta +1}{e}^{-t}$ and ${\tilde{g}}^{\prime }\left(t\right)=c\hspace{0.17em}{e}^{-t}\left[\left(\beta +1\right)t^{\beta }-t^{\beta +1}\right]$ with ${\tilde{g}}^{\prime }\left(t^{*}\right)=0$. The volumetric flow component is $D_{\mathrm{c}\mathrm{a}\mathrm{p}}\left(t\right)=-D_{t^{*}}^{\alpha }\tilde{g}\left(t\right)$, where $D_{t^{*}}^{\alpha }$ is the terminal-directed, orientation-corrected Caputo operator toward $t^{*}$ — left-sided for $t>t^{*}$, right-sided for $t<t^{*}$, with the orientation sign fixed so that its integer limit is the associated dilatancy on both sides (deviatoric component $\partial f/\partial q=1$; the constant $q$ is annihilated). Then

$$D_{\mathrm{c}\mathrm{a}\mathrm{p}}\left(t\right)=-{\displaystyle \frac{c}{\Gamma \left(1-\alpha \right)}}\left[\left(\beta +1\right)\hspace{0.17em}{I}_{\beta }\left(t\right)-I_{\beta +1}\left(t\right)\right], t>t^{*},$$

(6)

$$I_{\gamma }\left(t\right)={\left(t^{*}\right)}^{\gamma }{e}^{-t^{*}}\hspace{0.17em}{\displaystyle \frac{{\left(t-t^{*}\right)}^{1-\alpha }}{1-\alpha }}\hspace{0.17em}{\Phi }_1\left(1,\hspace{0.17em}-\gamma ;\hspace{0.17em}2-\alpha ;\hspace{0.17em}-{\displaystyle \frac{t-t^{*}}{t^{*}}},\hspace{0.17em}-\left(t-t^{*}\right)\right),$$

$$D_{\mathrm{c}\mathrm{a}\mathrm{p}}\left(t\right)=-{\displaystyle \frac{c}{\Gamma \left(1-\alpha \right)}}\left[\left(\beta +1\right)\hspace{0.17em}{J}_{\beta }\left(t\right)-J_{\beta +1}\left(t\right)\right], t<t^{*},$$

(7)

$$J_{\gamma }\left(t\right)=t^{\gamma }{e}^{-t}\hspace{0.17em}{\displaystyle \frac{{\left(t^{*}-t\right)}^{1-\alpha }}{1-\alpha }}\hspace{0.17em}{\Phi }_1\left(1-\alpha ,\hspace{0.17em}-\gamma ;\hspace{0.17em}2-\alpha ;\hspace{0.17em}-{\displaystyle \frac{t^{*}-t}{t}},\hspace{0.17em}-\left(t^{*}-t\right)\right),$$

with $\gamma \in \{\beta ,\beta +1\}$ and ${\Phi }_1$ the Humbert confluent function. For every fixed t > 0, the expressions are real-valued and finite in closed form for all $0<\alpha <1$, $\Psi >0$ (including $\beta =1/\Psi -1<0$); they diverge only in the endpoint limit t → 0⁺. Moreover $D_{\mathrm{c}\mathrm{a}\mathrm{p}}\left(t^{*}\right)=0$, so the flow is purely deviatoric at the critical state for every $\alpha$; and as $\alpha \to 1^{-}$, on both sides $D_{\mathrm{c}\mathrm{a}\mathrm{p}}\left(t\right)\to -{\tilde{g}}^{\prime }\left(t\right)$, the integer associated dilatancy in the $t$-conjugate representation — dilative on the dry side, contractive on the wet side. As $t\to 0^{+}$ ($\overline{p}\to 1$) both $D_{\mathrm{c}\mathrm{a}\mathrm{p}}$ and $-{\tilde{g}}_{\prime }$ diverge in step, reflecting the vertical surface tangent at the isotropic axis; this limit lies outside the loading paths, which remain in $\overline{p}\in \left[0.03,0.985\right]$. The $\overline{p}$-conjugate dilatancy is ${\overline{D}}_{\overline{p}}\left(t\right)=-e^t{D}_{\mathrm{c}\mathrm{a}\mathrm{p}}\left(t\right)$, an orientation-reversing rescaling that preserves the zero at $t^{*}$ and returns the contractive-positive convention; the physical dilatancy follows as D = M·D̅_p̅(t), with M ≡ 1 in the present normalized coordinates.

Proof. Since $f=q-\tilde{g}$ is linear in $q$, $D_{\psi }^{\alpha }f=-D_t^{\alpha }\tilde{g}$. Substituting ${\tilde{g}}^{\prime }\left(s\right)=c\hspace{0.17em}{e}^{-s}\left[\left(\beta +1\right)s^{\beta }-s^{\beta +1}\right]$ into the terminal-directed integrals and evaluating term by term with the master integrals (4)–(5) at $a=t^{*}$, $x=t$ and $\gamma \in \{\beta ,\beta +1\}$ gives the stated forms (6) and (7); the common leading sign on both sides is fixed by the orientation correction, so that ${\mathrm{l}\mathrm{i}\mathrm{m}}_{\alpha \to 1^{-}}{D}_{\mathrm{c}\mathrm{a}\mathrm{p}}=-{\tilde{g}}^{\prime }$ throughout.

Three structural properties were verified numerically (Figure 1). (i) As α → 1⁻, D_cap converges to the integer associated dilatancy at every tested state (seven states, agreement to the third decimal at α = 0.999). (ii) Associativity at the critical state, stated precisely: at t = t* both the ψ-Caputo volumetric component (its integration interval having shrunk to zero) and the integer volumetric gradient vanish, while the deviatoric components coincide identically; the full flow vector of the fractional rule therefore equals the associated flow vector — purely deviatoric — exactly at the critical state, for every α, and D_cap changes sign across it (contractive on the wet side, dilative on the dry side). In any neighbourhood of the critical state the two directions differ (the claim is exact coincidence at the terminal, not in its vicinity) — the defining property the fixed-window GL operator lacks, since its window remains finite at the critical state. (iii) The induced non-associativity is strongly non-uniform: at α = 0.641 the ratio D_cap/D_assoc ranges from ≈ 0.37 near the critical state to ≈ 1.13 far on the dry side, a richer structure than any uniform scaling.

Figure 2. Verification of the closed forms: relative error between the Humbert-Φ₁ expressions and singularity-aware weighted quadrature over 320 parameter cases (both sides, Ψ ∈ [1.05,2.0], α ∈ [0.3,0.95]). All cases fall below 10⁻¹¹; the full 1,240-case grid of Section 2.2 attains a worst relative error of 5.1×10⁻¹².

Figure 2. Verification of the closed forms: relative error between the Humbert-Φ₁ expressions and singularity-aware weighted quadrature over 320 parameter cases (both sides, Ψ ∈ [1.05,2.0], α ∈ [0.3,0.95]). All cases fall below 10⁻¹¹; the full 1,240-case grid of Section 2.2 attains a worst relative error of 5.1×10⁻¹².

3. Numerical Setup

3.1. Constitutive Equations and State Variables

For self-containment the governing relations of the companion study [21,22] are summarised here in the normalized coordinates used throughout, identical to those of Section 2: $\overline{p}=p^{\prime }/{p^{\prime }}_c$ and $\overline{q}=q/\left(M{p^{\prime }}_c\right)$; the normalized stress ratio and dilatancy are η̅ = η/M and D̅ = D/M (physical values η = Mη̅, D = MD̅; the canonical implementation sets M ≡ 1 in these coordinates). $\overline{q}=\overline{p}\hspace{0.17em}{\left[\Omega \mathrm{l}\mathrm{n}\left(1/\overline{p}\right)\right]}^{1/\Psi }$. Its integer flow is the integer-order embedded non-associated rule of [22] — not the associated gradient of the locus —

$$D_{\mathrm{i}\mathrm{n}\mathrm{t}}\left(\overline{p}\right)=1-(\Omega \mathrm{l}\mathrm{n}\left(1/\overline{p}\right){)}^{1/\Psi }.$$

The fixed-window Grünwald–Letnikov operator of [21] applies the same window to the locus function $\overline{g}\left(x\right)=x\hspace{0.17em}{\left[\Omega \mathrm{l}\mathrm{n}\left(1/x\right)\right]}^{1/\Psi }$,

$$D_{\mathrm{G}\mathrm{L}}\left(\overline{p}\right)=-h^{-\alpha }{\displaystyle \sum _{k=0}^N}{w}_k^{\left(\alpha \right)}\hspace{0.17em}\overline{g}\left(\overline{p}-kh\right), h=0.005,\hspace{0.25em}N=40,$$

an additive window in $\overline{p}$ with binomial weights $w_k^{\left(\alpha \right)}$; the critical-state Caputo dilatancy is that of Proposition 1. The AJOP compression law [23] sets the compression slope from the current preconsolidation pressure,

$${\lambda }_{\mathrm{A}\mathrm{J}\mathrm{O}\mathrm{P}}\left({p^{\prime }}_c\right)={\displaystyle \frac{a_c}{\mathrm{l}\mathrm{n}10}}\left[1+{\displaystyle \frac{a_c\hspace{0.17em}x}{\sqrt{\theta +{\left(a_{c}x\right)}^2}}}\right], x={\mathrm{l}\mathrm{o}\mathrm{g}}_{10}\left({p^{\prime }}_c/p_r\right),$$

with per-soil $\left(a_c,p_r,\theta \right)$ from [23] (p_r denotes the AJOP reference pressure, relabelled from the source symbol to avoid clashing with the response R). The effective hardening modulus is ${\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}=\lambda -\kappa$, with $\lambda ={\lambda }_{\mathrm{A}\mathrm{J}\mathrm{O}\mathrm{P}}\left({p^{\prime }}_c\right)$ for the AJOP-compression configurations and the constant $\lambda ={\lambda }_{\mathrm{L}\mathrm{F}}$ for the conventional ones; a configuration is hardening-singular — its response undefined and reported as such rather than clamped — when ${\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}\le 0$. The fractional order is state-dependent and frozen at the initial state,

$$\alpha \left(\mathrm{O}\mathrm{C}\mathrm{R}\right)={\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}}+\left(1-{\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\hspace{0.17em}{e}^{-k\left(\mathrm{O}\mathrm{C}\mathrm{R}-1\right)}, {\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}}=0.55,\hspace{0.25em}k=0.40,$$

a modeling assumption retained for comparability with [21]. Along a drained, constant-$p^{\prime }$, shear-strain-controlled path the state and the response $R$ (the M-normalized accumulated plastic volumetric response; the corresponding physical accumulated plastic volumetric strain is R_phys = MR) evolve as

$$\mathrm{d}\overline{p}=-{\displaystyle \frac{\overline{p}\hspace{0.17em}D}{{\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}}}\hspace{0.17em}\mathrm{d}{\epsilon }_s^p, R={\int }_0^{{\epsilon }_s^p}D\hspace{0.17em}\mathrm{d}\epsilon ,$$

while the approximately-undrained path integrates $\mathrm{d}\mathrm{l}\mathrm{n}\overline{p}=-D\hspace{0.17em}B\hspace{0.17em}\mathrm{d}{\epsilon }_s^p$ with $B=\left(1+e_0\right)/\kappa +1/{\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}$. Here ${\epsilon }_s^p$ is a dimensionless numerical loading budget, not a physical strain target: it is set to ${\epsilon }_s^p=6$ — far beyond any laboratory strain — solely to drive every configuration through the critical-state dwell, so that the accumulated-response comparison is taken at a common, fully-developed terminal. The dilatancy D entering these evolution equations is the M-normalized quantity of Section 2.1 (M ≡ 1), so the computed response is likewise M-normalized, R_phys = MR; because all factorial contrasts are normalized by the within-soil baseline, the common factor M cancels identically from every within-soil percentage index; it also leaves equilibrium locations unchanged, because a positive rescaling does not move any dilatancy zero. The drained factorial tables use 30,000 increments (the sweep of record); a 60,000-increment rerun verifies refinement, and the approximately undrained runs use 200,000 increments.

Beyond the constitutive relations of Section 3.1, the numerical machinery follows the companion study [21]: the locked Level-I parameter set (Table 2) [21,22], the incremental drained (constant-p′, shear-strain-controlled, budget εₛᵖ = 6) and approximately undrained paths, the state-dependent fractional order α(OCR) with α_min = 0.55 and k = 0.40 (a modeling assumption, retained for comparability with [21]), the hardening-singularity semantics in which singular configurations raise and are reported as undefined rather than clamped, and the 42-check self-testing reference implementation. The critical-state Caputo dilatancy enters as a third flow level alongside the integer and GL-window operators. Because Φ₁ evaluation per step is costly, D_cap($\overline{p}$) is precomputed on Hölder-aware per-side splines in the variable w = |t − t*|^1−α, which removes the (t−t*)^1−α kink at the critical state; the interpolation error is ≤ 2.4×10⁻⁸ against exact evaluation, and a full path integrated on the table agrees with a path integrated on exact per-step evaluation to 1.2×10⁻¹³ in accumulated strain. The classical reference geometry (the normalized MCC-type surface [19]), whose log-stress form is not a power law, receives its critical-state Caputo dilatancy (terminal t* = ln 2) by the same table machinery built on singularity-aware quadrature; the closed form is a property of the teardrop and is presented as such.

The design is presented as paired 2×2×2 factorials sharing the integer-flow baseline: the pairing {integer ↔ GL-window} reproduces the companion study [21] exactly, and the pairing {integer ↔ Caputo-CS} is new, with configurations C0–C7 mirroring F0–F7. This preserves one-to-one comparability of every index. For self-containment: with R denoting the M-normalized accumulated plastic volumetric response of a configuration and b = |R(C0)|, the indices are ΔG = (R(C1)−R(C0))/b, ΔF = (R(C2)−R(C0))/b, ΔC = (R(C3)−R(C0))/b, the pairwise interactions the standard 2×2 contrasts (e.g., I_GF = (R(C4)−R(C2)−R(C1)+R(C0))/b), and I_GFC the triple contrast; simple effects quoted below are geometry-conditional flow contrasts such as ΔF|teardrop = (R(C4)−R(C1))/b = ΔF + I_GF.

Table 1. Configurations of the Caputo pairing (mirroring F0–F7 of the companion).

Config	Yield geometry	Flow operator	Hardening/compression
C0	Classical	Integer-order	Conventional (λ)
C1	Teardrop	Integer-order (embedded rule of [22])	Conventional (λ)
C2	Classical	Caputo–CS	Conventional (λ)
C3	Classical	Integer-order	AJOP [23]
C4	Teardrop	Caputo–CS	Conventional (λ)
C5	Teardrop	Integer-order	AJOP [23]
C6	Classical	Caputo–CS	AJOP [23]
C7	Teardrop	Caputo–CS	AJOP [23]

Table 1. Configurations of the Caputo pairing (mirroring F0–F7 of the companion).

Config	Yield geometry	Flow operator	Hardening/compression
C0	Classical	Integer-order	Conventional (λ)
C1	Teardrop	Integer-order (embedded rule of [22])	Conventional (λ)
C2	Classical	Caputo–CS	Conventional (λ)
C3	Classical	Integer-order	AJOP [23]
C4	Teardrop	Caputo–CS	Conventional (λ)
C5	Teardrop	Integer-order	AJOP [23]
C6	Classical	Caputo–CS	AJOP [23]
C7	Teardrop	Caputo–CS	AJOP [23]

Table 2. Locked Level-I parameter set of the companion study [21] (single set, whole program; identical to its published Supplementary S2; soil sources as compiled in [21,22]). The λ and κ listed here are the source-calibration constants of [22]; the integrator takes their ratio κ/λ and applies it to the AJOP-consistent compression index of [21], giving λ = 0.187, κ = 0.0365 for Boston Blue Clay and λ = 0.0586, κ = 0.0223 for London Clay. The normalized constant-λ indices are independent of λ − κ. The factorial uses Boston Blue Clay and London Clay; kaolin and Lower Cromer Till complete the calibration record. The critical-state ratio M enters the analysis only through the normalization of Section 2.1 and the recovery of physical quantities (η = Mη̅, D = MD̅).

Soil	λ	κ	ν	M	e0	Ω	Ψ	n policy
Boston Blue Clay	0.184	0.036	0.10	1.353	2.059	0.93	1.4	n = 0.5+2/OCR (fixed at initial OCR)
Kaolin Clay	0.26	0.05	0.20	0.896	2.957	1.17	1.5	monotonic decreasing (n ≤ 4)
London Clay	0.168	0.064	0.25	0.827	1.843	1.23	1.4	n = 0.5+2/OCR (fixed at initial OCR)
Lower Cromer Till	0.063	0.018	0.30	1.20	0.747	0.97	1.05	monotonic decreasing

Table 2. Locked Level-I parameter set of the companion study [21] (single set, whole program; identical to its published Supplementary S2; soil sources as compiled in [21,22]). The λ and κ listed here are the source-calibration constants of [22]; the integrator takes their ratio κ/λ and applies it to the AJOP-consistent compression index of [21], giving λ = 0.187, κ = 0.0365 for Boston Blue Clay and λ = 0.0586, κ = 0.0223 for London Clay. The normalized constant-λ indices are independent of λ − κ. The factorial uses Boston Blue Clay and London Clay; kaolin and Lower Cromer Till complete the calibration record. The critical-state ratio M enters the analysis only through the normalization of Section 2.1 and the recovery of physical quantities (η = Mη̅, D = MD̅).

Soil	λ	κ	ν	M	e0	Ω	Ψ	n policy
Boston Blue Clay	0.184	0.036	0.10	1.353	2.059	0.93	1.4	n = 0.5+2/OCR (fixed at initial OCR)
Kaolin Clay	0.26	0.05	0.20	0.896	2.957	1.17	1.5	monotonic decreasing (n ≤ 4)
London Clay	0.168	0.064	0.25	0.827	1.843	1.23	1.4	n = 0.5+2/OCR (fixed at initial OCR)
Lower Cromer Till	0.063	0.018	0.30	1.20	0.747	0.97	1.05	monotonic decreasing

The last column of Table 2 states each soil’s hardening-exponent policy in words; Table 3 gives the numerical exponent n it produces, each entry written as n (OCR) — the exponent value with its overconsolidation ratio in parentheses — at the ratios used in the analysis. For Boston Blue Clay and London Clay n follows n = 0.5 + 2/OCR, fixed at the initial OCR (convention B); for kaolin and Lower Cromer Till the exponent is instead held monotonically decreasing with OCR under the per-soil cap n ≤ 4.

Table 3. Hardening exponent n used in the analysis (locked set), each entry n (OCR). Boston Blue Clay and London Clay follow the n–OCR convention of [22] at fixed initial OCR (convention B); kaolin and Lower Cromer Till are monotonically decreasing under the per-soil cap n ≤ 4.

Soil	n (OCR)
Boston Blue Clay	2.5(1.0)	1.5(2.0)	1.0(4.0)	0.74(8.5)
London Clay	2.5(1.0)	1.5(2.0)	0.83(6.0)	0.61(18.7)
Kaolin Clay	4.0(1.0)†	4.0(1.2)	4.0(1.4)	4.0(1.8)	3.1(2.5)	3.1(3.0)	2.25(4.5)	2.1(6.4)
Lower Cromer Till	4.0(1.0)†	4.0(2.0)	0.45(4.1)	0.45(9.7)

† At OCR = 1 the exponent is unidentifiable (the spacing ratio $ℛ$ → 1); the per-soil cap value is adopted.

Table 3. Hardening exponent n used in the analysis (locked set), each entry n (OCR). Boston Blue Clay and London Clay follow the n–OCR convention of [22] at fixed initial OCR (convention B); kaolin and Lower Cromer Till are monotonically decreasing under the per-soil cap n ≤ 4.

Soil	n (OCR)
Boston Blue Clay	2.5(1.0)	1.5(2.0)	1.0(4.0)	0.74(8.5)
London Clay	2.5(1.0)	1.5(2.0)	0.83(6.0)	0.61(18.7)
Kaolin Clay	4.0(1.0)†	4.0(1.2)	4.0(1.4)	4.0(1.8)	3.1(2.5)	3.1(3.0)	2.25(4.5)	2.1(6.4)
Lower Cromer Till	4.0(1.0)†	4.0(2.0)	0.45(4.1)	0.45(9.7)

† At OCR = 1 the exponent is unidentifiable (the spacing ratio $ℛ$ → 1); the per-soil cap value is adopted.

4. Results: Collapse and Migration of the Flow Effect

4.1. The Flow Main Effect Collapses

Table 4 reports the paired factorial at OCR = 5, with the GL-window column reproducing the companion study [21]. Replacing the GL-window flow by the critical-state Caputo flow leaves every flow-independent index unchanged to the reported precision (ΔG, ΔC, I_GC — an internal consistency check of the pairing) and collapses the flow main effect from −60.25% to +0.69% of |R(C0)| for Boston Blue Clay and from −60.32% to +0.10% for London Clay; the flow×compression interaction collapses equally (+41.70% → −0.60%; London: undefined → +0.13%). All tabulated values use n = 30,000 increments (the sweep of record); at n = 60,000 the Boston Blue Clay value moves only from +0.689% to +0.679%, and the near-zero London flow×compression contrast from +0.13% to +0.06%, i.e., within discretisation noise of contrasts that are numerically zero. Across the full OCR sweep (Table 5, Figure 4) the collapsed ΔF never exceeds 1.0% for either soil, against −37% to −115% for the GL window. The result is path-robust: under the approximately undrained path at OCR = 5, ΔF_cap = 0.001% (BBC) and 0.000% (London). It is also robust to the assumed fractional-order law: replacing α(OCR) by constant α = 0.55 or 0.75 across a 12-point grid (both soils, OCR = 3, 5, 10) leaves the collapse intact everywhere, with |ΔF_cap| ≤ 2.5% and the largest residual (+2.47%, Boston Blue Clay, OCR = 3, α = 0.55) occurring exactly where the mechanism predicts — the most strongly non-local operator on the longest approach path — and falling below 0.3% at α = 0.75.

4.2. The Collapse Is a Critical-State Dwell Effect

Figure 3 resolves the mechanism. At small imposed shear budgets — before the state reaches the critical state — the Caputo and integer flows differ strongly and ΔF_cap is large (+71% at budget 0.2 for BBC); as the budget grows and the state dwells at the critical state, where the Caputo flow is exactly associated and shares the integer terminal, ΔF_cap decays monotonically toward zero (+3.9% at budget 3, +0.7% at budget 6). London Clay, whose states reach the critical state faster, decays faster (below 1% already at budget 1.5). The companion study [21] showed that its classical configuration spends 96–100% of the imposed budget traversing the band between the two flow rules’ critical-state points, while its teardrop configuration reaches its own critical state early and remains there for the rest of the budget; in both cases the path dwells at or beside the critical state, where the GL window retains finite non-associativity and therefore accumulates its −60%, while the critical-state-anchored operator cannot. The dominant fractional flow main effect of the companion study is thereby identified as a critical-state-displacement effect of the window semantics — exactly the kind of operator-specific property that study cautioned about [21], now quantified.

4.3. The Flow Effect Migrates into the Geometry×Flow Interaction

The flow operator does not become irrelevant. The least interpretation-laden statement uses simple effects: at OCR = 5 for Boston Blue Clay, the simple effect of replacing the integer by the Caputo-CS flow is +0.69% on the classical geometry but −38.8% on the teardrop (London: +0.10% versus −10.7%), whereas under the GL window it is large on both geometries (−60.2% and −95.5% for BBC; −60.3% and −66.9% for London). Under the critical-state operator the flow effect is therefore strongly geometry-dependent — on the teardrop the Caputo flow (critical state at $\overline{p}$* = 0.490) and the embedded rule of [22] (crossing at 0.341) disagree about where dilation stops, while on the classical reference the two flows share the terminal and nearly coincide — and in factorial bookkeeping this expresses itself as a dominant geometry×flow interaction, I_GF = −39.5% (BBC) and −10.8% (London), of the same magnitude as and opposite sign to ΔG (Figure 5). We use ‘migration’ only in this bookkeeping sense: main effects and interactions depend on the coding and reference level of the design, and the invariant content is the geometry-conditional simple effects just quoted. The three-way interaction rises accordingly (I_GFC up to +18% for BBC). A practical corollary concerns the hardening singularity of the companion study [21]: because the Caputo-flow configurations follow different stress paths, London Clay AJOP-compression configurations remain regular through OCR = 6 where their GL counterparts were already undefined, and the singular exposure re-enters only at OCR ≥ 7 through the compression factor itself (Table 5).

Figure 4. Paired factorial across OCR = 3–10: ΔF collapses under the Caputo-CS operator (dashed) for both soils while ΔG is flow-independent by construction; the geometry×flow interaction absorbs the operator difference.

Figure 4. Paired factorial across OCR = 3–10: ΔF collapses under the Caputo-CS operator (dashed) for both soils while ΔG is flow-independent by construction; the geometry×flow interaction absorbs the operator difference.

Figure 5. Where the flow effect lives (OCR = 5): under the GL window it is a main effect; under the critical-state Caputo operator it migrates into the geometry×flow interaction.

Figure 5. Where the flow effect lives (OCR = 5): under the GL window it is a main effect; under the critical-state Caputo operator it migrates into the geometry×flow interaction.

5. Two Routes to State-Dependent Non-Associativity: A Clean Discriminant

State-dependent non-associativity has classically been introduced through explicit state variables — the state parameter of Been and Jefferies [33], the state-dependent dilatancy of Li and Dafalias [34] and its two-surface and bounding-surface realisations [35,36] — and, in overconsolidated-clay modelling, through interpolation constructions of the MIT type [37]. The program now contains two distinct mechanisms that generate state-dependent non-associated flow on the same surface family without an explicit state variable: the interpolation-exponent route of [22], in which the critical stress ratio is scaled as M_k = M·$ℛ$⁻ⁿ with n tied to overconsolidation, and the critical-state Caputo route [9], in which non-locality toward a fixed critical-state terminal modulates the flow direction. The results above yield a discriminant that is independent of calibration: the Caputo route holds the phase-transformation point fixed at $\overline{p}$* = e^−1/Ψ for every α and every OCR and modulates only the approach to it, whereas the interpolation route moves the phase-transformation point with state. For the Boston Blue Clay shape these anchor points are far apart (0.490 versus 0.341, Figure 1), so undrained effective-stress paths of heavily overconsolidated samples — whose contractive-to-dilative transition is expected, by analogy with the phase-transformation behaviour long resolved experimentally in sands since Ishihara et al. [38] and with the steady/critical state established there as a robust reference [39] — should distinguish the routes directly. A comparably systematic phase-transformation database for heavily overconsolidated clays is sparser, which we flag as a limitation of the proposed discriminant. A quantitative equivalence fit (whether an effective n(OCR) exists that reproduces the Caputo dilatancy within data scatter) is deliberately left to future work, since it must be conducted against laboratory data at the full-model level rather than at the component level isolated here.

6. Limitations

The scope restrictions of the companion study [21] carry over unchanged and are restated rather than resolved: one idealized shear-strain-controlled constant-p′ path with an approximate undrained variant; the α(OCR) form remains an uncalibrated modeling assumption (the collapse results were obtained under the same assumption as the companion for comparability, and were additionally checked at constant α = 0.55 and 0.75, Section 4.1, supporting a structural rather than α-specific mechanism); the compression law retains its calibration-domain idealisation; and the critical-state Caputo flow rule itself has not been validated against laboratory data — the present paper establishes its exact construction, its limits, and its factorial consequences, not its empirical superiority. The log-conjugate convention for the dilatancy magnitude is one of the admissible Jacobian conventions; the location of the zero, the associativity at the critical state, and the collapse results are convention-independent under the stated endpoint, loading, and convergence conditions, but reported magnitudes of D are not.

7. Conclusions

First, the teardrop bounding surface admits an exact critical-state Caputo fractional gradient in closed form once the operator is defined with respect to the log-stress function ψ = ln(1/$\overline{p}$) [24]: the critical-state terminal is exactly t* = 1/Ψ ($\overline{p}$* = e^−1/Ψ, independent of Ω), and the gradient reduces to incomplete-Beta–Kummer ₁F₁ (tip terminal) and Humbert Φ₁ (critical-state terminal, both sides) expressions, verified against singularity-aware quadrature over 1,240 cases to relative errors below 10⁻¹¹ (worst case 5.1×10⁻¹² on the far dry side; tip-terminal and shifted-terminal cases below 10⁻¹⁵); the resulting flow rule recovers integer associated flow (in the log-stress conjugate representation) as α → 1 and its flow vector coincides with the associated (purely deviatoric) vector exactly at the critical state for every fractional order, with a sign-changing two-sided zero.

Second, substituting this operator for the fixed-window GL flow collapses the dominant flow main effect of the companion factorial [21] from −60% to below 1% of the baseline response for both soils included in the factorial analysis at every overconsolidation ratio tested, under drained and approximately undrained paths alike, and the collapse persists across a 12-point constant-α grid (both soils; OCR = 3, 5, 10; α = 0.55, 0.75) with |ΔF| ≤ 2.5%, and under the $\overline{p}$-conjugate flow convention (|ΔF| < 0.1% at OCR = 5), supporting a structural, convention-independent mechanism: under the loading and convergence conditions verified here, an endpoint identity of the hardening law makes the converged response a function of the terminal state only, and the constant-p′ path makes the shared normalized terminal fix the absolute one.

Third, the collapse is a critical-state dwell effect: at small shear budgets the two operators differ strongly (up to +71% of baseline) and the difference decays monotonically as the state reaches and dwells at the shared critical state — identifying the flow main effect of [21] as a critical-state-displacement effect of the window semantics rather than an operator-independent consequence of fractional non-associativity; the reported −60% is not transferable across window semantics.

Fourth, under the critical-state operator the flow effect is strongly geometry-conditional (simple effects +0.7% on the classical surface versus −38.8% on the teardrop for Boston Blue Clay at OCR = 5), appearing in factorial bookkeeping as a dominant geometry×flow interaction; a numerical corollary under the present model assumptions is that the singular exposure of the embedded compression law is itself flow-dependent, London Clay AJOP-compression configurations remaining regular through OCR = 6 under the Caputo flow where their GL counterparts were undefined.

Fifth, the two available routes to state-dependent non-associativity on this surface family are cleanly distinguishable by a calibration-independent discriminant: the interpolation-exponent route of [22] moves the phase-transformation point with state, while the critical-state Caputo route holds it fixed at e^−1/Ψ (0.341 versus 0.490 for the Boston Blue Clay shape), making the discriminant experimentally accessible, in principle, on undrained paths of heavily overconsolidated samples [38,39]. All results inherit the reproducibility standard of the companion study: single locked parameter set, no-clamp singularity semantics, self-testing reference implementation, and full scripts as supplementary material.