Defining Function Dependence of Riemann–Liouville Plastic Flow on a Common Convex Yield Surface

1. Introduction

A yield boundary describes where plastic loading can occur. A flow rule describes the direction of plastic deformation. Classical associated flow follows the boundary normal. Non-associated flow requires additional directional information [1,2,3]. We use the standard level-set description of a smooth boundary [4]. A scalar defining function f is zero on that boundary and negative on the elastic side. Multiplication by a positive function h gives another description, hf, of the same region. On f = 0, the ordinary gradient is multiplied by h, so its normalized direction is unchanged. The same geometry therefore admits many equivalent scalar descriptions. We ask whether this equivalence survives when a specified nonlocal derivative replaces ordinary differentiation in stress space.

Stress-fractional plasticity can produce non-associated directions without a separate conventional plastic potential. Sumelka obtained a fractional construction that recovers a classical non-associated Drucker–Prager special case [5]. Sun and Shen introduced a fractional flow rule into a granular bounding-surface model [6]. Sun, Gao and Zhu developed state-dependent stress–dilatancy modeling without a plastic potential [7]. Each construction specifies a constitutive model. The boundary alone does not specify the fractional operator. Qu, Sun and Sumelka emphasize the stress length scale and distinguish operators that sample past stress from those directed toward a reference critical state [8]. A teardrop-surface study also found different factorial flow contrasts after replacing a fixed-window operator with a critical-state-directed Caputo construction [9]. Coordinates and terminal conventions therefore need to be stated together with fractional order.

The role of the complete model is clear in applications. Cyclic granular formulations require rules for loading reversal as well as a fractional direction [10]. Xu, Shen and Sun include material-state and fabric effects in a cyclic soil–structure interface model [11]. Wu, Cheng and Zhu combine fractional flow, subloading and multiple hardening rules, with comparisons against cyclic Karlsruhe sand tests [12]. Li and co-workers compare a temperature-dependent formulation with drained Ohya rock measurements [13]. A controlled component-interaction study further shows that geometry, fractional flow and compression hardening interact under its stated loading paths [14]. These results concern complete constitutive systems. They leave open whether a fixed fractional operator gives the same direction for different scalar descriptions of one boundary.

Three-dimensional models require further choices. Lu and co-workers define soil flow in a transformed stress space [15]. Sumelka and Nowak study numerical approximation of a fractional plastic flow rule [16]. Tai and Sun examine plastic dissipation for a modified Cam-clay construction [17]. Zhang and co-workers combine a fractional flow rule with thermodynamic structure and learned state-dependent parameters [18]. These studies motivate separate checks of flow direction, conversion to physical tensors, hardening consistency and admissibility. A successful check in one category does not establish the others.

The present question concerns differentiation in stress space. Time-fractional rheology addresses a different problem. Caputo places fractional differentiation in a stress–strain relation [19]. Standard treatments also distinguish Riemann–Liouville and Caputo operators through their endpoint terms and initial data [20,21,22]. Suzuki, Naghibolhosseini and Zayernouri develop return mapping for visco-elasto-plastic elements with time-dependent memory [23]. Feischl, Niederkofler and Wohlmuth study discretizations of a finite-support Riesz–Caputo elasto-plastic problem [24]. These operators and evolution problems differ from the one-sided stress-space comparison studied here.

The mathematical literature likewise requires a precise operator definition. Ortigueira and Machado discuss criteria for fractional derivatives [25]. Valério, Ortigueira and Lopes examine different definitions and their domains [26]. Tarasov and Cresson and Szafrańska explain why ordinary product and chain rules cannot simply be carried over to non-integer differentiation [27,28]. Tarasov also discusses nonlocality as a defining issue [29]. General fractional vector calculus and its Riesz-form extension use distinct multidimensional constructions [30,31]. Fernandez and Al-Refai establish properties for specified non-singular-kernel operators [32]. These studies motivate careful declarations of kernels, domains and coordinates. They do not establish the boundary-equivalence result below. A separate Caputo restart study concerns preservation of inherited history under reversals [33]. In our comparison, the lower terminals remain fixed. Only the scalar defining function changes.

The related Caputo study addresses terminal compatibility and admissibility [34]. The companion AJOP-T formulation provides a tensor map and complete hardening consistency [35]. That formulation is available as a preprint and has not been peer-reviewed. Our question comes before its mechanical implementation. With all operator settings fixed, does one zero boundary determine one normalized fractional vector? Nonlocality suggests that values away from the boundary matter. The stronger question is whether strict convexity, a closed boundary and agreement of any prescribed finite number of boundary derivatives can remove that dependence.

We show that these conditions do not restore equivalence for the specified componentwise left Riemann–Liouville operator. For each finite derivative order, we construct globally strictly convex polynomials with a common elliptical boundary. Their boundary derivatives agree through that order. A determinant with a strict sign proves that their fractional directions differ at a stated point for every fractional order between zero and one. For each fixed fractional order, continuity extends the difference to a nearby boundary arc. Exact formulas, independent quadrature and a convolution based only on function values check the angular separations. We then test the operator on a published modified Cam-clay geometry and examine a coupled mechanical diagnostic.

