PadeMuntzFractionalRiccati - crowlogic/arb4j GitHub Wiki

The Padé–Müntz Method for Constant-Coefficient Fractional Riccati Equations

Stephen Crowley April 2026

Abstract

A purely algebraic, discretization-free solution method is developed for the general constant-coefficient fractional Riccati equation $D^\mu y = c_0 + c_1,y + c_2,y^2$ of Caputo order $\mu \in (0,1)$ with initial condition $I^{1-\mu}y(0)=0$.

Because the Müntz lattice $\lbrace t^{k\mu}\rbrace_{k\ge1}$ is closed under the Caputo derivative, fractional integration, and pointwise multiplication, the Müntz spectral Tau substitution yields an explicit Puiseux expansion $y(t)=\sum_{k\ge1} a_k,t^{k\mu}$ whose coefficients satisfy a closed Gamma-ratio convolution recurrence. The induced power series in the variable $z=t^\mu$ has a finite, strictly positive radius of convergence $R$ controlled by movable singularities of the Riccati flow in the complex $z$-plane, so the recurrence by itself produces only a local representation and is not globally convergent in $t$.

Global-in-time convergence on the positive real axis is recovered by re-summing the Puiseux series via the diagonal $[M/M]$ Padé approximant in $z$, constructed from the Müntz coefficients through a single Hankel linear system $H_M\mathbf{q}=-\mathbf{b}$. Under the de Montessus de Ballore and Nuttall–Pommerenke theorems, the diagonal Padé sequence converges to the analytic continuation of $y$ along $z=t^\mu$ for all $t\ge0$ whenever the solution has no real-positive poles, and a computable a-posteriori error bound

$$|y(t)-y_M(t)|\le \frac{|\Delta_M(t^\mu)|^2}{|\Delta_{M-1}(t^\mu)|-|\Delta_M(t^\mu)|}$$

in terms of successive Padé differences $\Delta_k=R_k-R_{k-1}$ is derived from a geometric tail estimate.

A concluding section identifies the deeper structure: the substitution $z=t^\mu$ promotes the solution to a resurgent transseries in $z$, the Padé step is the computational instantiation of Borel–Padé resummation, and the global analytic continuation is governed by a scalar Riemann–Hilbert problem whose jump contours are the Stokes lines of the transseries. This perspective points to natural generalizations to higher-degree polynomial nonlinearities, fractional ODE systems via Hermite–Padé simultaneous approximation, and non-constant coefficients via isomonodromic Riemann–Hilbert problems of Painlevé type.

1. Fractional Calculus and the Müntz Basis

1.1 Riemann–Liouville integral and Caputo derivative

Definition (Riemann–Liouville integral). Let $r>0$ and $f\in L^1_{\mathrm{loc}}[0,T]$. Define

$$I^r f(t) := \frac{1}{\Gamma(r)}\int_0^t (t-s)^{r-1} f(s),ds,\qquad t\in[0,T].$$

Definition (Caputo derivative). Let $\mu\in(0,1)$ and $f$ satisfy $I^{1-\mu}f\in C^1[0,T]$. Define

$$D^\mu f(t) := \frac{d}{dt}!\left(I^{1-\mu}f(t)\right) - \frac{f(0)}{\Gamma(1-\mu)}t^{-\mu},\qquad t\in(0,T].$$

The fundamental relations are

$$I^\mu D^\mu f(t) = f(t) - f(0),\qquad D^\mu I^\mu f(t)=f(t).$$

Remark (Vanishing correction term for the Riccati solution). For the fractional Riccati equation below, the initial condition $I^{1-\mu}y(0)=0$ forces $y(0)=0$, since $I^{1-\mu}[t^{k\mu}]\big|_{t=0}=0$ for all $k\ge1$. Consequently the correction term $-f(0)t^{-\mu}/\Gamma(1-\mu)$ vanishes identically along the Riccati solution, and the Caputo derivative reduces to $D^\mu y(t)=\frac{d}{dt}I^{1-\mu}y(t)$ throughout.

