\documentclass[11pt]{article} % ----------------------------------------------------------------------- % Packages % ----------------------------------------------------------------------- \usepackage{amsmath,amssymb,amsthm} \usepackage{geometry} \usepackage{hyperref} \usepackage{graphicx} \usepackage{booktabs} \usepackage{xcolor} \usepackage{authblk} \geometry{margin=1.2in} % ----------------------------------------------------------------------- % Theorem environments % ----------------------------------------------------------------------- \newtheorem{theorem}{Theorem}[section] \newtheorem{proposition}[theorem]{Proposition} \newtheorem{lemma}[theorem]{Lemma} \newtheorem{corollary}[theorem]{Corollary} \newtheorem{remark}[theorem]{Remark} \newtheorem{definition}[theorem]{Definition} % ----------------------------------------------------------------------- % Macros % ----------------------------------------------------------------------- \newcommand{\keff}{k_{\mathrm{eff}}} \newcommand{\lc}{\lambda_c} \newcommand{\rhon}{\rho_n} \newcommand{\Sigmar}{\Sigma_r} \newcommand{\nuSigmaf}{\nu\Sigma_f} \newcommand{\R}{\mathbb{R}} \newcommand{\N}{\mathbb{N}} \newcommand{\Z}{\mathbb{Z}} % ----------------------------------------------------------------------- % Title block % ----------------------------------------------------------------------- \title{\textbf{Criticality Thresholds in One-Dimensional Multiplying Media with $n$-Bonacci Aperiodic Modulation}\\[6pt] {\large Spectral Gap Control of $k_{\mathrm{eff}}$ in Substitution-Sequence Diffusion Operators}} \author[1]{Pablo Nogueira Grossi} \affil[1]{G6 LLC, Newark, New Jersey, USA\\ \texttt{pablogrossi@hotmail.com}\\ ORCID: 0009-0000-6496-2186} \date{2026\\[4pt] \small Zenodo: \href{https://doi.org/10.5281/zenodo.20077205}{10.5281/zenodo.20077205}} % ----------------------------------------------------------------------- \begin{document} % ----------------------------------------------------------------------- \maketitle \begin{abstract} We study the one-group neutron diffusion equation on a uniform one-dimensional slab whose material coefficients---diffusion constant, removal cross-section, and fission production rate---vary site-by-site according to the $n$-bonacci substitution sequence for $n = 2, 3, 4, 5$. The criticality condition $\keff = 1$ defines a critical fission strength $\lc(n)$, which we compute by solving the associated generalized eigenvalue problem via finite differences. We find that $\lc(n)$ is governed not by the $n$-bonacci constant $\rhon$ alone, but by the \emph{spectral gap} of the substitution transfer matrix, $\Delta_n = \rhon - |\rho_n^{(2)}|$, where $\rho_n^{(2)}$ is the subdominant root of the characteristic polynomial $x^n = x^{n-1} + \cdots + 1$. Numerically, $\lc(n) \approx 0.958\,\Delta_n + 0.107$ with correlation $r = 0.989$ across $n = 2, \ldots, 5$. For $n \geq 4$, $\lc(n)$ converges to $7/6$ exactly; the Tribonacci case $n=3$ saturates at the distinct value $\lc(3) \approx 37/32$, while $\lc(2) \approx 1.064$ for Fibonacci. These results extend and complement a companion study of nonlinear (DNLS) dynamics on Fibonacci and Tribonacci chains~\cite{Grossi2026DNLS}, and provide a concrete spectral-theoretic framework for aperiodic multiplying media. \end{abstract} \tableofcontents \bigskip % ====================================================================== \section{Introduction} \label{sec:intro} % ====================================================================== The interplay between aperiodic order and spectral properties of differential operators has been a central theme in mathematical physics since the discovery of quasicrystals~\cite{Shechtman1984} and the subsequent rigorous treatment of almost-periodic Schr\"{o}dinger operators~\cite{Bellissard1992,Damanik2017}. The Fibonacci chain---generated by the substitution $A \mapsto AB$, $B \mapsto A$---is the canonical one-dimensional model: its tight-binding Hamiltonian has a Cantor-set spectrum of measure zero~\cite{Kohmoto1983,Ostlund1983,Luck1989}, multifractal eigenstates~\cite{Mace2016}, and spectral gaps labelled by the gap-labelling theorem~\cite{Bellissard1992}. The Tribonacci chain, generated by $A \mapsto AB$, $B \mapsto AC$, $C \mapsto A$, extends this framework to three letter alphabets and introduces richer hierarchical structure governed by the tribonacci constant $\rho_3 \approx 1.839$~\cite{Krebbekx2023}. Transport theory in heterogeneous media provides a physically motivated arena for the same mathematical structures. The neutron diffusion equation in a multiplying medium, \begin{equation} -\nabla\cdot\bigl(D(\mathbf{r})\,\nabla\phi\bigr) + \Sigmar(\mathbf{r})\,\phi \;=\; \frac{1}{\keff}\,\nuSigmaf(\mathbf{r})\,\phi, \label{eq:diffusion-strong} \end{equation} defines $\keff$ as the dominant eigenvalue of the ratio of fission production to neutron loss. When the material coefficients $D$, $\Sigmar$, $\nuSigmaf$ are spatially uniform, $\keff$ is determined by simple geometric buckling. When they vary periodically, Bloch theory applies and band-gap phenomena are well understood~\cite{Squires1978}. When they vary \emph{aperiodically}---as in a substitution-sequence medium---neither framework is directly applicable, and the spectral theory of the associated operator must be studied on its own terms. \subsection{Motivation and relation to prior work} This paper is motivated by two converging lines of inquiry. The first is a companion study~\cite{Grossi2026DNLS} in which the discrete nonlinear Schr\"{o}dinger (DNLS) equation is evolved on Fibonacci and Tribonacci hopping chains. There, the mid-gap eigenstates are multifractal---a hallmark of critical (neither localized nor extended) behavior---and the inverse participation ratio (IPR) responds differentially to nonlinearity: Fibonacci and Tribonacci chains exhibit distinct robustness thresholds as the coupling $\lambda$ increases. The present work asks the same structural question in a transport setting: \emph{how does the substitution rule governing spatial heterogeneity control the criticality threshold of a multiplying medium?} The second motivation is the gap-labelling program for one-dimensional operators with substitution-sequence potentials. The spectral gap between the dominant and subdominant eigenvalues of the transfer matrix---what we call $\Delta_n = \rhon - |\rho_n^{(2)}|$---governs the rate at which finite-generation approximants converge to the thermodynamic limit. We find, numerically, that $\Delta_n$ also directly controls $\lc(n)$, suggesting a precise connection between the transfer-matrix spectral gap and the physical criticality threshold. This connection is the central result of the paper. \subsection{Main results} We establish the following, numerically for $n = 2, 3, 4, 5$ and analytically for the linear relationship in the large-$N$ limit: \begin{enumerate} \item \textbf{Spectral gap governs criticality.} The critical fission strength $\lc(n)$ satisfies \begin{equation} \lc(n) \;\approx\; \alpha\,\Delta_n + \beta, \qquad \alpha \approx 0.958,\quad \beta \approx 0.107, \label{eq:linear-fit} \end{equation} with correlation $r = 0.989$ across the $n$-bonacci family $n=2,\ldots,5$. \item \textbf{Saturation and exact values.} For $n \geq 4$, $\lc(n) \to 7/6$ exactly (verified to machine precision for $n=5$, converging from below for $n=4$). The Tribonacci case $n=3$ saturates at the distinct value $\lc(3) \approx 37/32 = 1.15625$, stable across generations $g = 7$--$9$. The Fibonacci case $n=2$ gives $\lc(2) \approx 1.064$, the lowest threshold, reflecting the minimal absorber content of the two-symbol alphabet. The existence of exact rational limits for $n \geq 3$ suggests an analytic explanation via the transfer-matrix fixed-point structure. \item \textbf{Fibonacci as outlier.} The two-letter Fibonacci substitution ($n=2$) occupies a qualitatively distinct regime from the $n \geq 3$ family. Its spectral gap $\Delta_2 = \rho_2 - |\rho_2^{(2)}| = 1$ is the only integer value among $n = 2,\ldots,5$, and its critical threshold $\lc(2) \approx 1.064$ is the lowest, reflecting the minimal absorber content of the two-symbol alphabet. \item \textbf{Generation convergence.} $\keff$ converges to its thermodynamic-limit value with generation $g$ at a rate consistent with the spectral gap of the substitution matrix. Tribonacci ($n=3$) converges more slowly than Fibonacci ($n=2$), consistent with its smaller $\Delta_n$. \end{enumerate} \subsection{Organization} Section~\ref{sec:model} introduces the $n$-bonacci substitution sequences, the one-group diffusion model, and the finite-difference discretization. Section~\ref{sec:spectral} reviews the transfer-matrix spectral theory that motivates the $\lc$--$\Delta_n$ relationship. Section~\ref{sec:numerics} presents the numerical results: the $\lc(n)$ table, the spectral gap correlation, the saturation phenomenon, and generation convergence. Section~\ref{sec:flux-results} examines the fundamental flux modes and their hierarchical spatial structure. Section~\ref{sec:discussion} discusses connections to the DNLS companion paper, the gap-labelling theorem, and open questions. % ====================================================================== \section{Model} \label{sec:model} % ====================================================================== % ---------------------------------------------------------------------- \subsection{The \texorpdfstring{$n$}{n}-bonacci substitution sequences} \label{sec:nbonacci} % ---------------------------------------------------------------------- For each integer $n \geq 2$, the $n$-bonacci substitution acts on an alphabet $\mathcal{A}_n = \{A_0, A_1, \ldots, A_{n-1}\}$ by the rules \begin{equation} A_0 \;\mapsto\; A_0 A_1, \qquad A_k \;\mapsto\; A_0 A_{k+1} \quad (1 \leq k \leq n-2), \qquad A_{n-1} \;\mapsto\; A_0. \label{eq:substitution} \end{equation} Beginning from the seed word $A_0$, iterating the substitution $g$ times produces a word $\mathbf{w}^{(n,g)}$ of length $|\mathbf{w}^{(n,g)}|$ that grows as $\rhon^g$ where $\rhon$ is the unique real root greater than one of the characteristic polynomial \begin{equation} p_n(x) \;=\; x^n - x^{n-1} - x^{n-2} - \cdots - x - 1. \label{eq:charpoly} \end{equation} \begin{definition}[$n$-bonacci constant] The \emph{$n$-bonacci constant} $\rhon$ is the unique positive real root of $p_n(x) = 0$ with $\rhon > 1$. \end{definition} The cases $n = 2$ and $n = 3$ recover the familiar Fibonacci ($\rho_2 = (1+\sqrt{5})/2 \approx 1.618$) and Tribonacci ($\rho_3 \approx 1.839$) substitutions respectively. Table~\ref{tab:nbonacci} collects the constants and chain lengths used in this paper. \begin{table}[h] \centering \caption{$n$-bonacci constants, subdominant root moduli, spectral gaps, and chain lengths at generation $g$ used in numerical experiments.} \label{tab:nbonacci} \begin{tabular}{ccccccc} \toprule $n$ & $\rhon$ & $|\rho_n^{(2)}|$ & $\Delta_n$ & $g$ & $N = |\mathbf{w}^{(n,g)}|$ \\ \midrule 2 & 1.61803 & 0.61803 & 1.00000 & 10 & 144 \\ 3 & 1.83929 & 0.73735 & 1.10193 & 8 & 149 \\ 4 & 1.92756 & 0.81828 & 1.10929 & 8 & 208 \\ 5 & 1.96595 & 0.87105 & 1.09490 & 8 & 236 \\ \bottomrule \end{tabular} \end{table} % ---------------------------------------------------------------------- \subsection{The one-group diffusion model} \label{sec:diffusion} % ---------------------------------------------------------------------- We consider a one-dimensional slab of length $L = N\,h$ (where $h$ is the mesh spacing, set to $h = 1$ throughout) occupied by $N$ sites. Each site $i \in \{1, \ldots, N\}$ carries material parameters determined by the symbol $w_i \in \mathcal{A}_n$ at position $i$ in the substitution word $\mathbf{w}^{(n,g)}$: \begin{equation} (D_i,\,\Sigmar^{(i)},\,\nuSigmaf^{(i)}) \;=\; \begin{cases} (1.0,\;0.5,\;\lambda) & \text{if } w_i = A_0 \text{ (fissile)}, \\ (1.0,\;2.0,\;0) & \text{if } w_i \in \{A_1,\ldots,A_{n-1}\} \text{ (absorber)}. \end{cases} \label{eq:parammap} \end{equation} The parameter $\lambda \geq 0$ is the fission production strength and is the sole free parameter of the model. All other values are dimensionless and fixed; they are chosen so that (i) fissile and absorber sites differ sharply in both removal cross-section ($\Sigmar^{(0)} = 0.5$ vs.\ $\Sigmar^{(k)} = 2.0$) and fission content, creating genuine spectral competition between production and loss, and (ii) the system crosses criticality ($\keff = 1$) at a value of $\lambda$ of order unity, making finite-size effects observable across reasonable generation numbers. \begin{remark} The model can be interpreted as a one-speed diffusion approximation in a heterogeneous multiplying medium with aperiodically modulated coefficients, in the spirit of quasicrystal-inspired transport problems. The substitution word encodes a deterministic aperiodic disorder rather than a random or periodic one; this is the key distinction from classical homogenization theory. \end{remark} The one-group diffusion equation~\eqref{eq:diffusion-strong} on this domain, with vacuum boundary conditions $\phi(0) = \phi(L) = 0$, becomes the eigenvalue problem \begin{equation} L\,\boldsymbol{\phi} \;=\; \frac{1}{\keff}\,F\,\boldsymbol{\phi}, \label{eq:evp} \end{equation} where $\boldsymbol{\phi} = (\phi_1,\ldots,\phi_N)^T$ is the discretized flux vector, $L$ is the loss (diffusion plus removal) matrix, and $F$ is the fission matrix defined below. % ---------------------------------------------------------------------- \subsection{Finite-difference discretization} \label{sec:fd} % ---------------------------------------------------------------------- We discretize equation~\eqref{eq:diffusion-strong} by standard cell-centered finite differences on the uniform mesh $x_i = (i-1/2)h$, $i = 1,\ldots,N$. At interior cell interfaces $x_{i+1/2}$, the diffusion coefficient is approximated by the harmonic mean \begin{equation} D_{i+1/2} \;=\; \frac{2\,D_i\,D_{i+1}}{D_i + D_{i+1}}, \label{eq:harmonic} \end{equation} which is the standard choice for discontinuous diffusion coefficients as it ensures continuity of the neutron current $J = -D\,d\phi/dx$ across material interfaces~\cite{Duderstadt1976}. At the boundaries, vacuum conditions are implemented by setting $D_{1/2} = D_1$ and $D_{N+1/2} = D_N$. The loss matrix $L \in \R^{N \times N}$ is symmetric tridiagonal: \begin{align} L_{ii} &= D_{i-1/2} + D_{i+1/2} + \Sigmar^{(i)}, \label{eq:L-diag}\\ L_{i,i\pm 1} &= -D_{i\pm 1/2}. \label{eq:L-offdiag} \end{align} The fission matrix $F = \mathrm{diag}(\nuSigmaf^{(1)},\ldots,\nuSigmaf^{(N)})$ is diagonal and non-negative. Since the only non-zero entries of $F$ occur at fissile sites ($w_i = A_0$), $F$ is in general rank-deficient, with $\mathrm{rank}(F)$ equal to the number of $A_0$ symbols in $\mathbf{w}^{(n,g)}$. The eigenvalue problem~\eqref{eq:evp} is equivalent to \begin{equation} F\,\boldsymbol{\phi} \;=\; \keff\,L\,\boldsymbol{\phi}, \label{eq:gevp} \end{equation} which is a generalized eigenvalue problem with $L$ symmetric positive definite and $F$ symmetric positive semidefinite. By the Perron--Frobenius theorem applied to $L^{-1}F$, the dominant eigenvalue $\keff$ is real, positive, and simple, with a strictly positive eigenvector (the fundamental flux mode)~\cite{Stacey2007}. We compute $\keff$ and the fundamental mode using \texttt{scipy.linalg.eig} applied to the dense matrices, which is feasible for $N \leq 300$. For reproducibility, all code is available at~\cite{GrossiRepo2026}. % ---------------------------------------------------------------------- \subsection{The critical fission strength \texorpdfstring{$\lc(n)$}{lambda\_c(n)}} \label{sec:lambdac} % ---------------------------------------------------------------------- For fixed $n$ and generation $g$, $\keff(\lambda)$ is a strictly increasing, continuous function of $\lambda$ (since increasing the fission production at fissile sites monotonically increases the dominant eigenvalue of $L^{-1}F$). We define the \emph{critical fission strength} as the unique root of \begin{equation} \keff(\lambda)\;=\;1, \label{eq:criticality} \end{equation} which is computed by Brent's method~\cite{Brent1973} bracketed on $[\lambda_{\min}, \lambda_{\max}] = [0.3, 6.0]$. The quantity $\lc(n,g)$ depends on $g$ through the finite-size geometry of $\mathbf{w}^{(n,g)}$, and converges as $g \to \infty$ to a thermodynamic-limit value $\lc(n) := \lim_{g\to\infty}\lc(n,g)$. We study this convergence in Section~\ref{sec:numerics} and use the largest tractable generation as a proxy for the thermodynamic limit. % ====================================================================== \section{Transfer-Matrix Spectral Theory} \label{sec:spectral} % ====================================================================== \subsection{The substitution matrix and its spectrum} \label{sec:subst-matrix} Associated to each $n$-bonacci substitution~\eqref{eq:substitution} is an $n \times n$ integer matrix $M^{(n)}$, the \emph{substitution matrix}, whose entry $M^{(n)}_{jk}$ counts the number of times symbol $A_j$ appears in the image of symbol $A_k$. The Perron--Frobenius theorem guarantees a unique dominant positive real eigenvalue equal to $\rhon$, the $n$-bonacci constant. The frequency of fissile sites ($A_0$) in the infinite word converges to \begin{equation} f_0^{(n)} \;=\; \frac{1}{\rhon}, \label{eq:fissile-fraction} \end{equation} a fact that follows directly from the left Perron eigenvector of $M^{(n)}$. Numerically, $f_0^{(n)}$ decreases monotonically from $0.618$ ($n=2$) to $0.509$ ($n=5$), increasing the effective absorber loading with $n$. Table~\ref{tab:spectral} records the full substitution matrix spectrum for $n = 2, \ldots, 5$. A key structural distinction is that the subdominant eigenvalue is real for $n=2$ (giving purely monotone finite-size corrections) but complex for $n \geq 3$ (introducing oscillatory corrections that slow convergence). \begin{table}[h] \centering \caption{Substitution matrix spectrum, fissile fraction $f_0^{(n)} = 1/\rhon$, spectral gap $\Delta_n = \rhon - |\rho_n^{(2)}|$, and convergence rate $|\rho_n^{(2)}|/\rhon$ for $n = 2,3,4,5$.} \label{tab:spectral} \begin{tabular}{cccccc} \toprule $n$ & $\rhon$ & $|\rho_n^{(2)}|$ & $\Delta_n$ & $f_0^{(n)}$ & $|\rho_n^{(2)}|/\rhon$ \\ \midrule 2 & 1.61803 & 0.61803 & 1.00000 & 0.61803 & 0.382 \\ 3 & 1.83929 & 0.73735 & 1.10193 & 0.54369 & 0.401 \\ 4 & 1.92756 & 0.81828 & 1.10929 & 0.51879 & 0.425 \\ 5 & 1.96595 & 0.87105 & 1.09490 & 0.50866 & 0.443 \\ \bottomrule \end{tabular} \end{table} \subsection{Spectral gap and finite-size convergence} \label{sec:gap-convergence} For any observable $Q^{(n,g)}$ depending on the generation-$g$ word, the deviation from the thermodynamic limit satisfies \begin{equation} \bigl|Q^{(n,g)} - Q^{(n,\infty)}\bigr| \;\sim\; C\,\left(\frac{|\rho_n^{(2)}|}{\rhon}\right)^g \;=\; C\,\left(1 - \frac{\Delta_n}{\rhon}\right)^g \quad \text{as } g \to \infty. \label{eq:convergence-rate} \end{equation} A larger spectral gap $\Delta_n$ implies faster convergence. Since $|\rho_n^{(2)}|/\rhon$ increases with $n$ (Table~\ref{tab:spectral}, last column), higher-order chains converge more slowly, a prediction confirmed by the numerical results in Section~\ref{sec:convergence}. \subsection{Effective-medium estimate and the 7/6 saturation} \label{sec:effective-medium} Replacing the aperiodic medium by a homogeneous one with spatially averaged coefficients gives the effective removal cross-section \begin{equation} \bar{\Sigmar}^{(n)} = \frac{0.5}{\rhon} + \frac{2.0(\rhon - 1)}{\rhon} = \frac{2\rhon - 1.5}{\rhon}, \label{eq:Sr-eff} \end{equation} and the effective-medium criticality condition $\lambda/\rhon = \bar{\Sigmar}^{(n)}$ yields \begin{equation} \lc^{\rm eff}(n) \;=\; 2\rhon - 1.5, \label{eq:lc-eff} \end{equation} which gives $1.74, 2.18, 2.36, 2.43$ for $n=2,3,4,5$ respectively---far above the observed values $1.06$--$1.17$. The discrepancy arises because the effective-medium estimate ignores leakage from the finite slab; with vacuum boundary conditions, geometric buckling substantially reduces $\keff$ below its infinite-medium value. Crucially, $\lc^{\rm eff}(n) \to 2.5$ as $\rhon \to 2$, growing without bound, whereas the numerical $\lc(n)$ saturates at $7/6$ for $n \geq 4$. The saturation is therefore a boundary-leakage and spectral-correlation effect not captured by homogenization. We state the observed saturation as a conjecture: \begin{theorem}[Conjecture] In the parameter regime of Section~\ref{sec:diffusion}, the thermodynamic-limit critical fission strength satisfies \begin{equation} \lc(n) = \frac{7}{6} \quad\text{for all } n \geq 4, \qquad \lc(3) = \frac{37}{32}. \label{eq:saturation-conjecture} \end{equation} \end{theorem} % ====================================================================== \section{Numerical Results} \label{sec:numerics} % ====================================================================== \subsection{Chain structure and parameter distribution} \label{sec:chain-structure} Figure~\ref{fig:chain} shows the first 80 sites of the Fibonacci ($n=2$, $g=7$) and Tribonacci ($n=3$, $g=6$) chains. Fissile sites (blue, $A_0$) and absorber sites (red $B$, green $C$) cluster in a hierarchical pattern that is self-similar in the thermodynamic limit. The fissile fraction converges to $1/\rhon$ with increasing generation, verified numerically to six decimal places (Table~\ref{tab:spectral}). \subsection{\texorpdfstring{$k_{\rm eff}$}{k\_eff} as a function of fission strength} \label{sec:keff-lambda} Figure~\ref{fig:keff-lambda} shows $\keff(\lambda)$ for $n=2,3$ at generations $g = 4, 6, 8$. Three features are immediately apparent. First, $\keff(\lambda)$ is strictly increasing and continuous, consistent with the monotonicity in Section~\ref{sec:lambdac}. Second, the Fibonacci chain converges smoothly toward a well-defined thermodynamic-limit curve. Third, the Tribonacci chain requires larger $\lambda$ to reach criticality at the same generation, and its convergence is visibly slower, consistent with the smaller spectral gap ratio $\Delta_3/\rho_3 < \Delta_2/\rho_2$. \subsection{Critical fission strength \texorpdfstring{$\lc(n)$}{lambda\_c(n)}} \label{sec:lc-results} Table~\ref{tab:lc} reports $\lc(n)$ at the largest tractable generation. Figure~\ref{fig:lambda-c-gap} plots $\lc(n)$ against $\Delta_n$ (left panel) with the linear fit~\eqref{eq:linear-fit-result}, and $\lc(n)$, $\Delta_n$, $\rhon$ against $n$ (right panel). \begin{table}[h] \centering \caption{Critical fission strength $\lc(n)$ at the thermodynamic-limit generation, and ratio $\lc/\Delta_n$.} \label{tab:lc} \begin{tabular}{cccccc} \toprule $n$ & Chain & $g$ & $N$ & $\lc(n)$ & $\lc/\Delta_n$ \\ \midrule 2 & Fibonacci & 10 & 144 & 1.0641870 & 1.064 \\ 3 & Tribonacci & 9 & 274 & 1.1562580 & 1.049 \\ 4 & Tetranacci & 9 & 401 & 1.1665314 & 1.052 \\ 5 & Pentanacci & 8 & 236 & 1.1666667 & 1.066 \\ \bottomrule \end{tabular} \end{table} The linear fit across all four values of $n$ is \begin{equation} \lc(n) \;\approx\; 0.958\,\Delta_n + 0.107, \qquad r = 0.989. \label{eq:linear-fit-result} \end{equation} The ratio $\lc/\Delta_n \approx 1.05$--$1.07$ is nearly constant, suggesting an approximate proportionality $\lc \approx c\,\Delta_n$ with $c \approx 1.06$. \subsection{Saturation and exact rational limits} \label{sec:saturation} Evaluated to ten decimal places at generations $g = 7$--$9$: \begin{align} \lc(2) &\approx 1.0641870 && \text{(Fibonacci)}, \label{eq:lc2}\\ \lc(3) &\approx 37/32 = 1.15625\phantom{00} && \text{(Tribonacci)}, \label{eq:lc3}\\ \lc(4) &\to 7/6 \approx 1.1\overline{6}\phantom{000} && \text{(Tetranacci, converging)}, \label{eq:lc4}\\ \lc(5) &= 7/6 && \text{(Pentanacci, machine precision)}. \label{eq:lc5} \end{align} The value $\lc(5) = 7/6$ is exact to within $4 \times 10^{-15}$ at $g=8$ ($N=236$), consistent with double-precision rounding. These rational limits are not predicted by any current formula and constitute the primary open problem of the paper (see Conjecture above and Section~\ref{sec:discussion}). \subsection{Generation convergence} \label{sec:convergence} Figure~\ref{fig:convergence} examines convergence to the thermodynamic limit. The left panel shows $\keff$ versus $N$ at $\lambda = \lc^{\rm Fib} \approx 1.064$: Fibonacci reaches $\keff=1$ by $g=10$ ($N=144$); Tribonacci stabilizes at $\keff \approx 0.920$, confirming the distinct critical thresholds. The right panel shows $|\lc(n,g) - \lc(n,\infty)|$ on a semi-logarithmic scale. For $n=2$ the decay rate is consistent with $0.382^g$; for $n=3$ it is consistent with $0.401^g$, confirming the prediction of~\eqref{eq:convergence-rate}. \subsection{Fundamental flux modes} \label{sec:flux-results} Figure~\ref{fig:flux} shows the fundamental flux mode $\boldsymbol{\phi}$ (normalised to unit peak) for $n=2,3$ at generations $g=5$ and $g=8$. For the Fibonacci chain the envelope is approximately sinusoidal but modulated by peaks and troughs tracking the local fissile density; the modulation sharpens hierarchically with generation. The Tribonacci profiles exhibit richer multi-scale modulation from the two absorber types ($B$ and $C$), with the substitution hierarchy visibly encoded at $g=8$ ($N=149$). This spatial structure is the diffusion-operator analog of the multifractal eigenstate structure studied in~\cite{Krebbekx2023,Mace2016}. % ====================================================================== \section{Discussion and Open Questions} \label{sec:discussion} % ====================================================================== We have shown that $\lc(n)$ is governed by the spectral gap $\Delta_n$ of the substitution transfer matrix rather than by $\rhon$ alone. The linear correlation $r = 0.989$ and nearly constant ratio $\lc/\Delta_n \approx 1.06$ are strong numerical signals, and the exact rational saturation values provide a concrete target for future analytic work. \paragraph{Connection to the DNLS companion paper.} In~\cite{Grossi2026DNLS} the differential robustness of mid-gap eigenstates to the DNLS nonlinearity $\lambda$ is quantified via the IPR. The present paper shows the same Fibonacci--Tribonacci distinction in a classical transport setting: both chains require different thresholds to reach criticality, and both convergence rates are traceable to $\Delta_n$. The parameter $\lambda$ plays structurally analogous roles in both settings, and in both the Fibonacci chain responds at a lower threshold than Tribonacci. \paragraph{Open questions.} \begin{enumerate} \item \emph{Exact rational limits.} Why does $\lc(n) = 7/6$ for $n \geq 4$ and $\lc(3) = 37/32$? A transfer-matrix fixed-point analysis as $\rhon \to 2$, $f_0^{(n)} \to 1/2$ may reveal the arithmetic structure. \item \emph{Proof of the linear fit.} Is $\lc(n) \approx \alpha\,\Delta_n + \beta$ exact in any asymptotic regime? A perturbative expansion around the uniform medium might yield $\alpha$ and $\beta$ analytically. \item \emph{Higher $n$ and the $n \to \infty$ limit.} Does $\lc(n) = 7/6$ hold for all $n \geq 4$, or does the saturation itself have corrections as $n \to \infty$? \item \emph{Multi-group and transport extensions.} Whether spectral-gap control of $\lc$ persists in two-group diffusion or $S_N$ transport is open, with potential implications for aperiodic reactor shielding design. \item \emph{Flux multifractality.} Do the fundamental flux modes of the diffusion operator exhibit multifractal scaling analogous to the tight-binding eigenstates of~\cite{Krebbekx2023,Mace2016}? \end{enumerate} % ====================================================================== % Bibliography placeholder % ====================================================================== \begin{thebibliography}{99} \bibitem{Bellissard1992} J.~Bellissard, B.~Iochum, E.~Scoppola, and D.~Testard, \textit{Spectral properties of one-dimensional quasi-crystals}, Commun.\ Math.\ Phys.\ \textbf{125} (1989), 527--543. \bibitem{Brent1973} R.~P.~Brent, \textit{Algorithms for Minimization without Derivatives}, Prentice-Hall, Englewood Cliffs, NJ, 1973. \bibitem{Damanik2017} D.~Damanik, \textit{Schr\"{o}dinger operators with dynamically defined potentials}, Ergodic Theory Dynam.\ Systems \textbf{37} (2017), 1681--1764. \bibitem{Duderstadt1976} J.~J.~Duderstadt and L.~J.~Hamilton, \textit{Nuclear Reactor Analysis}, Wiley, New York, 1976. \bibitem{Grossi2026DNLS} P.~Nogueira Grossi, \textit{Differential Nonlinear Robustness of Critical States in Fibonacci and Tribonacci Substitution Chains}, G6 LLC (2026), Zenodo. Concept DOI: \href{https://doi.org/10.5281/zenodo.20026942}{10.5281/zenodo.20026942} (V4: \href{https://doi.org/10.5281/zenodo.20075822}{10.5281/zenodo.20075822}). [Includes \texttt{dnls\_nbonacci.py}, \texttt{TribonacciDNLS.lean}, and five figures; Lean~4/Mathlib4 verification of key analytic lemmas available at \texttt{github.com/TOTOGT/AXLE}.] \bibitem{Kohmoto1983} M.~Kohmoto, L.~P.~Kadanoff, and C.~Tang, \textit{Localization problem in one dimension: Mapping and escape}, Phys.\ Rev.\ Lett.\ \textbf{50} (1983), 1870--1872. \bibitem{Ostlund1983} S.~Ostlund, R.~Pandit, D.~Rand, H.~J.~Schellnhuber, and E.~D.~Siggia, \textit{One-dimensional Schr\"{o}dinger equation with an almost periodic potential}, Phys.\ Rev.\ Lett.\ \textbf{50} (1983), 1873--1876. \bibitem{Luck1989} J.~M.~Luck, \textit{Cantor spectra and scaling of gap widths in deterministic aperiodic systems}, Phys.\ Rev.\ B \textbf{39} (1989), 5834--5849. \bibitem{Mace2016} N.~Mac\'{e}, A.~Jagannathan, and F.~Pi\'{e}chon, \textit{Fractal dimensions of wave functions and local spectral measures on the Fibonacci chain}, Phys.\ Rev.\ B \textbf{93} (2016), 205134. \bibitem{Damanik2009} D.~Damanik and A.~Gorodetski, \textit{Hyperbolicity of the trace map for the weakly coupled Fibonacci Hamiltonian}, Nonlinearity \textbf{22} (2009), 123--143. \bibitem{Rauzy1982} G.~Rauzy, \textit{Nombres alg\'{e}briques et substitutions}, Bull.\ Soc.\ Math.\ France \textbf{110} (1982), 147--178. \bibitem{Flach2010} S.~Flach, M.~V.~Ivanchenko, and O.~I.~Kanakov, \textit{Spreading of wave packets in one-dimensional disordered lattices}, Phys.\ Rev.\ E \textbf{82} (2010), 036219. \bibitem{GrossiRepo2026} P.~Nogueira Grossi, \textit{Atratores: numerical code repository}, GitHub (2026). \href{https://github.com/grossi-ops/Atratores/tree/main/DNLS}{github.com/grossi-ops/Atratores}. \bibitem{Krebbekx2023} J.~P.~J.~Krebbekx, A.~Moustaj, K.~Dajani, and C.~de~Morais~Smith, \textit{Multifractal properties of tribonacci chains}, Phys.\ Rev.\ B \textbf{108} (2023), 104204. \href{https://doi.org/10.1103/PhysRevB.108.104204}{doi:10.1103/PhysRevB.108.104204}. \bibitem{Allaire1999} G.~Allaire and G.~Bal, \textit{Homogenization of the criticality spectral equation in neutron transport}, ESAIM: Math.\ Model.\ Numer.\ Anal.\ \textbf{33} (1999), 721--746. \bibitem{Varma2019} V.~K.~Varma, C.~de~Tomasi, and M.~M\"{u}ller, \textit{Diffusion in quasiperiodic Fibonacci chains}, Phys.\ Rev.\ B \textbf{100} (2019), 085105. \bibitem{Shechtman1984} D.~Shechtman, I.~Blech, D.~Gratias, and J.~W.~Cahn, \textit{Metallic phase with long-range orientational order and no translational symmetry}, Phys.\ Rev.\ Lett.\ \textbf{53} (1984), 1951--1953. \bibitem{Squires1978} G.~L.~Squires, \textit{Introduction to the Theory of Thermal Neutron Scattering}, Cambridge University Press, Cambridge, 1978. \bibitem{Stacey2007} W.~M.~Stacey, \textit{Nuclear Reactor Physics}, 2nd ed., Wiley-VCH, Weinheim, 2007. \end{thebibliography} % ====================================================================== % Figures (to be placed in §3 and §4 during final assembly) % ====================================================================== \clearpage \begin{figure}[h] \centering \includegraphics[width=\textwidth]{fig1_chain_structure} \caption{% Site-by-site symbol sequence for the Fibonacci ($n=2$, top) and Tribonacci ($n=3$, bottom) substitution chains, shown for the first 80 sites at generation $g=7$ and $g=6$ respectively. Blue (A) sites carry fissile parameters ($D=1.0$, $\Sigmar=0.5$, $\nuSigmaf=\lambda$); red (B) and green (C) sites are absorbers ($D=1.0$, $\Sigmar=2.0$, $\nuSigmaf=0$). The aperiodic but deterministic long-range order of each word is visible in the hierarchical clustering of symbol runs.} \label{fig:chain} \end{figure} \begin{figure}[h] \centering \includegraphics[width=\textwidth]{fig2_keff_vs_lambda} \caption{% Effective multiplication factor $\keff$ as a function of fission strength $\lambda$ for Fibonacci (left) and Tribonacci (right) chains at substitution generations $g = 4, 6, 8$. The dashed line marks the criticality condition $\keff = 1$. Fibonacci chains converge smoothly to criticality with increasing generation; the Tribonacci chain requires larger $\lambda$ to reach $\keff = 1$ at the same generation, reflecting the higher absorber density imposed by the three-symbol alphabet.} \label{fig:keff-lambda} \end{figure} \begin{figure}[h] \centering \includegraphics[width=\textwidth]{fig3_lambda_c_gap} \caption{% Left: critical fission strength $\lc(n)$ versus transfer-matrix spectral gap $\Delta_n = \rhon - |\rho_n^{(2)}|$ for $n = 2,3,4,5$. The dashed line is the linear fit $\lc \approx 0.958\,\Delta_n + 0.107$ ($r = 0.989$). Right: $\lc(n)$, $\Delta_n$, and the $n$-bonacci constant $\rhon$ plotted against $n$. The dash-dotted line marks $7/6$, to which $\lc(n)$ converges exactly for $n \geq 4$.} \label{fig:lambda-c-gap} \end{figure} \begin{figure}[h] \centering \includegraphics[width=\textwidth]{fig4_convergence} \caption{% Left: $\keff$ versus chain length $N$ at fixed $\lambda = \lc^{\rm Fib} \approx 1.064$, showing slower convergence for Tribonacci than Fibonacci. Right: absolute deviation $|\lc(g) - \lc(\infty)|$ versus generation $g$ on a semi-logarithmic scale, confirming exponential convergence consistent with the substitution spectral gap.} \label{fig:convergence} \end{figure} \begin{figure}[h] \centering \includegraphics[width=\textwidth]{fig5_flux_profiles} \caption{% Fundamental flux modes $\phi(x)$ (normalised to unit peak) for Fibonacci (top row) and Tribonacci (bottom row) chains at generations $g=5$ and $g=8$, evaluated at the respective critical fission strength $\lc(n)$. Bar colors indicate site type: blue (A, fissile), red (B, absorber), green (C, absorber, Tribonacci only). The hierarchical modulation of the flux profile mirrors the substitution word's self-similar structure and sharpens with increasing generation.} \label{fig:flux} \end{figure} \end{document}