2. Fixed Operator and Equivalence Requirement

The Riemann–Liouville (RL) operator contains an endpoint term as well as a Caputo integral [20,21,22,36]. We retain both terms. Removing the endpoint term would change the operator, even when the sampled function is constant. Holding this convention fixed is essential: the comparison must change the defining function without also changing the derivative.

We use dimensionless coordinates (x,y) on $\mathcal{D}=[0,\mathrm{\infty }{)}^2$. For x,y > 0 and 0 < α < 1, define the left RL component by the following endpoint-plus-Caputo identity:

$${\mathcal{G}}_{\alpha ,x}\left[f\right]\left(x,y\right)={\displaystyle \frac{f\left(0,y\right)x^{-\alpha }+{\int }_0^x{\partial }_{s}f\left(s,y\right){\left(x-s\right)}^{-\alpha }\hspace{0.17em}\mathit{ds}}{\Gamma \left(1-\alpha \right)}},$$

(1)

Define the y component in the same way. Each component differentiates one coordinate while holding the other fixed. Both lower terminals are zero. Coordinates, scales and internal variables remain fixed throughout the comparison. The functions below are polynomials, and the weakly singular integrals are finite. Write $g={\mathcal{G}}_{\alpha }\left[f\right]$ for the fractional vector and u = g/∥g∥ for its unit direction, defined whenever g ≠ 0.

A lower terminal specifies where a coordinate integration path starts. It is part of the operator definition. For example, a constant slice with value C has the RL derivative ${\mathcal{G}}_{\alpha ,x}\left[C\right]=C{x}^{-\alpha }/\Gamma \left(1-\alpha \right)$. The value need not be zero. This endpoint contribution helps explain why agreement only at the yield boundary is insufficient.

Let $\widehat{f}=\mathit{hf}$, with h > 0 throughout the stated domain. The two functions then have the same zero boundary. They also have the same sign at every point, so they describe the same elastic side. The ordinary product rule gives, on f = 0,

$$\nabla \widehat{f}=h\nabla f.$$

(2)

The ordinary normalized normals therefore agree. Fractional direction equivalence is a separate condition. It requires ${\mathcal{G}}_{\alpha }\left[\widehat{f}\right]=c{\mathcal{G}}_{\alpha }\left[f\right]$, where the scalar c is positive at the evaluation point. This condition means that the two fractional vectors have the same orientation. A common zero set does not by itself establish it.

A small change in the defining function provides a useful diagnostic. Set h = 1 + ϵϕ and $b={\mathcal{G}}_{\alpha }\left[\varphi f\right]$. Linearity of the RL operator and ordinary differentiation of the unit vector give

$${\left.{\displaystyle \frac{d}{dϵ}}{\displaystyle \frac{g+ϵb}{\parallel g+ϵb\parallel }}\right|}_{ϵ=0}={\displaystyle \frac{\left(I-u{u}^T\right)b}{\parallel g\parallel }}.$$

(3)

A positive constant multiplier leaves the normalized direction unchanged. A perturbation with a nonzero component perpendicular to g changes that direction. Equation (3) differentiates the normalization map in ordinary calculus. It does not assume a fractional product rule.

Example 1. For f = x + y − 1 and h = 1 + βx, evaluated on x + y = 1, direct calculation gives

$${\mathcal{G}}_{\alpha ,x}\left[\mathit{hf}\right]=\left(1+{\displaystyle \frac{\beta x}{2-\alpha }}\right)g_x, {\mathcal{G}}_{\alpha ,y}\left[\mathit{hf}\right]=\left(1+\beta x\right)g_y.$$

(4)

At x = y = 1/2, α = 1/2 and β = 1, the angle is exactly arctan(1/17) = 3.366460663^∘. This calculation already gives a counterexample. It does not depend on quadrature, a loading path or the mechanical model used later.

Figure 1 shows why Example 1 works. Positive multiplication leaves the ordinary boundary normal unchanged. In contrast, the RL operator samples the full defining function along two coordinate paths. The functions differ along those paths, even though both vanish at the endpoint. Equation (4) therefore gives different unit directions. The next examples also make higher derivatives agree on the boundary.

3. Exact Counterexamples

3.1. An affine Boundary

Endpoint identity. Fix a positive endpoint distance d and an integer n ≥ 1. The lower-terminal RL derivative of (s − d)ⁿ, evaluated at s = d, is

$$D_{0+}^{\alpha }{\left(s-d\right)}^n{|}_{s=d}={\displaystyle \frac{{\left(-1\right)}^{n-1}\alpha }{\left(n-\alpha \right)\Gamma \left(1-\alpha \right)}}{d}^{n-\alpha }.$$

(5)