Proposition (Action on powers). Let $s>-1$ and $r>0$. Then

$$I^r t^s = \frac{\Gamma(s+1)}{\Gamma(s+r+1)},t^{s+r},$$

$$D^\mu t^s = \frac{\Gamma(s+1)}{\Gamma(s+1-\mu)},t^{s-\mu}.$$

Proof. For $I^r t^s=\frac{1}{\Gamma(r)}\int_0^t(t-u)^{r-1}u^s,du$, the substitution $u=tx$ gives

$$I^r t^s = \frac{t^{s+r}}{\Gamma(r)}\int_0^1(1-x)^{r-1}x^s,dx = \frac{t^{s+r}}{\Gamma(r)}\cdot\frac{\Gamma(s+1)\Gamma(r)}{\Gamma(s+r+1)} = \frac{\Gamma(s+1)}{\Gamma(s+r+1)},t^{s+r}.$$

For $D^\mu t^s$, $I^{1-\mu}t^s=\frac{\Gamma(s+1)}{\Gamma(s+2-\mu)}t^{s+1-\mu}$. Differentiating and using $\Gamma(s+2-\mu)=(s+1-\mu)\Gamma(s+1-\mu)$ gives the result. Since $t^s$ is continuous at $0$ for $s>-1$, the Caputo correction is absent. ∎

1.2 The Müntz basis

Fix $\mu\in(0,1)$ and $T>0$.

Definition (Müntz system). For $k\in\mathbb{N}$, define $e_k(t):=t^{k\mu}$ on $[0,T]$. The span of $\lbrace 1,e_k:k\ge1\rbrace$ is denoted

$$\mathcal{M} := \mathrm{span}{1,t^{\mu},t^{2\mu},\ldots}.$$

Remark (Completeness). By the Müntz–Szász theorem, $\lbrace t^{k\mu}\rbrace_{k\ge1}$ is complete in $L^2[0,1]$ iff $\sum_{k\ge1}(k\mu)^{-1}=\infty$. Since $\sum_{k\ge1}k^{-1}=\infty$, completeness holds for every $\mu\in(0,1)$. This underwrites the Müntz–Tau spectral method.

Proposition (Closure of $\mathcal{M}$). For $k\ge1$ and $r>0$,

$$D^\mu t^{k\mu} = \frac{\Gamma(k\mu+1)}{\Gamma((k-1)\mu+1)},t^{(k-1)\mu},$$

$$I^r t^{k\mu} = \frac{\Gamma(k\mu+1)}{\Gamma(k\mu+r+1)},t^{k\mu+r},$$

$$t^{j\mu},t^{k\mu} = t^{(j+k)\mu}.$$

Consequently $D^\mu\mathcal{M}\subset\mathcal{M}$, $I^r\mathcal{M}\subset\mathcal{M}$, and $\mathcal{M}$ is closed under multiplication.

These three closure properties together — under $D^\mu$, under $I^r$, and under multiplication — are exactly what makes the Müntz lattice the canonical basis for any constant-coefficient fractional ODE with polynomial nonlinearity. The Müntz–Tau substitution preserves the entire algebraic structure of the equation.

2. The Constant-Coefficient Fractional Riccati Equation

2.1 Equation and Volterra representation

Let $\mu\in(0,1)$ and $c_0,c_1,c_2\in\mathbb{C}$. Consider

$$D^\mu y(t)=c_0+c_1,y(t)+c_2,y(t)^2,\qquad I^{1-\mu}y(0)=0.$$

Applying $I^\mu$ and using $I^\mu D^\mu f = f - f(0)$ (and $y(0)=0$) gives the equivalent Volterra integral equation

$$y(t) = \frac{1}{\Gamma(\mu)}\int_0^t (t-s)^{\mu-1}!\left(c_0+c_1,y(s)+c_2,y(s)^2\right)ds.$$

