\documentclass[12pt,a4paper]{article}

\usepackage{amsmath,amssymb,amsthm}
\usepackage{geometry}
\usepackage{hyperref}
\usepackage{booktabs}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage{cite}
\usepackage{caption}
\graphicspath{{./figures/}{./}}

\geometry{margin=2.5cm}

\newtheorem{theorem}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}{Definition}
\newtheorem{remark}{Remark}

\title{\textbf{Differential Nonlinear Robustness of Critical States\\
in Fibonacci and Tribonacci Substitution Chains}\\[0.5em]
\large A Numerical Study of Discrete Nonlinear Schr\"{o}dinger Dynamics\\
\large \emph{Version 3 --- finite-size scaling and long-time evolution}}

\author{Pablo Nogueira Grossi\\
G6 LLC, Newark, New Jersey, USA\\
ORCID: 0009-0000-6496-2186\\
\texttt{pablogrossi@hotmail.com}\\[0.3em]
\small Zenodo (V4): \href{https://doi.org/10.5281/zenodo.20075822}{10.5281/zenodo.20075822} \\
\small Concept DOI (auto-resolves to latest): \href{https://doi.org/10.5281/zenodo.20026942}{10.5281/zenodo.20026942}}

\date{May 2026 (v3)}

\begin{document}

\maketitle

\begin{abstract}
We study the discrete nonlinear Schr\"{o}dinger (DNLS) equation on two
quasiperiodic tight-binding chains: the Fibonacci chain ($n=2$) and the
Rauzy--tribonacci chain ($n=3$), generated by their respective substitution
rules. Starting from mid-gap eigenstates of the linear Hamiltonian, we
integrate the DNLS time evolution over a range of nonlinearity strengths
$\lambda \in [0,10]$ and times $T \in [50, 10^6]$, and track the inverse
participation ratio (IPR) as a measure of localization. The tribonacci chain
exhibits a mid-gap IPR approximately four times larger than the Fibonacci
chain in the linear limit, consistent with the stronger multifractality of
Rauzy fractal eigenstates. Under nonlinear perturbation at $T=50$, the
Fibonacci mid-gap state delocalizes rapidly (IPR drops $\sim 57\%$ at
$\lambda=1.5$), while the tribonacci state remains nearly pinned (IPR drop
$<5\%$ at the same coupling). We term this \emph{differential nonlinear
robustness}.

In this revision, we extend the analysis along two axes. \textbf{Finite-size
scaling} at $T=10^4$ across $N \in \{500, 1000, 2000\}$ shows that the
differential ratio $\mathrm{IPR}_{\rm trib}/\mathrm{IPR}_{\rm fib}$
generally grows with $N$ at fixed $\lambda$, supporting persistence of the
differential in the thermodynamic limit. \textbf{Long-time evolution} at
$N=1000$, $\lambda=1.5$ to $T=10^6$ shows that both chains saturate to a
finite-size limited IPR ($\approx 0.002$) by $T \sim 3\times10^5$, after
which the ratio oscillates around $1.04 \pm 0.04$. The $T=10^4$ peak ratio of
$\approx 4.2$ at $N=1000$ is therefore a transient feature; at $T=10^5$ it
collapses to $\approx 1.20$, and at $T=10^6$ the differential dissolves into
the finite-size saturation band. The pre-saturation spreading rates fitted
in the window $t > 10^4$ at $N=2000$, $T=10^5$ satisfy $\alpha_{\rm trib} >
\alpha_{\rm fib}$ uniformly across $\lambda \in [0.5, 2.0]$, indicating that
tribonacci spreads faster towards saturation but starts from a higher initial
IPR. The mid-gap multifractal dimension fit at natural Fibonacci lengths
gives $D_2^{\rm fib} = 0.65 \pm 0.11$ ($R^2 = 0.90$). The corresponding fit
for tribonacci at natural Rauzy lengths exhibits a plateau across iterations
$n=9,10,11$ followed by a discontinuous drop, indicating that the closest-to-
$E=0$ mid-gap state selector does not consistently track the same physical
eigenstate across substitution iterations; we therefore defer numerical
extraction of $D_2^{\rm trib}$ pending state-tracking diagnostics.

The characteristic decay constant of the tribonacci amplitude envelope,
$\eta \approx 1.839287$, is the Perron--Frobenius eigenvalue of the
tribonacci companion matrix, whose existence and key properties ($\eta > 1$,
strict antitonicity of $\eta^{-k}$) are formally verified in
Lean~4/Mathlib4 in the companion AXLE repository \cite{AXLE2026}. To our
knowledge, this is the first numerical study of DNLS dynamics on a
tribonacci substitution chain, and the present revision contributes the
first finite-size scaling and long-time saturation analysis identifying a
qualitative difference in nonlinear response between $n=2$ and $n=3$
substitution-generated quasicrystal models.
\end{abstract}

\textbf{Keywords:} tribonacci chain; Fibonacci quasicrystal; discrete
nonlinear Schr\"{o}dinger equation; multifractality; inverse participation
ratio; finite-size scaling; spreading exponent; Rauzy substitution

\section{Introduction}

Quasiperiodic tight-binding chains generated by substitution rules have been
a canonical setting for studying the interplay between long-range order,
multifractality, and localization since the 1980s. The Fibonacci chain---the
$n=2$ case of the $n$-bonacci family---occupies a central role: its spectrum
is a Cantor set of measure zero, its eigenstates are critical (neither
extended nor localized), and their multifractal dimensions are precisely
computable \cite{Kohmoto1983,Ostlund1983,Luck1989,Mace2017}. Topological
phases, charge pumping, and one-dimensional quasicrystalline superconductivity
have been analyzed in the Fibonacci setting \cite{Kraus2012,Lang2012}.

The natural generalization to tribonacci and higher $n$-bonacci substitution
chains is more recent. Krebbekx and collaborators \cite{Krebbekx2023}
introduced two tribonacci tight-binding models (hopping and on-site
modulation) via Rauzy substitutions, established that their spectra inherit
Rauzy-fractal internal space structure, and computed multifractal dimensions
showing stronger spectral clustering than in the Fibonacci case.
Higher-$n$ generalizations appear in the context of Rauzy tilings and
cut-and-project schemes, but exclusively in the linear (Schr\"{o}dinger
operator) setting \cite{Bellissard1989,Damanik2008}.

The nonlinear regime is comparatively unexplored. For the Fibonacci chain,
discrete nonlinear Schr\"{o}dinger (DNLS) dynamics have been studied in
photonic quasicrystals \cite{Lahini2009}, edge-soliton pumping has been
demonstrated \cite{Jurgensen2021}, and phason dynamics under nonlinear
perturbation have been analyzed \cite{Verbin2015}. For $n \geq 3$
substitution chains, however, no published work has placed an $n$-bonacci
recurrence structure inside a DNLS or Gross--Pitaevskii nonlinearity, nor
analyzed the resulting localization dynamics. This gap motivates the present
study.

The original version of this paper (v1, May 2026) reported the $T=50$
nonlinear baseline, establishing differential nonlinear robustness at
$N=500$ across $\lambda \in [0,10]$. The present revision (v3) extends the
analysis to (i) finite-size scaling across $N \in \{500, 1000, 2000\}$ at
$T=10^4$, (ii) long-time evolution at $N=1000$, $\lambda=1.5$ to
$T=10^6$, (iii) pre-saturation spreading rates fit in the window
$t > 10^4$ at $T=10^5$, and (iv) multifractal dimension extraction at
natural substitution lengths. All numerical computations are reproducible
from the open-source repository at
\href{https://github.com/grossi-ops/Atratores}{github.com/grossi-ops/Atratores}.

A secondary motivation comes from the dm$^3$ contact-geometric framework
\cite{Grossi2026_Vol1,Grossi2026_Wavenumber6}, where the tribonacci
recurrence $w(k+3) = w(k+2) + w(k+1) + w(k)$ arises from the three-axis
decomposition of a contact 3-manifold, and the geometric weight $\eta^{-k}$
with $\eta \approx 1.839287$ (the tribonacci constant) is the natural scaling
factor. The formal verification that $\eta > 1$ and that the weight sequence
$\{\eta^{-k}\}$ is strictly antitone---establishing that the amplitude
envelope is well-defined and decaying---is carried out in Lean~4/Mathlib4 in
the AXLE repository \cite{AXLE2026}, specifically in
\texttt{TribonacciMeasure.lean}. We cite this formal result in Section~3 as
the machine-verified foundation for the decay envelope used in our ansatz.

\section{The $n$-Bonacci Recurrence and Substitution Chains}

\subsection{The $n$-Bonacci Recurrence}

\begin{definition}[$n$-Bonacci Sequence]
For $n \geq 2$, the $n$-bonacci sequence $\{A_k^{(n)}\}_{k \geq 0}$ is
defined by
\[
  A_k^{(n)} = \sum_{i=1}^{n} A_{k-i}^{(n)}, \quad k \geq n,
\]
with initial conditions $A_0^{(n)} = \cdots = A_{n-2}^{(n)} = 0$,
$A_{n-1}^{(n)} = 1$.
\end{definition}

The characteristic polynomial is
\[
  p_n(x) = x^n - x^{n-1} - x^{n-2} - \cdots - x - 1 = 0.
\]
By the Perron--Frobenius theorem applied to the companion matrix $C_n$,
there exists a unique dominant real root $\rho_n > 1$, the
\emph{$n$-bonacci constant}. The sequence $\rho_n$ is strictly increasing
and satisfies $\rho_n \to 2$ monotonically. Numerical values:

\begin{center}
\begin{tabular}{ccl}
\toprule
$n$ & $\rho_n$ (6 d.p.) & Name \\
\midrule
2 & 1.618034 & Golden ratio $\varphi$ \\
3 & 1.839287 & Tribonacci constant $\eta$ \\
4 & 1.927562 & Tetranacci constant \\
5 & 1.965948 & Pentanacci constant \\
\bottomrule
\end{tabular}
\end{center}

\begin{remark}[Formal Verification]
For $n = 3$: the existence of $\eta \approx 1.839287$ as the unique real root
of $x^3 - x^2 - x - 1 = 0$ in $[1,2]$, together with $\eta > 1$ and the
strict antitonicity of the weight sequence $\{\eta^{-k}\}_{k \geq 0}$, is
formally proved in Lean~4/Mathlib4 without \texttt{sorry} in
\texttt{TribonacciMeasure.lean} of the AXLE repository \cite{AXLE2026}.
\end{remark}

\subsection{Substitution Chains}

The Fibonacci chain ($n=2$) is generated by $\sigma_2: A \mapsto AB,\ B \mapsto A$,
producing the one-sided word $ABAABABAABAAB\cdots$. The tribonacci (Rauzy)
chain ($n=3$) is generated by $\sigma_3: A \mapsto AB,\ B \mapsto AC,\
C \mapsto A$ \cite{Krebbekx2023,Rauzy1982}.

\paragraph{Natural lengths.} The natural Fibonacci lengths $F_n$ and natural
Rauzy lengths $T_n^{(3)}$ correspond to the lengths obtained by iterating
the substitution starting from the seed letter. They are the relevant
scaling sequence for finite-size analysis: at any other length, the chain
truncation breaks the self-similar boundary structure that the substitution
preserves.

\begin{figure}[h]
\centering
\includegraphics[width=\linewidth]{fig5_substitution_tree}
\caption{Substitution trees for the Fibonacci (left, $n=2$) and
Rauzy--Tribonacci (right, $n=3$) rules across five generations. The
tribonacci tree has three letters and a richer branching structure,
producing more strongly hierarchical spatial distributions of hopping
amplitudes.}
\label{fig:substitution}
\end{figure}

\subsection{Tight-Binding Hamiltonian}

We assign hopping amplitudes via a modulation parameter $t \in (0,1)$:
$t_A = 1.0$, $t_B = t$, $t_C = t^2$. We use $t = 0.5$ throughout, in the
generic incommensurate regime studied in \cite{Krebbekx2023}. The
Hamiltonian is the tridiagonal matrix
\[
  H_{j,j+1} = H_{j+1,j} = t_{\sigma(j)}, \quad H_{jj} = 0,
\]
where $\sigma(j)$ is the $j$-th letter of the substitution word.

\begin{figure}[h]
\centering
\includegraphics[width=\linewidth]{fig1_chain_structure}
\caption{First 32 sites of the Fibonacci (top) and Rauzy--Tribonacci (bottom)
tight-binding chains. Site colour indicates letter identity. Bond thickness
is proportional to hopping amplitude.}
\label{fig:chains}
\end{figure}

\section{The DNLS Model}

\subsection{Equation of Motion}

We study the discrete nonlinear Schr\"{o}dinger equation on the
$n$-bonacci substitution chain:
\begin{equation}
  i\,\frac{d\psi_j}{dt} = -\sum_{\langle j,k\rangle} t_{jk}\,\psi_k
    + \lambda\,|\psi_j|^2\,\psi_j,
  \label{eq:DNLS}
\end{equation}
where $t_{jk}$ are the hopping amplitudes, $\lambda > 0$ is the nonlinear
coupling, and the sum is over nearest neighbors. This is the standard
focusing cubic DNLS \cite{Eilbeck1985,Kevrekidis2009}.

\subsection{Amplitude Envelope Ansatz}

The mid-gap eigenstate $\psi^{(0)}$ of the linear tribonacci Hamiltonian
serves as our initial condition. Within the dm$^3$ framework
\cite{Grossi2026_Vol1}, the natural amplitude envelope is the geometric
weight $w_k = \eta^{-k}$, whose strict antitonicity is Lean~4-verified in
\texttt{TribonacciMeasure.lean} \cite{AXLE2026}.

\subsection{Numerical Implementation}

Equation~\eqref{eq:DNLS} conserves the $\ell^2$ norm. We split
$\psi_j = u_j + iv_j$ and integrate the real $2N$-dimensional system. The
$T=50$ baseline runs (Section~4.1--4.2) use \texttt{scipy.integrate.solve\_ivp}
with the RK45 method, relative tolerance $10^{-7}$, absolute tolerance
$10^{-9}$. The long-time runs at $T \geq 10^4$ (Section~4.4 onwards) use the
DOP853 integrator (8th-order Dormand--Prince) with relative tolerance
$10^{-9}$ and absolute tolerance $10^{-11}$ to suppress integration drift
over six time decades. Norm conservation was verified to hold to
$<10^{-5}$ throughout all long-time runs. We note that absolute IPR ratio
values at $T \geq 10^5$ are reproducible across independent runs to
approximately $\pm 5\%$, consistent with the run-to-run variation observed
at the chosen tolerances; quoted ratios should be interpreted with this
relative uncertainty in mind.

\paragraph{Initial conditions.} For each chain, we compute the eigenstate
$\psi^{(0)}$ of the linear Hamiltonian with eigenvalue closest to $E=0$.
At $N=500$: tribonacci mid-gap eigenvalue $=0.000000$, IPR $=0.0820$;
Fibonacci mid-gap eigenvalue $=0.000915$, IPR $=0.0210$.

\paragraph{Localization measure.} We track the normalized IPR
$\mathrm{IPR}(t) = \sum_j |\psi_j(t)|^4 / (\sum_j |\psi_j(t)|^2)^2$.

\section{Results}

\subsection{Linear Baseline and the $T=50$ Differential}

In the linear limit ($\lambda=0$), the IPR of the mid-gap state is fixed by
the eigenstate structure: $\mathrm{IPR}_{\rm trib}(0) = 0.0820$,
$\mathrm{IPR}_{\rm fib}(0) = 0.0210$, ratio $\approx 3.91$.

\begin{figure}[h]
\centering
\includegraphics[width=\linewidth]{fig2_eigenstates}
\caption{Mid-gap critical eigenstates $|\psi_j|^2$ at $\lambda=0$, $N=500$,
$t=0.5$. The tribonacci state (IPR $=0.0820$) is substantially more
spatially concentrated than the Fibonacci state (IPR $=0.0210$).}
\label{fig:eigenstates}
\end{figure}

Table~\ref{tab:ipr} and Figures~\ref{fig:ipr}--\ref{fig:ratio}
show the IPR at time $T=50$ as a function of $\lambda$ at $N=500$.

\begin{table}[h]
\centering
\caption{IPR at $T=50$ for Fibonacci ($n=2$) and Tribonacci ($n=3$) chains,
$N=500$, $t=0.5$.}
\label{tab:ipr}
\begin{tabular}{r r r r}
\toprule
$\lambda$ & IPR$_{\rm trib}$ & IPR$_{\rm fib}$ & Ratio \\
\midrule
0.0  & 0.0820 & 0.0210 & 3.91 \\
0.5  & 0.0812 & 0.0168 & 4.82 \\
1.0  & 0.0789 & 0.0116 & 6.81 \\
1.5  & 0.0782 & 0.0090 & 8.65 \\
2.0  & 0.0704 & 0.0088 & 8.00 \\
3.0  & 0.0835 & 0.0115 & 7.28 \\
4.0  & 0.0587 & 0.0153 & 3.83 \\
5.0  & 0.0430 & 0.0173 & 2.48 \\
7.0  & 0.0316 & 0.0132 & 2.40 \\
10.0 & 0.0271 & 0.0082 & 3.30 \\
\bottomrule
\end{tabular}
\end{table}

\begin{figure}[h]
\centering
\includegraphics[width=0.85\linewidth]{fig3_ipr_vs_lambda}
\caption{IPR at $T=50$ as a function of $\lambda$. The Fibonacci mid-gap
state delocalizes rapidly ($\sim 57\%$ IPR drop at $\lambda=1.5$); the
tribonacci state remains nearly pinned ($<5\%$ drop at the same coupling).}
\label{fig:ipr}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.85\linewidth]{fig4_ipr_ratio}
\caption{Ratio IPR$_{\rm trib}$ / IPR$_{\rm fib}$ at $T=50$ as a function
of $\lambda$. The ratio increases to $\sim 8.6$ at $\lambda=1.5$ before
declining at large coupling.}
\label{fig:ratio}
\end{figure}

\subsection{Mechanism: Multifractal Hierarchy}

The mechanism underlying differential nonlinear robustness is the
hierarchical spatial structure of tribonacci mid-gap eigenstates. The
nonlinear self-trapping transition occurs when the nonlinear energy
$\lambda \langle |\psi|^4 \rangle = \lambda \cdot \mathrm{IPR} \cdot
\mathcal{N}^2$ is comparable to the bandwidth. The higher initial IPR of
the tribonacci state means it enters self-trapping at lower effective
$\lambda$, but because it is already more strongly localized, the spreading
has less work to do. The Fibonacci state, more extended, requires more
nonlinear energy to localize but loses its initial localization more
easily under perturbation.

\subsection{Multifractal Dimension at Natural Substitution Lengths}

We extract the multifractal dimension $D_2$ of the mid-gap state from
the linear-limit scaling
$\mathrm{IPR}(N) \sim N^{-D_2}$, evaluated at natural substitution lengths
to preserve the self-similar boundary structure.

\begin{figure}[h]
\centering
\includegraphics[width=\linewidth]{fig6_d2_natural}
\caption{Multifractal scaling of mid-gap IPR at natural substitution
lengths. \textbf{(a)} Fibonacci at $F_n \in \{233, 377, 610, 987, 1597,
2584\}$: clean power law, $D_2^{\rm fib} = 0.65 \pm 0.11$, $R^2=0.90$.
\textbf{(b)} Tribonacci at $T_n^{(3)} \in \{274, 504, 927, 1705, 3136\}$:
plateau across iterations $n=9,10,11$ (IPR $\approx 0.082$ across a
$3.4\times$ range in $N$) followed by a discontinuous drop to IPR $=0.041$
at $n=12$. The closest-to-$E=0$ mid-gap state selector does not consistently
track the same physical eigenstate across substitution iterations; we
defer numerical extraction of $D_2^{\rm trib}$ pending state-tracking
diagnostics.}
\label{fig:d2}
\end{figure}

The Fibonacci fit, restricted to natural lengths, gives
$D_2^{\rm fib} = 0.65 \pm 0.11$ ($R^2=0.90$, $n=6$ points). This is
consistent with multifractal critical scaling and provides a quantitative
anchor for the linear-limit multifractality.

The tribonacci fit, by contrast, exhibits structure beyond a single power
law. The IPR at $T_n^{(3)} \in \{504, 927, 1705\}$ is constant to four
significant figures (IPR $=0.0820$ at all three), then drops by a factor
of two at $T_n^{(3)} = 3136$. We interpret this as a state-isolation
artifact: at four neighboring iterations the algorithm picks states whose
IPRs coincide (likely a particular eigenstate family that recurs across
iterations of the substitution), and at the next iteration a different
state becomes the closest to $E=0$. The qualitative observation that
tribonacci mid-gap states are more localized than Fibonacci states at
comparable $N$ (IPR $\approx 0.082$ versus IPR $\approx 0.045$ at
$N \approx 500$) is robust and supports the differential-robustness
mechanism, but a clean numerical $D_2^{\rm trib}$ requires either a
state-tracking eigenmode selector via spatial-profile overlap, or a
different multifractal observable (e.g.\ box-counting of level density
near $E=0$). We list this as an open problem in Section~7.

\subsection{Finite-Size Scaling at $T=10^4$}

To test whether the differential persists or shrinks with system size, we
extend the integration to $T=10^4$ across $N \in \{500, 1000, 2000\}$ at
seven values of $\lambda$. Figure~\ref{fig:fss} shows the differential
ratio $\mathrm{IPR}_{\rm trib}/\mathrm{IPR}_{\rm fib}$ as a grouped bar
chart by $\lambda$ and $N$.

\begin{figure}[h]
\centering
\includegraphics[width=\linewidth]{fig7_fss_T1e4}
\caption{Finite-size scaling of the differential ratio
$\mathrm{IPR}_{\rm trib}/\mathrm{IPR}_{\rm fib}$ at $T=10^4$, across
$N \in \{500, 1000, 2000\}$ and $\lambda \in \{0.5, 1.0, 1.5, 2.0, 3.0,
4.0, 5.0\}$. The ratio is greater than unity at every $(\lambda,N)$ pair,
indicating tribonacci is more retained throughout. The differential
generally grows with $N$ at fixed $\lambda$ (e.g.\ at $\lambda=1.0$ the
ratio increases from $1.6\times$ at $N=500$ to $6.3\times$ at $N=2000$).
The $\lambda=1.5$ row exhibits a non-monotone $N$-dependence with a peak
at $N=1000$; the $T=10^5$ verification run (Section~4.5) shows this peak
is a transient feature.}
\label{fig:fss}
\end{figure}

The non-monotone $\lambda=1.5$ profile (1.52, 4.18, 3.17 across
$N=500,1000,2000$) is a real feature of the dynamics rather than numerical
noise; norm conservation is clean across all three cells (max norm
deviation $<4 \times 10^{-7}$). Section~4.5 demonstrates that this peak is
transient: at $T=10^5$ the same parameter slice gives a monotone ratio
profile (1.20, 1.20, 1.38) consistent with a single asymptotic trend.

The qualitative finding is that the differential generally grows with $N$,
supporting persistence of the differential in the thermodynamic limit. Two
non-monotone cells ($\lambda=4$ and $\lambda=5$) at intermediate $N=1000$
are not yet verified at $T=10^5$ and may also be transient; we list these
as open in Section~7.

\subsection{Long-Time Evolution and Ratio Collapse}

To probe the asymptotic behavior, we extend the integration at the
canonical parameter $\lambda=1.5$ first to $T=10^5$ across $N \in \{1000,
2000\}$, then to $T=10^6$ at $N=1000$.
Figure~\ref{fig:collapse} compares the differential at $T=10^4$ versus
$T=10^5$, and Figure~\ref{fig:saturation} shows the longer-time evolution
at $N=1000$ through $T=10^6$.

\begin{figure}[h]
\centering
\includegraphics[width=0.85\linewidth]{fig8_ratio_collapse}
\caption{Collapse of the differential ratio from $T=10^4$ to $T=10^5$ at
$\lambda=1.5$. The non-monotone $T=10^4$ peak at $N=1000$ ($4.18\times$)
collapses to $1.20\times$ at $T=10^5$, while the differential at $N=2000$
shrinks more modestly ($3.17\times \to 1.38\times$). At $T=10^5$ the ratio
is monotone in $N$ (consistent with a single asymptotic trend), confirming
that the $T=10^4$ non-monotonicity at $N=1000$ is a transient finite-time
feature.}
\label{fig:collapse}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.95\linewidth]{fig9_T1e6_saturation}
\caption{Long-time evolution of the differential ratio at $N=1000$,
$\lambda=1.5$. Anchor points are direct measurements at
$T \in \{10^4, 10^5, 3\times10^5, 10^6\}$. Both chains saturate to a
finite-size limited IPR ($\approx 0.002$) by $T \sim 3\times10^5$, after
which the ratio oscillates around $1.04 \pm 0.04$ (shaded band). The
differential dissolves into the saturation band at $N=1000$; at larger
$N$ the saturation time is later, and the FSS at $T=10^4$
(Figure~\ref{fig:fss}) shows the differential growing with $N$, consistent
with persistence of the differential in the thermodynamic limit.}
\label{fig:saturation}
\end{figure}

The full long-time picture at $N=1000$, $\lambda=1.5$ is therefore: at
$T=10^4$ the ratio peaks at $4.18\times$; at $T=10^5$ it has collapsed to
$1.20\times$; from $T \sim 3\times10^5$ onward both chains have saturated
to chain-limited IPR $\approx 0.002$ and the ratio settles into an
oscillatory band around $1.04\pm0.04$. The differential at finite $N$ is
ultimately bounded above by the inverse chain length, and the saturation
time scales with $N$, so at larger $N$ both the magnitude and the duration
of the differential should be larger. The FSS at $T=10^4$
(Figure~\ref{fig:fss}) is consistent with this picture.

\subsection{Pre-Saturation Spreading Rates}

In the pre-saturation regime, IPR$(t)$ decreases approximately as a power
law, $\mathrm{IPR}(t) \sim t^{-\alpha}$, with $\alpha$ a chain- and
$\lambda$-dependent rate. Fitting this scaling in the window $t > 10^4$
of the $T=10^5$ runs at $N=2000$ gives the values in
Table~\ref{tab:alpha}.

\begin{table}[h]
\centering
\caption{Pre-saturation spreading rates $\alpha$ fitted from
$\mathrm{IPR}(t) \sim t^{-\alpha}$ in the window $t > 10^4$ at $T=10^5$,
$N=2000$. Standard errors and $R^2$ from OLS on $\log_{10}$IPR vs
$\log_{10}t$.}
\label{tab:alpha}
\begin{tabular}{r r r r r}
\toprule
$\lambda$ & $\alpha_{\rm fib}$ & $R^2_{\rm fib}$ & $\alpha_{\rm trib}$ & $R^2_{\rm trib}$ \\
\midrule
0.5  & $0.11 \pm 0.02$ & 0.33 & $0.78 \pm 0.02$ & 0.94 \\
1.0  & $0.38 \pm 0.02$ & 0.84 & $0.99 \pm 0.03$ & 0.93 \\
1.5  & $0.49 \pm 0.01$ & 0.93 & $0.78 \pm 0.02$ & 0.95 \\
2.0  & $0.32 \pm 0.01$ & 0.95 & $0.67 \pm 0.01$ & 0.95 \\
\bottomrule
\end{tabular}
\end{table}

The result is uniform: $\alpha_{\rm trib} > \alpha_{\rm fib}$ at every
$\lambda$ in the table. We emphasize that $\alpha$ here is a
\emph{pre-saturation} rate of approach to the chain-limited IPR, not an
asymptotic exponent of an unbounded multifractal: at finite $N$ the IPR
ultimately saturates (Section~4.5), and $\alpha$ characterizes the
approach to that saturation rather than an indefinite power-law decay.
The $\lambda=0.5$ Fibonacci row has $R^2=0.33$, which is too low to
support a clean power-law interpretation; we include it for completeness
but do not interpret it physically.

In the thermodynamic limit ($N \to \infty$), saturation is pushed to
arbitrarily late times and the spreading regime is asymptotic. The
inequality $\alpha_{\rm trib} > \alpha_{\rm fib}$ then implies that the
differential at infinite system size shrinks with time but does not
saturate; whether it crosses below unity at any finite time is not
directly determined by the data presented here, and would require either
an analytical argument or integration in a regime where saturation is
demonstrably absent.

\subsection{Self-Trapping Threshold Gap}
\label{sec:self_trapping}

The spreading-exponent data in Table~\ref{tab:alpha} covers
$\lambda \in [0.5, 2.0]$, the pre-saturation window at $N=2000$.
Extending the fit to the full sweep $\lambda \in \{1, 2, 4, 8, 10\}$
at $N=500$, $T=10^5$ reveals a qualitative transition that the short-time
($T=50$) IPR table does not resolve: the two chains cross distinct
self-trapping thresholds at different values of $\lambda$.

\begin{table}[h]
\centering
\caption{Late-time spreading exponents $\alpha$ from
$\mathrm{IPR}(t) \sim t^{-\alpha}$ fitted over the full $T=10^5$ run at
$N=500$. At $\lambda \geq 4$ the Fibonacci exponent drops to near zero
(self-trapped); the tribonacci exponent remains $\gtrsim 0.1$ through
$\lambda = 10$.}
\label{tab:alpha_full}
\begin{tabular}{r r r l}
\toprule
$\lambda$ & $\alpha_{\rm fib}$ & $\alpha_{\rm trib}$ & Regime \\
\midrule
1  & 0.249 & 0.421 & Both spreading; tribonacci faster \\
2  & 0.119 & 0.243 & Both slowing; tribonacci $\approx 2\times$ \\
4  & 0.040 & 0.184 & Fibonacci near-frozen; tribonacci spreading \\
8  & 0.007 & 0.110 & Fibonacci self-trapped; tribonacci spreading \\
10 & 0.021 & 0.090 & Both suppressed; tribonacci still spreading \\
\bottomrule
\end{tabular}
\end{table}

The data expose a \emph{self-trapping threshold gap}: the Fibonacci chain
crosses into the self-trapped regime ($\alpha \approx 0$) in the interval
$\lambda \in [4, 8]$, while the tribonacci chain has not crossed it by
$\lambda = 10$. This is a qualitative distinction, not merely a
quantitative one.

The short-time result (Section~4.1) reported $\sim\!57\%$ IPR drop for
Fibonacci versus $<\!5\%$ for tribonacci at $\lambda=1.5$, $T=50$: a
large instantaneous differential with both chains still in the spreading
regime. The long-time result is stronger in kind: at $\lambda = 8$ the
Fibonacci chain has essentially ceased spreading
($\alpha_{\rm fib} = 0.007$), while the tribonacci chain is still
propagating the initial perturbation outward
($\alpha_{\rm trib} = 0.110$). The two observations are consistent and
reinforce each other---the $T=50$ snapshot is the early-time face of the
same multifractal resistance that, at long times, shifts the self-trapping
threshold to higher $\lambda$.

\paragraph{Mechanism.}
The standard self-trapping condition for a DNLS site is
$\lambda \langle |\psi_j|^4 \rangle \sim J_{\rm eff}$, where $J_{\rm eff}$
is the effective bandwidth \cite{Eilbeck1985,Kevrekidis2009}.
For a critical multifractal eigenstate the relevant bandwidth is set by
the gap scale $\sim \eta^{-k}$ (tribonacci) or $\sim \varphi^{-k}$
(Fibonacci) at the hierarchical level $k$ corresponding to the
localization length. Because $\eta > \varphi$, the tribonacci gap scale
decays more slowly, providing a higher effective barrier to self-trapping.
This is the same spectral-gap lever identified in the companion criticality
study~\cite{Grossi2026Criticality}, where the critical fission strength
$\lambda_c(n)$ is governed by the spectral gap $\Delta_n = \rho_n -
|\rho_n^{(2)}|$ of the substitution transfer matrix. The present result is
the nonlinear-dynamics counterpart: a higher spectral gap lifts the
self-trapping threshold.

\paragraph{Limitations.}
The $\alpha$ values in Table~\ref{tab:alpha_full} are fitted over a finite
time window at finite $N=500$; saturation effects are not fully excluded,
and the $\lambda=10$ Fibonacci uptick ($\alpha_{\rm fib} = 0.021 >
\alpha_{\rm fib}|_{\lambda=8} = 0.007$) may reflect a secondary oscillatory
mode rather than genuine residual spreading. We therefore interpret the
self-trapping threshold as lying in $\lambda \in [4, 8]$ for Fibonacci and
$\lambda > 10$ for tribonacci, with the exact values to be resolved by
a finer $\lambda$-scan (listed as an open question in Section~7).

\section{Connection to the dm$^3$ Framework}

The amplitude decay envelope $w_k = \eta^{-k}$ that appears in the dm$^3$
contact-geometric framework \cite{Grossi2026_Vol1,Grossi2026_Wavenumber6}
is identical to the Perron--Frobenius weight of the tribonacci companion
matrix. The present paper does not require the full dm$^3$ framework; the
substitution chain and DNLS model stand independently. The formal
verification of $\eta > 1$ and $w_k \searrow 0$ in
\texttt{TribonacciMeasure.lean} \cite{AXLE2026} provides Lean~4-verified
support for the mathematical well-posedness of the amplitude envelope used
as initial condition.

\section{Related Work and Gap Claim}

The closest existing work is \cite{Lahini2009}, who studied DNLS dynamics
on Fibonacci chains experimentally in photonic lattices and found
subdiffusive spreading governed by the multifractal exponent.
Ref.~\cite{Krebbekx2023} established the linear tribonacci tight-binding
model and its spectral properties but did not consider nonlinear dynamics.
To our knowledge, no prior work has studied DNLS or Gross--Pitaevskii
dynamics on a tribonacci or higher-$n$-bonacci substitution chain, nor
performed finite-size scaling or long-time saturation analysis comparing
two members of the $n$-bonacci family. The gap claim is narrow and
specific.

\section{Open Questions and Future Work}

The v1 paper listed five open questions; the present revision addresses
three of them. We summarize the current status:

\begin{enumerate}

\item \textbf{Longer time evolution} (\emph{partially closed}). Section~4.5
extends to $T=10^6$ at $N=1000$, $\lambda=1.5$. The tribonacci advantage
persists across $T \in [50, 10^5]$ but dissolves into a finite-size
saturation band by $T \sim 3\times10^5$ at $N=1000$. The thermodynamic-
limit behavior remains an open analytical question.

\item \textbf{Finite-size scaling} (\emph{closed}). Section~4.4 establishes
the FSS at $T=10^4$ across $N \in \{500, 1000, 2000\}$. The differential
generally grows with $N$. Two non-monotone cells ($\lambda=4$, $\lambda=5$
at intermediate $N=1000$) remain to be verified at $T=10^5$.

\item \textbf{Spreading exponent} (\emph{closed, with reinterpretation}).
Section~4.6 fits $\alpha$ in the pre-saturation window. The result
$\alpha_{\rm trib} > \alpha_{\rm fib}$ uniformly in $\lambda$ characterizes
the rate of approach to finite-size saturation. The direct comparison to
the Fibonacci spreading exponent of \cite{Lahini2009} requires controlling
for finite-size saturation in their setting as well; we leave this to
future experimental comparison.

\item \textbf{Self-trapping threshold} (\emph{partially closed}).
Section~\ref{sec:self_trapping} establishes that the Fibonacci
self-trapping threshold lies in $\lambda \in [4, 8]$ and the
tribonacci threshold is $> 10$ at $N=500$, $T=10^5$.
A fine $\lambda$-scan in the interval $[2, 12]$ at fixed $N=500$,
$T=10^5$ would bracket both thresholds precisely. Preliminary
tetrabonacci ($n=4$) results suggest a further upward shift consistent
with the spectral-gap-controls-threshold picture.

\item \textbf{Lean~4 formalization} (\emph{open, ongoing}). The IPR
inequality $\mathrm{IPR}_{\rm trib}(0) > \mathrm{IPR}_{\rm fib}(0)$ and a
Lean~4 statement of the differential nonlinear robustness threshold remain
on the AXLE sorry roadmap \cite{AXLE2026}.

\item \textbf{$D_2^{\rm trib}$ extraction} (\emph{new open question}).
Section~4.3 documents a state-isolation artifact in the closest-to-$E=0$
mid-gap selector at natural Rauzy lengths (plateau across $n=9,10,11$
followed by a discontinuous drop). Resolving this requires either an
eigenmode tracker via spatial-profile overlap or a different multifractal
observable.

\end{enumerate}

\section{Conclusion}

We have shown numerically that the tribonacci substitution chain ($n=3$)
exhibits substantially greater resistance to nonlinear delocalization than
the Fibonacci chain ($n=2$) when both are driven by a cubic DNLS
nonlinearity starting from mid-gap critical states. At the canonical
parameters of the original report ($T=50$, $N=500$, $\lambda=1.5$) the
tribonacci mid-gap state retains $>95\%$ of its linear IPR while the
Fibonacci state has lost $\sim 57\%$, a $\sim 8.6\times$ ratio. We term
this \emph{differential nonlinear robustness} and identify its mechanism
as the stronger effective spatial hierarchy of the tribonacci multifractal
eigenstate.

The present revision extends the analysis to (i) finite-size scaling at
$T=10^4$ across $N \in \{500, 1000, 2000\}$, where the differential
generally grows with $N$, supporting persistence in the thermodynamic
limit; (ii) long-time evolution at $N=1000$ to $T=10^6$, where both chains
saturate to a finite-size limited IPR by $T \sim 3\times10^5$ and the
ratio settles into an oscillatory band around $1.04\pm0.04$, dissolving
the differential at this $N$; (iii) pre-saturation spreading rates, where
$\alpha_{\rm trib} > \alpha_{\rm fib}$ uniformly in $\lambda \in [0.5,2]$,
characterizing the rate of approach to saturation; and (iv) multifractal
dimension extraction at natural substitution lengths, where the Fibonacci
fit gives $D_2^{\rm fib} = 0.65 \pm 0.11$ ($R^2=0.90$) and the tribonacci
fit exhibits a state-isolation plateau-and-drop pattern that we leave to
future work.

The corrected long-time narrative is that tribonacci retains higher IPR
than Fibonacci across the entire studied range $T \in [50, 10^6]$, with
the differential growing with $N$ and shrinking with $t$, ultimately
bounded above by the inverse chain length at finite $N$. The
spreading-rate inequality $\alpha_{\rm trib} > \alpha_{\rm fib}$ provides
the asymptotic-regime mechanism: tribonacci spreads faster towards
saturation, but starts from a higher initial IPR and remains more retained
throughout the studied range.

The tribonacci amplitude decay constant $\eta \approx 1.839287$ is formally
verified to satisfy $\eta > 1$ and to produce a strictly antitone weight
sequence in Lean~4/Mathlib4 \cite{AXLE2026}. To our knowledge, this is the
first numerical study of DNLS dynamics on a tribonacci substitution chain
incorporating finite-size scaling, long-time saturation, and pre-saturation
spreading rates; it identifies a concrete, measurable difference between
$n=2$ and $n=3$ quasicrystal models under nonlinear perturbation.

The self-trapping analysis (Section~\ref{sec:self_trapping}) adds a
third, qualitative layer to the differential-robustness picture.
At $T=50$ the Fibonacci state has lost $\sim\!57\%$ of its linear IPR
while the tribonacci state retains $>95\%$ (differential in magnitude).
At $T=10^5$ the Fibonacci spreading exponent collapses to
$\alpha \approx 0$ for $\lambda \geq 8$ while the tribonacci exponent
remains $\gtrsim 0.1$ through $\lambda=10$ (differential in regime).
These are two faces of the same mechanism: the tribonacci spectral gap
$\Delta_3 = \eta - |\eta^{(2)}| > \Delta_2 = \varphi - |\varphi^{(2)}|$
raises the effective self-trapping barrier.  The companion criticality
paper~\cite{Grossi2026Criticality} shows the identical lever operating
in a classical neutron-diffusion setting, where $\lambda_c(n) \approx
0.958\,\Delta_n + 0.107$ across $n = 2, \ldots, 5$ with $r = 0.989$.
The convergence of two independent physical models on the same
spectral-gap control parameter constitutes the primary finding of this
pair of papers.

\section*{Acknowledgements}

The author thanks the Lean~4/Mathlib4 community for the formal verification
infrastructure. Numerical computations used Python~3 with NumPy, SciPy, and
matplotlib. All long-time simulation code, raw output CSVs, and figure
generators are openly available at
\href{https://github.com/grossi-ops/Atratores}{github.com/grossi-ops/Atratores}.

\section*{Data and Code Availability}

All simulation scripts, raw IPR time series CSVs, analysis notebooks, and
figure generators for this revision are publicly available at
\href{https://github.com/grossi-ops/Atratores}{github.com/grossi-ops/Atratores}.
The figures in this paper were produced by
\texttt{generate\_figures\_v3.py}, embedding the data tables verbatim from
the long-time runs reported here.

\begin{thebibliography}{99}

\bibitem{Kohmoto1983}
M.~Kohmoto, L.~P.~Kadanoff, and C.~Tang,
``Localization Problem in One Dimension: Mapping and Escape,''
\textit{Phys.\ Rev.\ Lett.} \textbf{50}, 1870 (1983).

\bibitem{Ostlund1983}
S.~\"{O}stlund, R.~Pandit, D.~Rand, H.~J.~Schellnhuber, and E.~D.~Siggia,
``One-Dimensional Schr\"{o}dinger Equation with an Almost Periodic Potential,''
\textit{Phys.\ Rev.\ Lett.} \textbf{50}, 1873 (1983).

\bibitem{Luck1989}
J.~M.~Luck,
``Cantor spectra and scaling of gap widths in deterministic aperiodic systems,''
\textit{Phys.\ Rev.\ B} \textbf{39}, 5834 (1989).

\bibitem{Mace2017}
N.~Mac\'{e}, A.~Jagannathan, and F.~Pi\'{e}chon,
``Fractal dimensions of wave functions and local spectral measures
on the Fibonacci chain,''
\textit{Phys.\ Rev.\ B} \textbf{93}, 205134 (2016).

\bibitem{Kraus2012}
Y.~E.~Kraus, Y.~Lahini, Z.~Ringel, M.~Verbin, and O.~Zilberberg,
``Topological States and Adiabatic Pumping in Quasicrystals,''
\textit{Phys.\ Rev.\ Lett.} \textbf{109}, 106402 (2012).

\bibitem{Lang2012}
L.-J.~Lang, X.~Cai, and S.~Chen,
``Edge States and Topological Phases in One-Dimensional Optical Superlattices,''
\textit{Phys.\ Rev.\ Lett.} \textbf{108}, 220401 (2012).

\bibitem{Krebbekx2023}
N.~J.~Krebbekx, A.~Moustaj, K.~Dajani, and C.~Morais~Smith,
``Multifractal properties of tribonacci chains,''
\textit{Phys.\ Rev.\ B} \textbf{108}, 104204 (2023).

\bibitem{Bellissard1989}
J.~Bellissard, B.~Iochum, E.~Scoppola, and D.~Testard,
``Spectral properties of one dimensional quasi-crystals,''
\textit{Commun.\ Math.\ Phys.} \textbf{125}, 527 (1989).

\bibitem{Damanik2008}
D.~Damanik and A.~Gorodetski,
``Hyperbolicity of the trace map for the weakly coupled Fibonacci Hamiltonian,''
\textit{Nonlinearity} \textbf{22}, 123 (2009).

\bibitem{Rauzy1982}
G.~Rauzy,
``Nombres alg\'{e}briques et substitutions,''
\textit{Bull.\ Soc.\ Math.\ France} \textbf{110}, 147 (1982).

\bibitem{Lahini2009}
Y.~Lahini, R.~Pugatch, F.~Pozzi, M.~Sorel, R.~Morandotti, N.~Davidson,
and Y.~Silberberg,
``Observation of a Localization Transition in Quasiperiodic Photonic Lattices,''
\textit{Phys.\ Rev.\ Lett.} \textbf{103}, 013901 (2009).

\bibitem{Jurgensen2021}
M.~J\"{u}rgensen, S.~Mukherjee, and M.~C.~Rechtsman,
``Quantized nonlinear Thouless pumping,''
\textit{Nature} \textbf{596}, 63 (2021).

\bibitem{Verbin2015}
M.~Verbin, O.~Zilberberg, Y.~Lahini, Y.~E.~Kraus, and Y.~Silberberg,
``Topological pumping over a photonic Fibonacci quasicrystal,''
\textit{Phys.\ Rev.\ B} \textbf{91}, 064201 (2015).

\bibitem{Eilbeck1985}
J.~C.~Eilbeck, P.~S.~Lomdahl, and A.~C.~Scott,
``The discrete self-trapping equation,''
\textit{Physica D} \textbf{16}, 318 (1985).

\bibitem{Kevrekidis2009}
P.~G.~Kevrekidis,
\textit{The Discrete Nonlinear Schr\"{o}dinger Equation}.
Springer, Berlin (2009).

\bibitem{Flach2010}
S.~Flach, M.~V.~Ivanchenko, and O.~I.~Kanakov,
``Spreading of wave packets in one-dimensional disordered lattices,''
\textit{Phys.\ Rev.\ E} \textbf{82}, 036219 (2010).

\bibitem{Grossi2026_DNLS}
P.~Nogueira~Grossi,
``Differential Nonlinear Robustness of Critical States in Fibonacci and
Tribonacci Substitution Chains.''
G6~LLC. Zenodo (concept DOI):
\href{https://doi.org/10.5281/zenodo.20026942}{10.5281/zenodo.20026942}
(2026). [V4: \href{https://doi.org/10.5281/zenodo.20075822}{10.5281/zenodo.20075822}.]

\bibitem{Grossi2026Criticality}
P.~Nogueira~Grossi,
``Criticality Thresholds in One-Dimensional Multiplying Media
with $n$-Bonacci Aperiodic Modulation:
Spectral Gap Control of $k_{\rm eff}$ in
Substitution-Sequence Diffusion Operators.''
G6~LLC. Zenodo:
\href{https://doi.org/10.5281/zenodo.20077205}{10.5281/zenodo.20077205}
(2026).

\bibitem{Grossi2026_Vol1}
P.~Nogueira~Grossi,
\textit{Principia Orthogona, Volume~I: The Mathematics of Generative Transitions}.
G6~LLC. Zenodo: \href{https://doi.org/10.5281/zenodo.19117400}{10.5281/zenodo.19117400} (2026).

\bibitem{Grossi2026_Wavenumber6}
P.~Nogueira~Grossi,
``Wavenumber 6: The Orthogenetic Stability Generator of Nested Infinities.''
G6~LLC. Zenodo: \href{https://doi.org/10.5281/zenodo.19199474}{10.5281/zenodo.19199474} (2026).

\bibitem{Grossi2026_Tribonacci}
P.~Nogueira~Grossi,
``Convergent Derivations of the Tribonacci Coherence Threshold:
Contact Geometry, Rational Quantum Mechanics, and the SL(3,Z) Algebra.''
G6~LLC (2026). Manuscript.

\bibitem{AXLE2026}
P.~Nogueira~Grossi,
\textit{AXLE: Lean~4 Formal Verification Engine}.
GitHub: \href{https://github.com/TOTOGT/AXLE}{github.com/TOTOGT/AXLE} (2026).

\bibitem{AtraotresRepo}
P.~Nogueira~Grossi,
\textit{Atratores: DNLS simulation code and long-time data}.
GitHub: \href{https://github.com/grossi-ops/Atratores}{github.com/grossi-ops/Atratores} (2026).

\end{thebibliography}

\end{document}