Substitute v(s) = (s − d)ⁿ into (1). The endpoint contribution is (−1)ⁿd^n−α, and the integral contribution is n(−1)ⁿ⁻¹d^n−α/(n − α). Their sum, divided by Γ(1 − α), gives (5). The order of operations matters. We first differentiate with the other coordinate fixed, and then evaluate on the boundary. Setting f to zero before differentiation would replace the sampled function by a different one.

The polynomial evaluations below reduce to powers and beta integrals. We use the beta–gamma identity given in the NIST Digital Library of Mathematical Functions [37]. The endpoint identity above follows directly from (1).

Proposition 1. Fix 0 < α < 1 and any finite nonnegative integer K. There exist convex polynomial defining functions on $\mathcal{D}$ with the same affine boundary x + y = 1 and the same sign regions. Their boundary derivatives agree through order K, but their normalized RL directions differ at boundary points whose coordinates are positive and unequal. The boundary is a line segment in $\mathcal{D}$. It is not a closed curved surface.

Proof. Choose an even integer n ≥ 4 with n > K, take 0 < c < 1/n, and set

$$f=x+y-1, \widehat{f}=f+c{f}^n=\left(1+c{f}^{n-1}\right)f.$$

(6)

On $\mathcal{D}$, we have f ≥ −1. Hence 1 + cfⁿ⁻¹ ≥ 1 − c > 0 and 1 + cnfⁿ⁻¹ ≥ 1 − cn > 0. The positive multiplier preserves the boundary and the elastic side. Both functions are convex. The base function f is affine. For the other function, with a = (1,1)^T, the Hessian is

$${\nabla }^2\widehat{f}=\mathit{cn}\left(n-1\right)f^{n-2}a{a}^T\succcurlyeq 0.$$

(7)

Every derivative of cfⁿ with total order below n vanishes on f = 0. Thus the required finite collection of boundary derivatives agrees. Now take (x,y) with x + y = 1 and x,y > 0. Each coordinate slice has the form (s − d_i), with d_x = x and d_y = y. Divide (5) by its n = 1 case to obtain

$${\displaystyle \frac{{\mathcal{G}}_{\alpha ,i}\left[\widehat{f}\right]}{{\mathcal{G}}_{\alpha ,i}\left[f\right]}}=1-c{\displaystyle \frac{1-\alpha }{n-\alpha }}{d}_i^{n-1}, i=x,y.$$

(8)

The base components and both multiplying factors are positive. If x ≠ y, the factors differ. The two vectors therefore cannot be collinear, which proves the proposition. □ For each prescribed finite K, the construction chooses a suitable polynomial power. It does not produce one nonidentical analytic pair whose derivatives agree at every order.

For n = 4, c = 1/8, α = 1/2 and (x,y) = (1/3,2/3), (8) gives factors 1511/1512 and 1504/1512. Hence

$$\Delta \theta =\arctan \left({\displaystyle \frac{7\sqrt{2}}{4519}}\right)={0.125514134}^{\circ }.$$

(9)

This pair preserves convexity and all boundary derivatives through order three. The Hessian is the matrix of second partial derivatives. Here it is positive semidefinite, which establishes convexity but not strict convexity [38]. The result is restricted to $\mathcal{D}$. No claim is made about additional zeros outside that domain.

3.2. A Closed Strictly Convex Boundary

Proposition 2. Fix positive constants a,b and a finite nonnegative integer K. There exist globally strictly convex polynomial functions on ${\mathbb{R}}^2$ with the same bounded convex sublevel set and the same closed elliptical zero boundary. Their partial derivatives agree through order K everywhere on that boundary. For every 0 < α < 1, their normalized componentwise left RL directions with zero lower terminals differ at (3a,2b). For each fixed fractional order, they also differ on a neighboring boundary arc.

Proof. Choose an even n ≥ 4 with n > K, and 0 < c < 1/n. Define

$$F={\left({\displaystyle \frac{x-2a}{a}}\right)}^2+{\left({\displaystyle \frac{y-2b}{b}}\right)}^2-1, \widehat{F}=F+c{F}^n.$$

(10)

The base function satisfies F ≥ −1 on all of ${\mathbb{R}}^2$. The inequalities 1 + cFⁿ⁻¹ ≥ 1 − c > 0 and 1 + cnFⁿ⁻¹ ≥ 1 − cn > 0 preserve the zero set and sign regions. They also give the following positive-definite Hessian:

$${\nabla }^2\widehat{F}=\left(1+\mathit{cn}{F}^{n-1}\right){\nabla }^{2}F+\mathit{cn}\left(n-1\right)F^{n-2}\nabla F\nabla F^T\succ 0.$$

(11)

To see the strict sign, note that ∇²F = 2diag(a⁻²,b⁻²)≻0. In (11), the first term is therefore positive definite and the second is positive semidefinite. The standard Hessian criterion implies global strict convexity [38]. Both functions describe the same elastic region, the filled ellipse F ≤ 0, which lies in the positive quadrant. Every partial derivative of cFⁿ with total order below n vanishes on F = 0. This proves the required boundary matching.