Remark (Hölder regularity at the origin). Evaluating the Volterra equation at $t=0$ gives $y(0)=0$. Near the origin, $y(t)=\frac{c_0}{\Gamma(\mu+1)},t^\mu + \mathcal{O}(t^{2\mu})$, so $y$ is Hölder continuous of order $\mu$ but not $C^1$ at $0$. The Hölder exponent $\mu$ is the natural regularity scale of the problem and the reason ordinary Taylor series in $t$ fail.

2.2 The Müntz–Tau Puiseux series solution

Theorem (Series solution in $\lbrace t^{k\mu}\rbrace$). There exists $R>0$ and coefficients $(a_k)_{k\ge1}\subset\mathbb{C}$ such that

$$y(t)=\sum_{k=1}^{\infty}a_k,t^{k\mu},\qquad |t|<R^{1/\mu},$$

with

$$a_1 = \frac{c_0}{\Gamma(\mu+1)},$$

$$a_{k+1} = \frac{\Gamma(k\mu+1)}{\Gamma((k+1)\mu+1)} \left(c_1,a_k+c_2!!\sum_{\substack{j+\ell=k\1\le j,\ell\le k-1}}!!a_j,a_\ell\right),\qquad k\ge1.$$

Remark (Quadratic convolution at $k=1$). For $k=1$, the convolution sum is empty (no pairs $(j,\ell)$ with $j,\ell\ge1$ and $j+\ell=1$), so $a_2 = \frac{\Gamma(\mu+1)}{\Gamma(2\mu+1)},c_1,a_1$. The quadratic term first contributes at $k=2$, giving $a_3 = \frac{\Gamma(2\mu+1)}{\Gamma(3\mu+1)}(c_1,a_2 + c_2,a_1^2)$.

Proof. Assume $y(t)=\sum_{k\ge1}a_k,t^{k\mu}$ with radius of convergence $R^{1/\mu}>0$. Then by closure of the Müntz lattice under multiplication,

$$y(t)^2=\sum_{m=2}^{\infty}!\Big(\sum_{\substack{j+\ell=m\1\le j,\ell\le m-1}}!a_j,a_\ell\Big)t^{m\mu},$$

which we abbreviate $b_m:=\sum_{j+\ell=m}a_j,a_\ell$. Substituting into the Volterra equation and splitting the integral, the first integral equals $t^\mu/\Gamma(\mu+1)$. For the others, the action-on-powers identity with $s=k\mu$ gives

$$\int_0^t (t-s)^{\mu-1}s^{k\mu},ds = \Gamma(\mu),\frac{\Gamma(k\mu+1)}{\Gamma((k+1)\mu+1)},t^{(k+1)\mu}.$$

Reindexing with $n=k+1$ and $n=m+1$ and matching the coefficient of $t^{n\mu}$ for each $n\ge1$ yields the recurrence. Existence of $R>0$ follows from a standard fixed-point or majorant-series argument. ∎

Theorem (Finite Puiseux radius). The Puiseux series $\sum_{k\ge1}a_k,z^k$ in $z=t^\mu$ has a finite, strictly positive radius of convergence $R$. Specifically, with $M_0=|c_0|$,

$$|a_k|;\le;\frac{M_0,C^{k-1}}{\Gamma(k\mu+1)},$$

where $C>0$ depends on $|c_1|$ and $|c_2|$, so that on any finite truncation

$$\sum_{k=1}^{N}|a_k|,|z|^k ;\le; \frac{M_0}{C}\big(E_\mu(C|z|)-1\big),$$

where $E_\mu$ is the Mittag-Leffler function. The radius is finite (not infinite) because the function $g(z)=y(z^{1/\mu})$ generically develops movable singularities — the Riccati poles in the complex $z$-plane.

Proof. Let $B_k := M_0,C^{k-1}/\Gamma(k\mu+1)$. For the base case, $B_1=M_0/\Gamma(\mu+1)\ge|a_1|$. The convolution majorant satisfies

$$\sum_{\substack{j+\ell=k\j,\ell\ge1}}!\frac{1}{\Gamma(j\mu+1),\Gamma(\ell\mu+1)} ;\le;\frac{2^k}{\Gamma(k\mu+1)},$$

via

$$\frac{\Gamma(k\mu+1)}{\Gamma(j\mu+1),\Gamma(\ell\mu+1)}\le\binom{k}{j}\quad\text{for }\mu\le1\text{ and }j+\ell=k,$$

a Beta-function consequence of $B(j\mu+1,\ell\mu+1)\ge B(j+1,\ell+1)$. Substituting the inductive hypothesis into the recurrence gives

$$|a_{k+1}|;\le;\frac{M_0,C^{k-1}}{\Gamma((k+1)\mu+1)}\big(|c_1|+|c_2|,M_0,2^k/C\big).$$

Choose $C$ (depending on the truncation order $N$) with $C\ge\max(|c_1|,2|c_2|M_0,2^N)$; the induction closes. The bound on the truncation follows from $\sum_k z^k/\Gamma(k\mu+1)=E_\mu(z)$. The finiteness of $R$ is geometric: the complexified Riccati flow develops poles at $z$ with $|z|=R$, so the radius is bounded above by the distance to the nearest such pole. ∎

Remark (Why this is only a local representation). This is the central obstruction: even though every $a_k$ is computed in closed form by the recurrence, the series does not converge for all $t\ge0$. The complex $z$-plane carries movable poles of the Riccati flow, and the Puiseux radius $R$ is the distance to the nearest such pole. On the positive real $t$-axis the solution $y(t)$ may exist for all $t$, yet the Müntz series ceases to converge once $t^\mu$ exceeds $R$. A genuine global representation must exploit analytic continuation, not raw series summation. This is the role played by Padé approximation in the next section.

3. Padé Resummation in $z=t^\mu$

3.1 Power series in $z$ and the Hankel system

Let $z:=t^\mu$ and define $g(z):=y(z^{1/\mu})$. Then

$$g(z)=\sum_{k=1}^{\infty}a_k,z^k,\qquad |z|<R.$$

Fix $M\in\mathbb{N}$ and seek polynomials

$$P_M(z):=\sum_{k=0}^{M}p_k,z^k,\qquad Q_M(z):=1+\sum_{k=1}^{M}q_k,z^k,$$

satisfying the Padé matching condition

$$Q_M(z),g(z)-P_M(z)=\mathcal{O}(z^{2M+1})\quad(z\to 0).$$

Lemma (Hankel form of the denominator system). Substituting the series into the matching condition and matching coefficients of $z^n$ for $0\le n\le 2M$ gives $p_0=0$,

$$p_n = a_n+\sum_{j=1}^{\min(n,M)}q_j,a_{n-j},\qquad 1\le n\le M,$$

$$0 = a_n+\sum_{j=1}^{M}q_j,a_{n-j},\qquad M+1\le n\le 2M.$$

The denominator system is the Hankel linear system

$$H_M,\mathbf{q}=-\mathbf{b},\qquad H_M=\begin{pmatrix} a_M & a_{M-1} & \cdots & a_1\ a_{M+1} & a_M & \cdots & a_2\ \vdots & \vdots & \ddots & \vdots\ a_{2M-1} & a_{2M-2} & \cdots & a_M \end{pmatrix},\quad \mathbf{q}=\begin{pmatrix}q_1\ \vdots\ q_M\end{pmatrix},\quad \mathbf{b}=\begin{pmatrix}a_{M+1}\ \vdots\ a_{2M}\end{pmatrix}.$$

Proof. Direct expansion of $Q_M,g$ and identification of coefficients. ∎