The direction calculation uses a second identity. Let the polynomial v vanish to at least second order at a positive endpoint d. This assumption makes the endpoint terms integrable after integration by parts in (1), giving

$$D_{0+}^{\alpha }v\left(d\right)=-{\displaystyle \frac{\alpha }{\Gamma \left(1-\alpha \right)}}{\int }_0^d{\displaystyle \frac{v\left(s\right)}{{\left(d-s\right)}^{1+\alpha }}}\hspace{0.17em}\mathit{ds}.$$

(12)

The boundary contribution at d is zero. The contribution at the lower terminal cancels the explicit RL endpoint term. At z_* = (3a,2b), let $g={\mathcal{G}}_{\alpha }\left[F\right]\left(z_{*}\right)$. Direct polynomial differentiation then gives

$$\begin{array}{r}{g}_x={\displaystyle \frac{\alpha a^{-\alpha }{3}^{1-\alpha }\left(1+\alpha \right)}{\left(1-\alpha \right)\left(2-\alpha \right)\Gamma \left(1-\alpha \right)}}>0\\ g_y=-{\displaystyle \frac{\alpha b^{-\alpha }{2}^{2-\alpha }}{\left(2-\alpha \right)\Gamma \left(1-\alpha \right)}}<0.\end{array}$$

(13)

Write $w={\mathcal{G}}_{\alpha }\left[F^n\right]\left(z_{*}\right)$ for the vector contributed by the added power. Its coordinate slices are [(s/a − 3)(s/a − 1)]ⁿ and (s/b − 2)²ⁿ. Both are nonnegative and are positive on sets of nonzero length. Applying (12) with the corresponding endpoint substitutions gives

$$\begin{array}{r}{w}_x=-{\displaystyle \frac{\alpha a^{-\alpha }}{\Gamma \left(1-\alpha \right)}}{\int }_0^3{u}^{n-\alpha -1}{\left(u-2\right)}^n\hspace{0.17em}\mathit{du}<0\\ w_y=-{\displaystyle \frac{\alpha b^{-\alpha }{2}^{2n-\alpha }}{\left(2n-\alpha \right)\Gamma \left(1-\alpha \right)}}<0.\end{array}$$

(14)

Consequently,

$$\det \left(g,{\mathcal{G}}_{\alpha }\left[\widehat{F}\right]\right)=c\left(g_x{w}_y-g_y{w}_x\right)<0.$$

(15)

Equation (13) gives a positive first base component and a negative second component. Both added components in (14) are negative. The two terms in (15) therefore make the determinant strictly negative. The vectors are nonzero and noncollinear for every allowed α. Their RL components are continuous for x,y > 0. For each fixed fractional order, the determinant remains nonzero on a sufficiently small boundary arc around z_*. This proves the proposition. □ The same polynomial pair works for every α ∈ (0,1). A larger prescribed K may require a different polynomial power.

For the circle a = b = 1, n = 4, c = 1/8, and α = 1/2, the two vectors at (3,2) are $g=\left(\sqrt{3},-2\sqrt{2}/3\right)/\sqrt{\pi }$ and $\widehat{g}=\left(34631\sqrt{3}/40040,-26\sqrt{2}/15\right)/\sqrt{\pi }$. Hence the exact separation is

$$\Delta \theta =\arctan \left({\displaystyle \frac{416838\sqrt{6}}{1767869}}\right)={30.008768325}^{\circ }.$$

(16)

Proposition 2 replaces the flat boundary by a smooth closed one and retains a bounded, strictly convex elastic region. The result applies to the stated ellipse family and RL operator. It does not assert the same theorem for every closed surface or every terminal choice. Nor does it cover all modified Cam-clay, teardrop or Riesz–Caputo formulations. The ellipse is a mathematical counterexample, rather than a calibrated model of a particular material.

Figure 2 explains the integral mechanism for the half-order circle example in (16). Let F_x(s) = F(s,2) and F_y(s) = F(3,s), with endpoint distances d_x = 3 and d_y = 2. The slices approach the same boundary values. Away from the endpoint, however, the difference cF_i(s)⁴ remains. The substitution s = d_it gives the nonnegative integrands $H_i\left(t\right)=c\hspace{0.17em}{d}_i{F}_i{\left(d_{i}t\right)}^4/{\left[d_i\left(1-t\right)\right]}^{1+\alpha }$. These include the Jacobian factor d_i and have continuous zero limits at t = 1. Equation (12) gives a negative correction to each RL component. Together with the opposite signs in (13), these corrections give the strict determinant sign in (15).

3.3. Cubic Comparison and Complete Constitutive Equivalence

The globally sign-preserving cubic map used in the host is $\widehat{f}=f+c{f}^3$, c > 0. Its surface gradient and Hessian equal those of f. Formula (5) gives

$${\displaystyle \frac{{\mathcal{G}}_{\alpha ,i}\left[f+c{f}^3\right]}{{\mathcal{G}}_{\alpha ,i}\left[f\right]}}=1+c{\displaystyle \frac{1-\alpha }{3-\alpha }}{d}_i^2.$$

(17)

For c = 1/2, α = 1/2 and (x,y) = (1/3,2/3), the factors are 91/90 and 94/90. The resulting angle is $\Delta \theta =\arctan \left(\sqrt{2}/93\right)={0.871206706}^{\circ }$. The cubic function is not convex throughout the negative-f half-space. The quartic example on $\mathcal{D}$ and the cubic example therefore serve different purposes. The quartic proves dependence within a convex class. The cubic matches the relabeling used in the mechanical diagnostic.

Figure 3 compares the symmetric quartic and cubic affine examples. Their directions agree at the midpoint because the two component factors coincide there. The nonsymmetric multiplier in Figure 1 does not have that cancellation. Agreement at one boundary point therefore does not imply agreement along the boundary.

The failure of an ordinary product rule does not by itself prove Proposition 2. That proposition also requires global strict convexity, a closed boundary and finite-order derivative matching. These boundary conditions leave the nonlocal slice integrals in (14) nonzero. This is the additional content beyond the general algebraic restrictions discussed in [27,28] and the nonlocality principle in [29].

For the affine family f + cfⁿ, the condition x = y makes the two component factors equal. This cancellation belongs to that symmetric family. The multiplier h = 1 + x in (4) gives a nonzero angle even at x = y = 1/2. The meridional coordinates used here are also different from Lode coordinates. Equality of these two coordinates has no general implication for a Lode-axis state.

A flow direction is only one part of a constitutive system. The following calculation uses the standard active-loading equations of rate-independent plasticity [2,3]. It shows which quantities must change together when the defining function changes.

Full-system scaling. Let z denote the internal state, m the plastic direction and a the internal evolution per unit plastic multiplier. During active plastic loading, a rate-independent system satisfies

$$\dot{\sigma }=C_e\left(\dot{\epsilon }-\dot{\gamma }m\right), \dot{z}=\dot{\gamma }a, f_{\sigma }:\dot{\sigma }+f_z\cdot \dot{z}=0.$$

(18)

Suppose the defining function changes according to $\widehat{f}=\mathit{hf}$. On the surface, also impose $\left(\widehat{m},\widehat{a}\right)=c\left(m,a\right)$, with h,c > 0. The consistency numerator then scales by h, and the denominator scales by hc. Consequently, $\dot{\widehat{\gamma }}=\dot{\gamma }/c$. The multiplier change cancels the change in the complete evolution direction. Stress and internal-variable rates remain unchanged wherever the active solve is nonsingular. This is a sufficient condition for equivalence. It is not a necessity theorem for every possible path. In particular, the same factor must multiply the plastic direction and every internal evolution component. Ordinary hardening consistency alone does not guarantee this relation for two distinct fractional directions.

4. Numerical Verification

4.1. Exact Formula and Independent Discretization Tests

Supplementary Material S1 v14 records the complete v11 computation, which started from empty result directories. It includes the polynomial, closed-boundary, GL, modified Cam-clay, mechanical and refinement checks. Later revisions change the presentation and references. They do not substitute historical v02–v04 mechanical outputs or claim another full campaign run. Source hashes, script logs, timings and output hashes identify the recorded computation. The publication-artwork script reads those CSV results. A separate script generates the explanatory panels. PROVENANCE.csv links each numerical claim and figure to its script and output. The numerical tests check the implementation and reported magnitudes. The propositions supply the mathematical proof.

Figure 2 explains the mechanism at one point. Figure 4 adds a scan around the common circle. We use the boundary parameter (x,y) = (2 + cost,2 + sint) and 721 samples for each displayed order. Both fractional vectors have nonzero norm at every sampled point. The scan illustrates the theorem near the specified point. A finite scan cannot prove that the vectors are nonzero at every point of every closed surface.

We evaluated (5) in 240 combinations: powers 1–12, five fractional orders and four endpoint distances. A 70-digit coefficient calculation and independent adaptive quadrature agree with the formula to scaled errors below 6 × 10⁻¹⁶. The closed-ellipse test has 30 combinations, using five even powers from 4 to 12 and six orders from 0.1 to 0.99. Two 85-digit polynomial calculations agree to scaled error below 5 × 10⁻⁷⁰. Independent adaptive quadrature agrees to below 5 × 10⁻¹⁵. A common positive gamma factor is removed in the high-precision comparison. These are arithmetic checks, not physical error bounds. The two evaluation routes for the circle angle differ by about 4 × 10⁻¹¹ degrees. One uses high-precision coefficients, and the other uses double-precision monomial/Gamma calculations. Small differences between arithmetic routes are expected in finite precision [39]. Both values round to 30.008768 degrees.

We use three independent forms of calculation. Polynomial evaluation applies the RL formula term by term. Adaptive quadrature evaluates the endpoint-plus-integral identity. GL convolution samples function values on successively finer grids. QUADPACK provides established adaptive integration procedures [40]. Convolution consistency [41] and practical fractional-equation algorithms [42] motivate the separate discretization check. Diethelm, Garrappa and Stynes discuss why numerical accuracy must be tested for the actual fractional problem [43]. Predictor–corrector methods [44], time-fractional diffusion schemes [45,46] and Caputo approximations [36] address different problems. Their error estimates are not error bounds for the present stress-slice calculation.