Remark (Invertibility of $H_M$). $H_M$ is invertible iff $g(z)$ is not a rational function of degree less than $M$. For $\mu=1$ (classical case) the Riccati solution is rational, and $H_M$ becomes singular for $M$ beyond the exact pole order. For $\mu\in(0,1)$, $g(z)$ admits no finite-order rational representation; $\det H_M\ne 0$ generically for all $M\ge1$. If a particular order $M$ yields $\det H_M=0$, use $M\pm 1$.

Lemma (Solution). If $\det H_M\neq0$, then the system has the unique solution $\mathbf{q}=-H_M^{-1}\mathbf{b}$, and the numerator coefficients are given by the formula above.

The diagonal Padé approximant is then

$$R_M(z) := \frac{P_M(z)}{Q_M(z)} = \frac{\sum_{k=0}^{M}p_k,z^k}{1+\sum_{k=1}^{M}q_k,z^k}.$$

Definition (Padé domain). $\mathcal{D}_M:=\lbrace z\in\mathbb{C}:\det H_M\neq0,\ Q_M(z)\neq0\rbrace$.

3.2 Remainder representation and successive-difference bound

Define the remainder

$$\varepsilon_M(z):=g(z)-R_M(z).$$

By construction $\varepsilon_M(z)=\mathcal{O}(z^{2M+1})$, and standard Padé theory gives the structural representation

$$\varepsilon_M(z)=\frac{\Pi_{M+1}(z)}{Q_M(z)^2},z^{2M+1}$$

for some polynomial $\Pi_{M+1}$ of degree at most $M+1$.

For the a-posteriori bound, define

$$\Delta_k(z):=R_k(z)-R_{k-1}(z),\qquad k\ge1.$$

Assumption (Geometric tail). At a fixed $z\in\mathcal{D}_M\cap\mathcal{D}_{M-1}$:

$g$ has no poles in $\lbrace |\zeta|\le|z|\rbrace$;
$Q_M(z)\ne0$ and $Q_{M-1}(z)\ne0$;
the sequence $|\Delta_k(z)|$ is geometrically decreasing for $k\ge M$, i.e. $|\Delta_M(z)|<|\Delta_{M-1}(z)|$.

Theorem (Computable error bound). Under the assumption,

$$|g(z)-R_M(z)|;\le;\frac{|\Delta_M(z)|^2}{|\Delta_{M-1}(z)|-|\Delta_M(z)|}.$$

Proof. The telescoping identity $g(z)-R_M(z)=\sum_{k=M+1}^{\infty}\Delta_k(z)$ holds wherever the Padé sequence converges. Setting $\rho:=|\Delta_M(z)|/|\Delta_{M-1}(z)|<1$, the geometric contraction $|\Delta_{j+1}|\le\rho,|\Delta_j|$ for $j\ge M$ yields $|\Delta_{M+k}|\le\rho^k|\Delta_M|$, hence

$$|g(z)-R_M(z)|\le\sum_{k=1}^{\infty}\rho^k|\Delta_M(z)|=\frac{\rho|\Delta_M(z)|}{1-\rho}=\frac{|\Delta_M(z)|^2}{|\Delta_{M-1}(z)|-|\Delta_M(z)|}.\quad ∎$$

3.3 Global convergence on the positive real axis

Theorem (de Montessus de Ballore). Let $g$ be meromorphic in $\lbrace |z|<\rho\rbrace$ with exactly $n$ poles (counted with multiplicity) in that disc and no poles on the boundary. The diagonal Padé approximants $R_M$ of order $(M,M)$ converge to $g$ uniformly on compact subsets of $\lbrace |z|<\rho\rbrace\setminus\lbrace\text{poles}\rbrace$ as $M\to\infty$.

Theorem (Nuttall–Pommerenke). Let $g$ be meromorphic in $\mathbb{C}$ with at most countably many poles. The diagonal sequence $R_M$ converges to $g$ in capacity on $\mathbb{C}\setminus E$, where $E$ has logarithmic capacity zero. In particular, convergence is quasi-everywhere, though not necessarily pointwise at every individual point.

Remark (Which theorem applies). The de Montessus theorem gives pointwise uniform convergence on compacta but requires knowing the exact pole count. The Nuttall–Pommerenke theorem requires no such count but yields only quasi-everywhere convergence. When $g(z)=y(z^{1/\mu})$ is meromorphic in $\mathbb{C}$ with no poles on $(0,\infty)$ — the canonical case for globally bounded Riccati solutions on the real time axis — both theorems apply, and $R_M(t^\mu)\to y(t)$ for every $t>0$.

Theorem (Global Padé–Müntz representation). Let $y(t)$ solve the fractional Riccati equation and assume $y$ is bounded on $[0,\infty)$ (so that $g(z)$ has no poles on the positive real $z$-axis). Then

$$y(t)=\lim_{M\to\infty}R_M(t^\mu)=\lim_{M\to\infty}\frac{P_M(t^\mu)}{Q_M(t^\mu)}\quad\text{for every }t\ge0.$$

Proof. By the finite-radius theorem the singularities of $g$ at $|z|=R$ are movable poles in the complex $z$-plane. Since $y$ is bounded on $[0,\infty)$, $g$ has no poles on the positive real axis. $g$ is meromorphic in $\mathbb{C}$ with poles accumulating only at $\infty$, so Nuttall–Pommerenke applies to give convergence on $(0,\infty)$, and de Montessus de Ballore gives uniform convergence on compacta away from poles. ∎

3.4 Algorithmic summary

Given $(c_0,c_1,c_2,\mu)$ and Padé order $M$:

Step 1. Compute the Puiseux coefficients $a_1,\ldots,a_{2M}$ by the recurrence. Cost: $\mathcal{O}(M^2)$.
Step 2. Form the Hankel matrix $H_M$ and right-hand side $\mathbf{b}$.
Step 3. Solve $H_M\mathbf{q}=-\mathbf{b}$ for the denominator coefficients. Cost: $\mathcal{O}(M^3)$.
Step 4. Recover the numerator coefficients via the explicit formula.
Step 5. Evaluate $y_M(t)=P_M(t^\mu)/Q_M(t^\mu)$ at any $t\ge0$. Cost: $\mathcal{O}(M)$ per query.

Steps 1–4 are independent of $t$; the entire family $\lbrace y_M(t):t\ge0\rbrace$ is encoded in a single rational function in $z=t^\mu$. There is no grid, no time stepping, no Newton iteration, no fixed-point loop.

4. Resurgence, Borel–Padé Resummation, and Riemann–Hilbert

The Padé–Müntz construction of the previous sections solves the constant-coefficient fractional Riccati equation in closed algebraic form by exploiting two independent structural facts: the Müntz lattice closes the differential algebra of $D^\mu,,I^r,,\cdot$; and the diagonal Padé in $z=t^\mu$ re-sums the resulting Puiseux series past its finite radius of analyticity. These two facts are not coincidences. Both are special cases of a deeper and far more general principle drawn from the theory of resurgent functions and the Riemann–Hilbert correspondence.

4.1 The substitution $z=t^\mu$ is an algebraization

A constant-coefficient fractional ODE in $t$ is not algebraic in $t$ — the operator $D^\mu$ is non-local, and the solution is generically only Hölder continuous of exponent $\mu$ at the origin. After the substitution $z=t^\mu$, however, the Volterra equation becomes equivalent to an algebraic ODE in $z$, in the precise sense that its solutions are Puiseux-expandable at the origin and the Cauchy–Kovalevskaya theorem applies. The substitution $z=t^\mu$ thus performs an algebraization: it lifts the problem from the Hölder category into the analytic category, where the full machinery of complex analysis becomes available.

4.2 Solutions are resurgent transseries