The Grünwald–Letnikov (GL) approximation uses weighted function values, a standard approach in fractional finite differences [47]. Here it provides a check that uses neither the polynomial RL evaluation nor derivative quadrature. For a coordinate slice, we compute

$$G_N\left[v\right]\left(d\right)=h^{-\alpha }{\displaystyle \sum _{j=0}^N}{w}_{j}v\left(d-\mathit{jh}\right), h=d/N, w_0=1, w_j={\displaystyle \frac{j-1-\alpha }{j}}{w}_{j-1}.$$

(19)

Each coordinate slice uses its own positive endpoint distance. We take N = 64,128,…,8192 at four fractional orders, giving 32 paired direction tests. Figure 5 shows the absolute angular error relative to the exact circle result. At half order, the finest-grid angle is 30.018319 degrees. The exact value rounds to 30.008768 degrees. At orders 0.25, 0.5, 0.65 and 0.9, the finest angular errors are 0.004091, 0.009551, 0.011804 and 0.005324 degrees, respectively. These are finite-grid errors. We report them explicitly rather than describing the GL calculation as accurate to machine precision.

4.2. A Benchmark Using Published Modified Cam Clay Geometry

Sun and Zheng use a critical-state-directed Caputo construction on a modified Cam-clay function [48]. We reuse their polynomial geometry with our stated RL operator. This gives a test on a published yield surface. It also keeps three issues separate: the surface geometry, the derivative definition and validation against material observations.

We take the modified Cam-clay function from Equation (10) of Sun and Zheng [48]. Nondimensionalization and positive constant scaling give F_M = M²x(x − L) + y². We then apply the zero-terminal RL operator from Section 2, including all coordinate-independent terms. This procedure does not reproduce the published Caputo constitutive model. On F_M = 0, for x,y > 0, direct evaluation gives

$$\begin{array}{r}{g}_x={\displaystyle \frac{\alpha M^2{x}^{1-\alpha }}{\Gamma \left(2-\alpha \right)}}\left[{\displaystyle \frac{\left(3-\alpha \right)x}{2-\alpha }}-L\right]\\ g_y={\displaystyle \frac{\alpha \left(3-\alpha \right)y^{2-\alpha }}{\Gamma \left(3-\alpha \right)}}.\end{array}$$

(20)

A term that is constant along one coordinate slice still contributes to its RL derivative. Removing that term would change the test. We use five fractional orders, M = 0.8,1.2,1.6, L = 1 and x = 0.2,0.4,0.6,0.8, with $y=M\sqrt{x\left(1-x\right)}$. This gives 60 combinations. Independent monomial evaluation and endpoint-plus-integral quadrature agree with (20). Their maximum scaled errors are 6.75 × 10⁻¹⁶ and 4.78 × 10⁻¹⁶, respectively. S1 also compares the surface function with $F_M+0.5F_M^3$ at each state. This benchmark checks the operator calculation on a published geometry. It includes no laboratory calibration or material validation.

5. Mechanical Consequences

The mechanical diagnostic converts the fractional vector into a physical plastic-flow tensor. We call the complete constitutive system used for this comparison the host model. Such a conversion must be specified in transformed-stress formulations [15]. We use the SMP map, work-conjugate tensor conversion and complete two-size consistency equations from the companion AJOP-T formulation [35]. Appendix A and the executable S1 sources define the host without requiring access to that preprint. We compare direct RL differentiation of f and f + 0.5f³ at identical operator settings. Each comparison also fixes internal variables and material parameters. Nonlinear convergence and tangent consistency are checked separately [49]. The compression-map hardening comes from the published, rate-independent AJOP-T model [50]. The companion preprint supplies its fractional tensor formulation. The loading surface and size-evolution equations are

$$f=\Omega \ln {\displaystyle \frac{p}{P_{0}R}}+{\left({\displaystyle \frac{\tilde{q}}{\mathit{Mp}}}\right)}^{\Psi }=0, h_A={\displaystyle \frac{v}{{\lambda }_A-\kappa }},$$

(21)

$$\dot{\ln P_0}=\dot{\gamma }{h}_A{m}_v^d, \dot{\ln R}=\dot{\gamma }{h}_A\left(m_v^k-m_v^d\right), H=\Omega h_A{m}_v^k.$$

(22)

For the canonical-graph control, solve the positive-shear branch of (21) for $y=\tilde{q}/p_r$ at fixed P₀,R. With x = p/p_r, define

$$f_{\mathrm{graph}}\left(x,y;L\right)=y-\mathit{Mx}{\left[\Omega \ln \left(L/x\right)\right]}^{1/\Psi }, L=P_{0}R/p_r.$$

(23)