Algebraic ODEs with polynomial nonlinearities have solutions that are resurgent in the sense of Écalle: their formal expansions at any singular point can be promoted to transseries

$$\widetilde y(z);=;\sum_{n\ge0}e^{-n,A/z},z^{n\beta_n}\sum_{k\ge0}a_{n,k},z^k,$$

where the leading sector $n=0$ is the Puiseux series and the higher sectors $n\ge1$ encode non-perturbative contributions controlled by the action $A$ (the Riccati pole structure). The alien derivatives $\Delta_A\widetilde y$ of Écalle measure the obstruction to analytic continuation across each singularity, and the Stokes constants $S_A$ record the discontinuities of the lateral Borel sums when crossing Stokes lines. For the constant-coefficient fractional Riccati, the leading-sector coefficients are exactly the Müntz $a_k$; the higher sectors capture the multivaluedness of $y$ around the movable poles.

4.3 Padé is computational Borel–Padé resummation

The diagonal Padé approximant in $z$ is the canonical computational instantiation of Borel–Padé resummation. Recall that the Borel transform sends a divergent or finite-radius series $\sum a_k z^k$ to its Borel transform $\widehat g(\zeta)=\sum a_k\zeta^{k-1}/\Gamma(k)$, then resums by Laplace-transforming $\widehat g$ along a ray. The Padé approximant $R_M$ computes essentially the same object: the rational function $P_M/Q_M$ is the unique meromorphic continuation that interpolates the formal series and respects the underlying singularity structure. Concretely, the de Montessus and Nuttall–Pommerenke theorems say that $R_M$ recovers exactly the lateral Borel–Laplace sum on the real $z$-axis whenever no Stokes ray crosses $(0,\infty)$; this is precisely the boundedness hypothesis.

4.4 The Riemann–Hilbert structure

The deeper geometric content is this. The poles of $R_M$ for finite $M$ are not arbitrary — as $M\to\infty$, they cluster along Stokes lines in the $z$-plane, which are the branch cuts of the analytic continuation of $y$. Nuttall–Pommerenke is a weak version of this fact: convergence in capacity means the spurious poles distribute themselves along the natural boundary of analyticity. The strong version, due to Stahl, Gonchar, and others, identifies the equilibrium measure on those Stokes lines as the limiting pole distribution of the diagonal Padé sequence. This equilibrium problem is precisely a scalar Riemann–Hilbert problem for the resolvent of a logarithmic potential, and the diagonal Padé asymptotics are obtained by the Deift–Zhou steepest descent / Fokas–Its–Kitaev machinery applied to the associated jump matrix.

The unifying statement of this paper, in maximum generality, is therefore:

For a constant-coefficient fractional ODE $D^\mu y=F(y)$ with polynomial $F$, the substitution $z=t^\mu$ promotes the solution to a resurgent transseries in $z$. The diagonal Padé in $z$ is a Borel–Padé resummation of that transseries, and its global analytic continuation along the positive real axis is governed by a Riemann–Hilbert problem whose jump contours are the Stokes lines of the transseries.

4.5 Generalizations

The framework of this paper extends in three natural directions, each obtained by relaxing one structural hypothesis.

Higher-degree polynomial nonlinearities (Bernoulli, Abel, and beyond). For $D^\mu y=F(y)$ with $F$ any polynomial of degree $d$, the Müntz lattice remains closed under the nonlinearity by the closure-under-multiplication property, so the Tau substitution yields a closed Gamma-ratio recurrence of the same form but with a $d$-fold convolution sum. The Padé step is unchanged, the Hankel system is unchanged, and the de Montessus / Nuttall–Pommerenke theorems apply identically. Cubic Riccati (Abel of the second kind), quartic, and arbitrary polynomial $F$ are all immediate.