The real branch requires L > 0, M,Ω,Ψ > 0, 0 < x < L and y ≥ 0. Numerical evaluation further requires 0 < a_x < x < L and y > 0. The two RL integration intervals are [a_x,x] and [0,y], with L fixed. All sampled logarithms and fractional powers are then real. We differentiate the full graph, including terms that are constant in one coordinate. Before differentiation, both labels f and f + 0.5f³ are mapped to the same function (23). The parameters and hardening closure are also identical. Their equality is therefore imposed by the graph convention. It is a control calculation, rather than evidence of invariance under arbitrary scalar relabeling.

For direct RL flow, the component ratio is D^d = g_x/g_y, with g_y > 0. The overconsolidated (OC) closure is $D^k=D^d+\left(D_6^k-D_6^d\right)$. This is a diagnostic choice. In general, it does not recover the non-associated AJOP-T seed at α = 1. The inherited branch difference vanishes at R = 1. None of these mechanical choices is needed for the exact results in Section 2 and Section 3. Those sections test equivalence of fractional directions. The present host illustrates the consequences within one stated constitutive system.

The recorded campaign uses α = 0.65, a_x = 0.1, a_y = 0 and p_r = 91.6 kPa. It includes three initial Lode states and four initial R values. Appendix A gives the tensor conversion, parameters and implicit residual. S1 records every path. Table 1 reports the intermediate-Lode case with R₀ = 0.6 at imposed deviatoric strain norm 0.01, using 80 steps.

Figure 6 compares the responses at the same imposed strain. Positive plastic volume means contraction. The cubic member changes sign for the selected terminal. This behavior does not occur for every admissible terminal. We selected a_x = 0.1 after exploratory probes, and S1 retains all four choices. At the shorter strain norm 0.005, the base/cubic volume pairs are (0.184809, -0.027455), (0.206751, 0.206976), (0.208251, 0.210010) and (0.187628, 0.187757) percent for a_x = 0.1,0.25,0.5,1, respectively. Both members remain contractive for the last three terminal choices.

Admissibility is checked separately. At R₀ = 0.125, all three cubic cases reach the imposed boundary of nonnegative plastic power. Their consistency denominators remain positive. Independent continuous integrations locate the events at strain norms 0.263991%, 0.237932% and 0.265314% for initial Lode indicators 1, 0 and -1, respectively. The baseline cases complete the prescribed loading. Beyond an event, the cubic path fails the chosen acceptance test. It should not be interpreted as an admissible alternative prediction there. This result concerns the direct-RL diagnostic. It does not establish failure of the companion AJOP-T seed or a general thermodynamic dissipation theorem.

Common-state Lode dependence. We next compare the two physical plastic-flow directions at the same state. This removes the effect of different path histories. The angle uses the Frobenius tensor norm, which is represented exactly by the host’s Mandel vectors. We scan 401 values of ξ = cos(3θ) ∈ [−1,1]. The common settings are p = 150 kPa, R = 0.6, v = 2.4, α = 0.65, a_x = 0.1 and cubic coefficient 0.5. Family A fixes physical q = 60 kPa. It sets P₀ on the common surface at each ξ. Family B fixes $\tilde{q}=64.749764630$ kPa and P₀ = 290.672614416 kPa, then solves for physical q. Both representations have identical settings within each pair.

Table 2. Physical-flow separation at three common-state Lode anchors.

Lode indicator ξ	Angle in family A in degrees	Angle in family B in degrees
-1	32.392889	16.992010
0	17.810246	17.810246
1	6.907904	19.190149

Table 2. Physical-flow separation at three common-state Lode anchors.

Lode indicator ξ	Angle in family A in degrees	Angle in family B in degrees
-1	32.392889	16.992010
0	17.810246	17.810246
1	6.907904	19.190149

Figure 7 shows the two state families. The largest sampled angle occurs at extension (ξ = −1) in family A and compression (ξ = 1) in family B. These maxima apply to the stated samples. They are not global maxima over all stresses, terminals or material states. Both flow directions have positive stress power and positive consistency denominator at every sampled state. We also check the three family A anchors with 256-node and independent adaptive quadrature. The change in the location of the maximum shows why a Lode comparison must state which invariants are held fixed. This is a local direction test, not a new evolving loading-path model.

5.1. Mechanical Verification and Refinement

The campaign combines four flow rules, three initial Lode indicators and four initial mapping ratios. This gives 48 nominal variants. Forty-five variants complete each of the 20-, 40- and 80-step grids, giving 6300 accepted main-grid steps. Three cubic cases at initial R = 0.125 are rejected and remain in the rejection ledger. Six additional volume-changing paths use only the 80-step grid. Their 480 accepted steps are counted directly in paths.csv. The companion AJOP-T seed and canonical-graph rules are controls. In the graph control, both labels are converted to the same function before the same operator is applied. Equality of the resulting pair is therefore a check of that construction.