Fractional ODE systems via Hermite–Padé. For a system $D^\mu \mathbf{y} = \mathbf{F}(\mathbf{y})$ with polynomial $\mathbf{F}$, each component admits a Müntz Puiseux expansion, and the scalar Hankel system is replaced by a matrix Padé (or simultaneous Hermite–Padé) system involving multiple denominators. The convergence theory generalizes: the scalar Riemann–Hilbert problem becomes matrix-valued, with jump contours determined by the eigenvalue structure of the linearization at the movable singularities. Multi-factor rough Heston, multi-asset rough volatility, and multi-component fractional reaction–diffusion all fall in this class.

Non-constant coefficients: isomonodromy and Painlevé. When $D^\mu y = c_0(t) + c_1(t) y + c_2(t) y^2$ with non-constant $c_i$ (the natural setting for the rough Heston model when leverage and feedback are time-varying), the Müntz lattice is generally no longer preserved (the coefficients can introduce new exponents), the explicit Gamma-ratio recurrence breaks down, and the Puiseux structure may need to be enlarged to include logarithmic and Fuchsian terms. The associated Riemann–Hilbert problem is then isomonodromic rather than equilibrium-type, and the nonlinear special functions that arise — notably the Painlevé transcendents and their $q$- and fractional analogs — replace the elementary rational Padé approximants. This is also where the rough volatility literature stalls, and where new mathematics is genuinely required.

4.6 Why this combination was missed

The community gap that this construction fills is sociological as much as mathematical. Fractional-calculus practitioners optimized for numerical Volterra integrators (Adams–Bashforth–Moulton, multi-factor Markovian approximants, hybrid Euler–Richardson schemes) because the convolution structure makes them natural; they did not push into resummation theory. The Padé and resurgence communities developed deep results for integer-order ODEs and divergent-series problems but did not extend systematically into the fractional setting. The rough-volatility community needed a fast empirical answer and obtained the Gatheral–Radoičić two-point Padé fix, which works numerically but is not derived from the canonical algebraic structure of the equation. No one sat at the intersection of all three and asked the structural question — what is the canonical algebraic substitution that linearizes the differential algebra of a constant-coefficient fractional ODE? The answer, as this paper makes clear, is $z=t^\mu$, and once that substitution is made, the rest of the construction is completely determined: Müntz lattice closure forces the Tau recurrence, Puiseux finiteness forces the need for resummation, and meromorphy forces diagonal Padé.

References

O. El Euch and M. Rosenbaum, The characteristic function of rough Heston models, Mathematical Finance, 29(1):3–38, 2019.
G. A. Baker and P. Graves-Morris, Padé Approximants, 2nd ed., Cambridge University Press, 1996.
S. Esmaeili and M. Shamsi, The Müntz–Legendre Tau method for fractional differential equations, Applied Mathematical Modelling, 40(2):671–684, 2015.
R. de Montessus de Ballore, Sur les fractions continues algébriques, Bull. Soc. Math. France, 30:28–36, 1902.
J. Nuttall, The convergence of Padé approximants of meromorphic functions, J. Math. Anal. Appl., 31(1):147–153, 1970.
Ch. Pommerenke, Padé approximants and convergence in capacity, J. Math. Anal. Appl., 41:775–780, 1973.
H. Stahl, The convergence of Padé approximants to functions with branch points, J. Approx. Theory, 91:139–204, 1997.
J. Écalle, Les fonctions résurgentes, Vols. I–III, Publ. Math. d'Orsay, 1981–1985.
O. Costin, Asymptotics and Borel Summability, Chapman & Hall/CRC, 2008.
P. Deift, Orthogonal Polynomials and Random Matrices: A Riemann–Hilbert Approach, Courant Lecture Notes 3, AMS, 1999.
G. Callegaro, M. Grasselli, and G. Pagès, Fast hybrid schemes for fractional Riccati equations (rough is not so tough), Mathematics of Operations Research, 46(1):221–254, 2021.
J. Gatheral and R. Radoičić, Rational approximation of the rough Heston solution, International Journal of Theoretical and Applied Finance, 22(3):1950010, 2019.