The 540 quadrature comparisons have maximum scaled differences of 1.843 × 10⁻¹⁰ against adaptive integration and 6.431 × 10⁻¹⁰ against 256-node quadrature. Across 80 rotated states, the relative work-conjugacy and rotation-covariance discrepancies are below 4.416 × 10⁻¹⁵ and 7.882 × 10⁻¹⁴, respectively. These tests check tensor conversion. They do not test invariance under a change in the defining function. Paired nonlinear solves for positive constant scaling and integer order differ by less than 5 × 10⁻¹³ in their largest endpoint component. The canonical and AJOP-T seed pairs coincide by construction for the selected relabeling. Their agreement does not independently validate a material model.

We also integrate the representative paths with DOP853, an explicit Runge–Kutta method described by Hairer, Nørsett and Wanner [51]. The independent end-stress separation is 50.682605 kPa. The sum of the observed 80-step errors relative to those continuous solutions is 0.225291 kPa. The separation is therefore approximately 225 times that observed error sum. This comparison is not a certified error bound. Observed orders on the coarse grids range from 0.701 to 0.997. The two slowly converging cubic compression cases reach orders 0.966 and 0.970 on the 80/160/320 grids. The nearly stationary extension case reaches 0.867, with a final-grid stress difference of 5.62 × 10⁻⁸ kPa. We therefore do not claim uniform first-order convergence. DOP853 also checks the power-boundary events at two tolerances. S1 retains the detailed errors, failed cases and terminal probes.

6. Discussion

The main result states that a yield boundary need not determine the fractional flow direction. Nonlocality alone suggests a dependence on values away from the boundary. Proposition 2 establishes more. The two functions remain globally strictly convex and describe the same bounded elastic region. Their derivatives agree everywhere on the boundary through any prescribed finite order. Their RL vectors still differ at the stated point for every fractional order between zero and one. Thus a finite amount of boundary data cannot determine the direction of this operator. The construction may require a larger polynomial power when the prescribed matching order increases.

A constitutive model may legitimately specify the complete scalar function, coordinates, scale and terminals. Changing that function then changes constitutive data, even when the boundary remains the same. If a model is intended to depend only on the boundary, equivalence must be proved for its chosen operator. Positive constant scaling passes our unit-direction test. The graph convention in (23) removes the particular input ambiguity by assigning both labels to one function. It introduces additional information and does not establish invariance under a change of coordinates.

The kernel, integration domain and definition of a vector operator are separate choices [26,30,31,32]. Our comparison fixes all of them. It shows that the full scalar defining function can be an additional constitutive input. It does not settle objectivity or equivalence for other fractional vector calculi.

Equal unit directions are also insufficient for complete constitutive equivalence. The internal evolution and plastic multiplier must transform consistently, as shown in (18). In the two-size diagnostic, the plastic-volume sign change depends on the lower terminal. Three cubic paths at small R leave the imposed nonnegative-power domain. Their rejection applies to that acceptance rule. It is not a universal thermodynamic conclusion.

Other studies examine plastic dissipation in particular fractional Cam-clay and thermodynamically structured models [17,18]. Our rejection test concerns the interaction between a scalar relabeling and the stated power condition. A complete stored-energy inequality would also need the energies and powers associated with all internal variables.

The contribution is distinct from the terminal and coordinate questions in the related Caputo study [34]. We attribute the tensor and hardening structure to the companion AJOP-T formulation [35]. The fractional-flow construction of Sumelka and Nowak [16] and the finite-support Riesz–Caputo operator of Feischl et al. [24] are different. Extending the strict determinant proof to either would require another derivation. Likewise, the modified Cam-clay test uses the surface from Sun and Zheng [48], while applying our RL operator rather than their Caputo rule.

Verification checks whether code implements the stated equations. Material validation instead requires independent physical observations [52]. Here, exact identities, independent discretizations and refinement checks isolate sensitivity to the defining function. Choosing a physically preferable representation would require further material evidence. This would include independently determined parameters and observations withheld from calibration, as illustrated by the types of sand and rock comparisons in [12,13].

The mathematical comparison uses a fixed positive-coordinate domain and stated terminals. The mechanical diagnostic adds monotonic active loading and excludes hydrostatic singularities. It does not include evolving fabric or chronological memory. These restrictions limit the mechanical illustration. They do not change the closed-boundary counterexample for the RL operator defined in Section 2.

7. Conclusions

For the stated componentwise left RL operator, a common yield boundary does not determine a unique normalized fractional direction. This remains true when the defining functions are globally strictly convex and their boundary derivatives agree through any prescribed finite order. The constructed pair has the same closed ellipse and bounded elastic region. Its determinant has a strict sign for every fractional order between zero and one. At half order, the circle example gives an exact angular separation of 30.008768 degrees.

Independent polynomial evaluation, quadrature and GL refinement support the analytical result. The published modified Cam-clay geometry provides an additional operator check. The recorded three-dimensional campaign illustrates a mechanical separation that persists under refinement. It also shows a terminal-dependent plastic-volume sign change and three power-boundary rejections. The complete defining function must therefore be stated whenever it forms part of the constitutive data. These results do not establish a calibrated material model or invalidate fractional flow in general.