\documentclass[pdflatex,sn-mathphys-num]{sn-jnl}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{lmodern,microtype}
\usepackage{amsmath,amssymb,mathtools,bm}
\usepackage{graphicx,booktabs,tabularx,array}
\usepackage{enumitem}
\usepackage{placeins,xurl,flafter,needspace}
\geometry{reset,a4paper,margin=25mm,bindingoffset=0mm,
          headheight=5.5pt,headsep=5.6mm,footskip=10mm}

% Typography follows the author's preferred 11-point article layout.
\makeatletter
\renewcommand\Titlefont{\normalfont\fontsize{16bp}{21bp}\selectfont\centering}
\renewcommand\normalsize{%
  \@setfontsize\normalsize{10.95pt}{13.6pt}%
  \abovedisplayskip=11pt plus 3pt minus 6pt
  \abovedisplayshortskip=0pt plus 3pt
  \belowdisplayshortskip=6.5pt plus 3.5pt minus 3pt
  \belowdisplayskip=\abovedisplayskip
  \let\@listi\@listI}
\renewcommand\small{%
  \@setfontsize\small{10pt}{12pt}%
  \abovedisplayskip=10pt plus 2pt minus 5pt
  \abovedisplayshortskip=0pt plus 3pt
  \belowdisplayshortskip=6pt plus 3pt minus 3pt
  \belowdisplayskip=\abovedisplayskip}
\renewcommand\abstractfont{\normalfont\small\unboldmath
  \leftskip=24pt\rightskip=24pt\parfillskip=0pt plus 1fil}
\renewcommand\keywordfont{\normalfont\small
  \leftskip=24pt\rightskip=24pt plus .5fil}
\renewcommand\figurecaptionfont{\normalfont\normalsize}
\renewcommand\tablecaptionfont{\normalfont\normalsize}
\AtBeginDocument{%
  \renewenvironment{proof}[1][\proofname]{%
    \par\pushQED{\qed}\normalfont\normalsize
    \topsep=6pt plus 6pt\relax\trivlist
    \item[\hskip\labelsep\itshape #1\@addpunct{.}]\ignorespaces
  }{\popQED\endtrivlist\@endpefalse}%
  \normalsize
  \setlength{\parindent}{1.2em}%
  \setlength{\parskip}{0pt}%
  \allowdisplaybreaks[1]\raggedbottom}
\makeatother
\normalsize
\emergencystretch=1.5em
\setlist{itemsep=2pt,topsep=5pt}
\graphicspath{{figures/}}

\theoremstyle{definition}
\newtheorem{definition}{Definition}[section]
\newtheorem{assumption}[definition]{Assumption}
\theoremstyle{plain}
\newtheorem{theorem}[definition]{Theorem}
\newtheorem{proposition}[definition]{Proposition}
\newtheorem{lemma}[definition]{Lemma}
\newtheorem{corollary}[definition]{Corollary}
\theoremstyle{remark}
\newtheorem{remark}[definition]{Remark}
\newtheorem{example}[definition]{Example}

\newcommand{\Z}{\mathbb Z}
\newcommand{\R}{\mathbb R}
\newcommand{\N}{\mathbb N}
\newcommand{\C}{\mathbb C}
\newcommand{\cA}{\mathcal A}
\newcommand{\cD}{\mathcal D}
\newcommand{\cH}{\mathcal H}
\newcommand{\cN}{\mathcal N}
\newcommand{\dd}{\mathrm d}
\newcommand{\e}{\mathrm e}
\newcommand{\ket}[1]{\lvert#1\rangle}
\newcommand{\norm}[1]{\lVert#1\rVert}
\newcommand{\abs}[1]{\lvert#1\rvert}
\newcommand{\round}{\mathcal R}
\newcommand{\cen}{\operatorname{cen}_N}
\newcommand{\cn}{\operatorname{cn}}
\newcommand{\sn}{\operatorname{sn}}
\newcommand{\dn}{\operatorname{dn}}
\newcommand{\supp}{\operatorname{supp}}
\newcommand{\diag}{\operatorname{diag}}
\newcommand{\tr}{\operatorname{tr}}
\hypersetup{pdftitle={Nonlinear Wave Limits of Reversible Finite-State Cellular Automata},pdfauthor={Zhenyu Sun},pdfsubject={Reversible cellular automata and semilinear wave equations}}

\begin{document}
\title[Nonlinear Wave Limits of Reversible Finite-State Cellular Automata]{Nonlinear Wave Limits of Reversible Finite-State Cellular Automata}
\author*[1]{\fnm{Zhenyu} \sur{Sun}}\email{sunzhenyu@neuq.edu.cn}
\affil*[1]{\orgdiv{School of Mathematics and Statistics}, \orgname{Northeastern University at Qinhuangdao}}
\abstract{We construct reversible finite-state cellular automata whose integer lifts approximate smooth solutions of one-dimensional semilinear wave equations. The kick--drift update is bijective for every odd modulus. Separate rounding of spatial edge contributions and local nonlinear forces gives an exact one-cell support bound for the full local operator algebra. Under a strict Courant condition, a discrete energy estimate bounds the field and staggered-velocity errors by $C_T(h^2+\tau^2+\delta/\tau^2)$ on a fixed physical time interval, where $h$, $\tau$, and $\delta$ are the spatial, temporal, and field spacings. An explicit modulus bound prevents modular wrapping throughout that interval. The choices $\tau=\nu h$, $0<\nu<1$, and $\delta=O(h^4)$ yield second-order convergence. Coupled scalar fields have a common-speed, locally Lorentz-covariant classical limit. Equal bare masses and equal quartic self-couplings preserve a symmetric invariant manifold. Unequal quartic couplings introduce a cubic source into the antisymmetric field equation, and we derive the first variation of the antisymmetric field. Numerical experiments compare exact integer trajectories with a nonlinear traveling wave and an independent continuum solver, measure light and massive dispersion branches, and examine continuously interacting wave packets. The accompanying programs reproduce the reported calculations.}
\keywords{Reversible cellular automata, Finite-state dynamics, Semilinear wave equation, Integer integration, Continuum limit, Causality}
\maketitle

\section{Introduction}
\label{sec:introduction}

Reversible cellular automata provide finite-state models with local, invertible updates. Approximation of a nonlinear field equation requires control of accumulated discretization and rounding errors as the spatial mesh is refined. Locality bounds the support of evolved observables, while convergence relates the values of the discrete fields to a specified continuum solution. We establish both properties for a family of rounded integer wave equations and exhibit an exact time-reversal involution.

The relation between reversible automata and field dynamics has been studied in the cellular automaton interpretation of quantum mechanics~\cite{tHooft2016}. A direct scalar-field antecedent is 't Hooft's integer automaton in one space and one time dimension~\cite{tHooft2013Duality}. Its massless recurrence is the centered wave equation at Courant number one, with explicit left- and right-moving solutions. Adding masses and nonlinear forces requires stability and resolution estimates and a finite register bound for the integer trajectory.

The reversible update used here belongs to the established memory construction for cellular automata~\cite{ToffoliMargolus1990,Kari2018}. Related rounded second-order recurrences appear in the exactly reversible molecular dynamics of Levesque and Verlet~\cite{LevesqueVerlet1993} and in integer integration methods such as JANUS~\cite{ReinTamayo2018}. For a spatially discretized wave equation, the Lipschitz constant of the full force grows as $h^{-2}$, where $h$ is the mesh spacing. Directly iterating this Lipschitz bound does not give the uniform stability estimate needed for the continuum limit. Discrete wave energy and a strict Courant condition instead give stability estimates uniform in $h$~\cite{Strikwerda2004,HairerLubichWanner2006}.

Locality also requires an operator-level argument. A reversible classical rule acts by pullback on diagonal observables; its permutation representation acts on the complete local matrix algebra. Off-diagonal operators can acquire a larger neighborhood than classical observables~\cite{ArrighiNesmeWerner2017}. Quantum cellular automata provide the corresponding algebraic setting~\cite{SchumacherWerner2004,Arrighi2011,Farrelly2020}. We decompose the spatial force into edge contributions and keep the nonlinear force local to each site. The resulting commuting kick operators give an exact one-cell support bound for the complete local algebra in both time directions.

The continuum equation is
\begin{equation}
 \partial_t^2 U-\partial_x^2 U+f(U)=0,
 \qquad U:(\mathbb R/\ell\mathbb Z)\times[-\tau_0,T]\longrightarrow\mathbb R^s,
 \label{eq:intro-pde}
\end{equation}
where $f$ is $C^1$ and $U$ is smooth on the stated interval. For a gradient force $f=\nabla V$, the common wave speed and scalar potential yield a locally Lorentz-covariant classical field equation. Periodic boundary conditions give a finite-volume problem; Lorentz boosts act on the covering spacetime.

Theorem~\ref{thm:finite-alphabet-wave-convergence} quantifies the approximation by finite-alphabet automata. Let $W^n$ be the physical integer field, let $U_h^n=(U(jh,n\tau))_j$, and write $d_t^-z^n=(z^n-z^{n-1})/\tau$. With either second-order initialization specified in that theorem,
\begin{equation}
 \max_{0\le n\tau\le T}
 \left(\|W^n-U_h^n\|_{1,h}
       +\|d_t^-W^n-d_t^-U_h^n\|_h\right)
 \le C_T\left(h^2+\tau^2+
                      \frac{\delta}{\tau^2}\right).
 \label{eq:intro-bound}
\end{equation}
Here $\delta$ is the field spacing. The discrete norms are defined in Section~\ref{subsec:discrete-wave-energy}. Under a strict Courant condition, $C_T$ is independent of the mesh and modulus. A one-dimensional Sobolev estimate controls the local nonlinearity. An explicit register bound then ensures agreement of the modular and integer trajectories. Second-order initialization, $\tau=\nu h$ with $0<\nu<1$, and $\delta=O(h^4)$ give second-order convergence. The local alphabet grows with the field resolution.

For two coupled scalar fields with equal bare masses and equal quartic self-couplings, exchange symmetry preserves equal field and memory registers exactly. In the zero-bare-mass case, the symmetric continuum field obeys a scalar quartic wave equation, while the antisymmetric branch has a mass determined by the quadratic coupling. Unequal quartic self-couplings introduce a cubic source into the antisymmetric field equation. We derive the first variation of the antisymmetric field and compare it with independent continuum solutions. For continuously interacting wave packets, matched single-packet solutions distinguish single-packet generation from the failure of nonlinear superposition.

The numerical experiments test the same integer update against an exact nonlinear traveling wave and an independent Fourier discretization of the coupled continuum equations. They also measure dispersion under refinement, nonlinear field generation, energy variation, and high-frequency power. The accompanying programs generate complete integer trajectories and reproduce the reported results.

Section~\ref{sec:reversibleca} establishes reversibility and operator support. Section~\ref{sec:quantized-model} defines the rounded force and its integer lift. Section~\ref{sec:controlled-wave-limit} proves the approximation theorem and constructs the exact traveling wave. Sections~\ref{sec:coupled} and~\ref{sec:numerics} treat the coupled model and numerical experiments. Section~\ref{sec:discussion} discusses the continuum scaling and its relation to quantum observables.

\Needspace{8\baselineskip}
\section{Finite-state reversibility and causal propagation}
\label{sec:reversibleca}

\subsection{Classical observables and the full operator algebra}

Let $\Lambda$ be a finite graph, with graph distance $d$, and let
$\Sigma$ be a finite alphabet.  The periodic lattices used below are
special cases.  For $R\subseteq\Lambda$, write
$R^{(a)}=\{x:d(x,R)\leq a\}$ and
$\mathcal N_a(x)=\{x\}^{(a)}$.  The configuration space and its Hilbert
space representation are
\begin{equation*}
 \Omega=\Sigma^{\Lambda},\qquad
 \cH=\ell^2(\Omega)=\bigotimes_{x\in\Lambda}\ell^2(\Sigma).
\end{equation*}
A bijection $\mathsf F:\Omega\to\Omega$ defines the permutation unitary
$U_{\mathsf F}\ket{\omega}=\ket{\mathsf F(\omega)}$.  We distinguish two local algebras:
\begin{equation*}
 \cA_R=\mathcal B\!\left(\bigotimes_{x\in R}\ell^2(\Sigma)\right)
       \otimes I_{R^c},
 \qquad
 \mathcal D_R=\{M_f:f:\Omega\to\mathbb C
                  \text{ depends only on }\omega|_R\},
\end{equation*}
where $M_f\ket{\omega}=f(\omega)\ket{\omega}$.  Thus
$\mathcal D_R$ is the algebra of classical, diagonal observables and
$\cA_R$ also contains operators that change basis configurations.
The full operator algebra may have a larger dependence neighborhood than its diagonal subalgebra~\cite{ArrighiNesmeWerner2017}.
A map has classical dependence radius at most $r$ if its output at each site $x$ is determined by the input on $\mathcal N_r(x)$.

\begin{proposition}[Diagonal causal cone]
\label{prop:diagonalcone}
Suppose $\mathsf F$ and $\mathsf F^{-1}$ have classical dependence radii $r_+$ and
$r_-$, respectively.  For every integer $n\geq0$,
\begin{align}
 U_{\mathsf F}^{-n}\mathcal D_R U_{\mathsf F}^n
   &\subseteq\mathcal D_{R^{(nr_+)}},\label{eq:diagonalforward}\\
 U_{\mathsf F}^{n}\mathcal D_R U_{\mathsf F}^{-n}
   &\subseteq\mathcal D_{R^{(nr_-)}}.\label{eq:diagonalbackward}
\end{align}
Consequently, if $A\in\mathcal D_R$, $B\in\cA_S$, and
$d(R,S)>nr_+$, then $[U_{\mathsf F}^{-n}AU_{\mathsf F}^n,B]=0$.
\end{proposition}

\begin{proof}
On a basis vector,
$U_{\mathsf F}^{-1}M_fU_{\mathsf F}\ket{\omega}=f(\mathsf F(\omega))\ket{\omega}$.
The pullback $f\circ \mathsf F$ depends only on $R^{(r_+)}$.  Iteration proves
\eqref{eq:diagonalforward}; the same argument applied to $\mathsf F^{-1}$
proves \eqref{eq:diagonalbackward}.  Operators on disjoint tensor
factors commute, which gives the last statement.
\end{proof}

If $B$ is also diagonal, this commutator vanishes at every separation,
including inside the cone.  A nontrivial support test therefore uses
classical dependence directly or an off-diagonal probe.  The next example exhibits the additional support of an off-diagonal operator.

\begin{proposition}[Different classical and operator radii]
\label{prop:toffolicounterexample}
There is an involutive classical automaton of radius one whose
permutation unitary takes a one-site operator to an operator with
support at distance two in one step.
\end{proposition}

\begin{proof}
Take a periodic chain of length $J\geq7$, with two bits $(a_x,b_x)$
per site, and arithmetic in $\mathbb Z_2$.  Define
\begin{equation*}
 a'_x=a_x,\qquad b'_x=b_x\oplus(a_{x-1}a_{x+1}).
\end{equation*}
The $a$ registers do not change, so $\mathsf F^2=I$, and both $\mathsf F$ and $\mathsf F^{-1}$
have radius one.  Let $X_{a_x}$ flip bit $a_x$.  Computing the
conjugation on basis states gives
\begin{equation}
 U_{\mathsf F}^\dagger X_{a_0}U_{\mathsf F}
 =X_{a_0}\,
   \operatorname{CNOT}_{a_{-2}\to b_{-1}}\,
   \operatorname{CNOT}_{a_2\to b_1},
 \label{eq:toffolispread}
\end{equation}
where a controlled-not adds its control bit to its target bit.  Indeed,
flipping $a_0$ between the two applications of $\mathsf F$ changes the second
update of $b_{-1}$ by $a_{-2}$ and that of $b_1$ by $a_2$; all other
updates cancel.  Denote the operator in \eqref{eq:toffolispread} by
$C$.  Applied to the all-zero state, $CX_{a_2}$ produces nonzero bits
$a_0,a_2,b_1$, whereas $X_{a_2}C$ produces only $a_0,a_2$.
Thus $[C,X_{a_2}]\neq0$, proving that $C$ is not supported in
$\{0\}^{(1)}$.
\end{proof}

This example illustrates the quantum-neighbourhood distinction studied in
Ref.~\cite{ArrighiNesmeWerner2017}. For the automata constructed here,
the sharper full-algebra bound follows from the local gate structure.

\subsection{Kick--drift dynamics and its inverse}

Let $N\geq2$, $s\geq1$, and $G=(\mathbb Z_N)^s$, so that every site has a field register
$Q_x\in G$ and an increment register $P_x\in G$.
For arbitrary functions $\mathcal G_x$ depending only on
$Q|_{\mathcal N_r(x)}$, define the two shears
\begin{align}
 K:(Q,P)&\longmapsto(Q,P+\mathcal G(Q)),\nonumber\\
 D:(Q,P)&\longmapsto(Q+P,P),\qquad T=D\circ K.
 \label{eq:structuralkickdrift}
\end{align}
All arithmetic in this subsection is in $G$. The functions may be specified by rounding or arbitrary finite tables.

\begin{proposition}[Reversibility and time reversal]
\label{prop:structuralinverse}
The map $T$ is a bijection, with inverse
\begin{equation}
 Q=Q'-P',\qquad P=P'-\mathcal G(Q).
 \label{eq:structuralinverse}
\end{equation}
Both $T$ and $T^{-1}$ have classical radius at most $r$.  Moreover,
\begin{equation}
 \mathcal J(Q,P)=(Q-P,-P),\qquad
 \mathcal J^2=I,\qquad \mathcal JT\mathcal J=T^{-1}.
 \label{eq:structuralreversal}
\end{equation}
The finite-state dynamics preserves counting measure and is
represented by a permutation unitary whose eigenvalues are roots of
unity.
\end{proposition}

\begin{proof}
The drift equation gives $Q=Q'-P'$ on each site.  Evaluating
$\mathcal G$ on this recovered field gives the second equation in
\eqref{eq:structuralinverse}.  These equations define a two-sided inverse
and use at most the stated neighbourhood.  To establish the reversal
identity, change onsite coordinates to
$(q,r)=(Q,Q-P)$.  In these coordinates,
\begin{equation}
 \widehat T(q,r)=(2q-r+\mathcal G(q),q).
 \label{eq:memoryrecurrence}
\end{equation}
The swap $S(q,r)=(r,q)$ satisfies
$S\widehat T S=\widehat T^{-1}$, as direct substitution shows.
Transforming $S$ back to $(Q,P)$ gives $\mathcal J$ in
\eqref{eq:structuralreversal}.  A bijection of a finite set preserves
counting measure.  Its cycles give permutation blocks with eigenvalues
$\exp(2\pi i k/m)$ for cycle lengths $m$, proving the final claim.
\end{proof}

The increment $P$ is staggered in time. Along a lifted trajectory, $P^n=Q^n-Q^{n-1}$, and time reversal exchanges the two stored field levels. This gives the involution $\mathcal J$ in~\eqref{eq:structuralreversal}.

The memory construction is standard in reversible cellular automata
\cite{ToffoliMargolus1987,ToffoliMargolus1990}.  Equation
\eqref{eq:memoryrecurrence} also makes its connection to integer
integration explicit.  The rounded second-order recurrence of
Levesque and Verlet is of this form \cite{LevesqueVerlet1993}; related
exact integer lattice maps were studied by Earn and Tremaine
\cite{EarnTremaine1992}.  JANUS develops bitwise reversible integer
integration further, including higher-order compositions for
$N$-body dynamics \cite{ReinTamayo2018}.  We use this established reversible construction to study spatial locality and the approximation of nonlinear waves.

\subsection{Full operator cones from commuting kicks}

\begin{theorem}[Full-algebra locality for finite-range kicks]
\label{thm:generaloperatorcone}
For the dynamics \eqref{eq:structuralkickdrift}, the permutation unitary
$U_T$ obeys, for every $n\geq0$,
\begin{equation}
 U_T^{\mp n}\cA_R U_T^{\pm n}\subseteq\cA_{R^{(2nr)}}.
 \label{eq:generaloperatorcone}
\end{equation}
The bound holds for arbitrary functions $\mathcal G_x$ of range $r$.
\end{theorem}

\begin{proof}
Write $K=\prod_x K_x$, where $K_x$ translates only $P_x$ by
$\mathcal G_x(Q)$ and reads only the $Q$ registers in
$\mathcal N_r(x)$.  These gates commute, including when their
supports overlap, because every gate leaves all $Q$ registers fixed
and translations of the $P$ registers commute.  The same is true of
their permutation unitaries.  Also $D$ is a product of onsite gates.

For $A\in\cA_R$, the operator $B=U_D^\dagger A U_D$ remains in
$\cA_R$.  In $U_K^\dagger B U_K$, reorder the mutually commuting
kick gates so that all gates with
$\mathcal N_r(x)\cap R=\varnothing$ cancel across $B$.  Only centres
$x\in R^{(r)}$ remain, and the union of their supports lies in
$R^{(2r)}$.  This proves the one-step bound for
$U_T^\dagger A U_T=U_K^\dagger U_D^\dagger A U_D U_K$.
For $U_T A U_T^\dagger$, first use the same cancellation argument for
$U_K A U_K^\dagger$, then apply the onsite drift.  Iterating the two
one-step bounds proves \eqref{eq:generaloperatorcone}.
\end{proof}

The factor two comes from the diameter bound on the support of a general range-$r$ kick.
It can be improved when the force decomposes over smaller interaction
sets.  In particular, a pairwise decomposition gives the radius-one
bound needed for the models in this paper.

\begin{theorem}[Radius-one full-algebra cone for pairwise forces]
\label{thm:pairwisecone}
Suppose the force has the form
\begin{equation}
 \mathcal G_x(Q)=H_x(Q_x)+\sum_{y\sim x}H_{xy}(Q_x,Q_y),
 \label{eq:pairwiseforce}
\end{equation}
with arbitrary functions valued in $G$.  Then both $U_T$ and its
inverse have full operator radius at most one:
\begin{equation}
 U_T^{\mp n}\cA_R U_T^{\pm n}\subseteq\cA_{R^{(n)}}
 \quad(n\geq0).
 \label{eq:pairwiseoperatorcone}
\end{equation}
In particular, for $A\in\cA_R$, $B\in\cA_S$ and $d(R,S)>n$,
\begin{equation}
 [U_T^{\mp n}AU_T^{\pm n},B]=0.
 \label{eq:pairwisecommutator}
\end{equation}
\end{theorem}

\begin{proof}
Split the kick into onsite translations
$P_x\mapsto P_x+H_x(Q_x)$ and directed-edge translations
$P_x\mapsto P_x+H_{xy}(Q_x,Q_y)$.  Each edge gate has support
$\{x,y\}$, and all these gates commute for the same reason as in
Theorem~\ref{thm:generaloperatorcone}.  In the conjugation of an
operator supported on $R$, all gates whose support is disjoint from
$R$ can therefore be reordered and cancelled.  The remaining gates
act onsite in $R$ or on edges incident to $R$, so their support is contained
in $R^{(1)}$.  The onsite drift does not enlarge support.  This proves
both one-step inclusions.  Iteration gives
\eqref{eq:pairwiseoperatorcone}, and the commutator follows from
disjoint tensor supports.
\end{proof}

For example, take an odd modulus $N$, let $\operatorname{cen}_N:G\to\mathbb Z^s$ select centered representatives, and let $\round$ denote componentwise nearest-integer rounding, with ties away from zero.
The force used below is the reduction modulo $N$ of
\begin{equation}
 \begin{split}
 \mathcal G_x(Q)\equiv{}&
  \sum_{y=x\pm1}\round\!\left(
     \nu^2[\operatorname{cen}_N(Q_y)-\operatorname{cen}_N(Q_x)]\right)\\
  &-\round\!\left(
     \frac{\tau^2}{\delta}
     f(\delta\operatorname{cen}_N(Q_x))\right)\pmod N.
 \end{split}
 \label{eq:pairwisequantizedforce}
\end{equation}
Here $f:\mathbb R^s\to\mathbb R^s$ may couple the species at the
same site.  Each spatial term is an edge function and the final term
is onsite.  Thus the full-algebra radius-one theorem applies for arbitrary $f$ at every
configuration, including configurations that wrap modulo $N$.
With lattice spacing $h$ and timestep
$\tau=\nu h$, the resulting microscopic support speed is $h/\tau$
in physical units.

These support bounds hold at fixed lattice spacing and modulus. The next section introduces the physical scales needed to compare integer trajectories with solutions of~\eqref{eq:intro-pde}.

\section{Integer lifts and rounded wave dynamics}
\label{sec:quantized-model}

\subsection{Representatives and force rounding}

Let $J\ge3$, $h=\ell/J$, and $\tau=\nu h$, where $0<\nu<1$. Spatial indices belong to $\Z/J\Z$. Let $N$ be odd and write
\begin{equation*}
 I_N=\{-(N-1)/2,\ldots,(N-1)/2\}.
\end{equation*}
The map $\operatorname{cen}_N:\Z_N\to I_N$ selects the unique centered representative. All differences inside a rounded force are ordinary integer differences of these representatives. This convention defines the force on every configuration, including configurations at the modular boundary.

The odd nearest-integer map is
\begin{equation}
 \round(z)=\begin{cases}
       \lfloor z+1/2\rfloor,&z\ge0,\\
       -\lfloor-z+1/2\rfloor,&z<0.
             \end{cases}
 \qquad
 \abs{\round(z)-z}\le\frac12.
 \label{eq:rounding}
\end{equation}
Thus ties are rounded away from zero and $\round(-z)=-\round(z)$. For rational inputs, rounding can be performed by signed integer division without floating-point comparisons.

At each site store $Q_j,P_j\in(\Z_N)^s$. For a field resolution $\delta>0$, define the integer force, component by component, by
\begin{equation}
 \begin{split}
 G_j^a(Q)={}&
 \round\!\left(\nu^2(\bar Q_{j+1}^a-\bar Q_j^a)\right)
 +\round\!\left(\nu^2(\bar Q_{j-1}^a-\bar Q_j^a)\right)\\
 &-\round\!\left(\frac{\tau^2}{\delta}
               f_a(\delta\bar Q_j)\right),
 \qquad \bar Q_j=\operatorname{cen}_N Q_j.
 \end{split}
 \label{eq:rounded-force}
\end{equation}
The finite-state update is
\begin{equation}
 P_j'=P_j+G_j(Q)\pmod N,
 \qquad Q_j'=Q_j+P_j'\pmod N.
 \label{eq:new-update}
\end{equation}
For polynomial $f$ with rational coefficients and rational $h,\tau,\delta$, the entire update can be evaluated with integer arithmetic. Its inverse recovers $Q=Q'-P'$ and then subtracts the force evaluated at the recovered fields. Each site has $N^{2s}$ states.

Each of the two spatial summands in \eqref{eq:rounded-force} depends on one edge, and the last summand depends on one site. The full-algebra support result of Section~\ref{sec:reversibleca} therefore applies with one-cell radius in both time directions. In physical units the microscopic support bound is $h/\tau=1/\nu$. The limiting wave speed in \eqref{eq:intro-pde} is one.

\subsection{The physical integer recurrence}

Starting from ordinary integers $Q^{-1},Q^0$, define an unrestricted integer recurrence using \eqref{eq:rounded-force} without the bars, and put
\begin{equation*}
 P^n=Q^n-Q^{n-1},\qquad W^n=\delta Q^n.
\end{equation*}
As long as $Q^n,P^n\in I_N$, reduction modulo $N$ intertwines this recurrence with \eqref{eq:new-update}.

\begin{lemma}[Quantization residual]
\label{lem:quantization-residual}
The unrestricted physical lift satisfies
\begin{equation}
 d_t^2W_j^n-\Delta_hW_j^n+f(W_j^n)=\rho_j^n,
 \qquad
 \abs{\rho_j^{a,n}}\le\frac{3\delta}{2\tau^2},
 \label{eq:physical-recurrence}
\end{equation}
where $d_t^2W^n=(W^{n+1}-2W^n+W^{n-1})/\tau^2$ and
$\Delta_hW_j=(W_{j+1}-2W_j+W_{j-1})/h^2$. Consequently,
\begin{equation*}
 \norm{\rho^n}_{2,h}\le\frac{3\sqrt{s\ell}}{2}\frac{\delta}{\tau^2},
 \qquad \norm{z}_{2,h}^2=h\sum_j\abs{z_j}^2.
\end{equation*}
\end{lemma}
\begin{proof}
Each rounded summand differs from its argument by at most $1/2$. Multiplication by $\delta$ and division by $\tau^2$ give the componentwise bound. The identity $\nu^2/\tau^2=h^{-2}$ identifies the spatial term. Summing over $J$ sites and $s$ components proves the norm estimate.
\end{proof}


\subsection{A priori and trajectory-based register bounds}

A sufficient no-wrap condition on a finite interval is
\begin{equation}
 \max_{n,j,a}\bigl(\abs{Q_j^{a,n}},\abs{P_j^{a,n}}\bigr)<\frac N2.
 \label{eq:nowrap-new}
\end{equation}
The convergence theorem below supplies an a priori bound. In the computations, \eqref{eq:nowrap-new} is also checked at every saved step, and an independently reduced modular trajectory is compared with the integer trajectory.

\subsection{The unrounded reference and its stability}

For a linear force $f(w)=\mu^2w$, the unrounded physical recurrence has the normal-mode relation
\begin{equation}
 \frac{4}{\tau^2}\sin^2\frac{\omega\tau}{2}
 =\mu^2+\frac{4}{h^2}\sin^2\frac{kh}{2}.
 \label{eq:stable-reference-dispersion}
\end{equation}
For a symmetric positive semidefinite mass matrix, replace $\mu^2$ by each of its eigenvalues. All modes with positive squared frequency are in the oscillatory range if
\begin{equation}
 4\nu^2+\tau^2\mu_{\max}^2<4.
 \label{eq:full-band-cfl}
\end{equation}
The zero-frequency mode is treated with its specified initial velocity. At the endpoints of the dimensionless stability interval $0\le\zeta\le4$, the recurrence $q_{n+1}-(2-\zeta)q_n+q_{n-1}=0$ has repeated roots and admits solutions with linear growth in $n$. The strict inequality in~\eqref{eq:full-band-cfl} excludes the upper endpoint.

Choose the principal frequency $0\leq\omega\tau<\pi$. For fixed physical $k,\mu$ and $\tau=\nu h$, expansion of \eqref{eq:stable-reference-dispersion} gives
\begin{equation}
 \omega^2=\mu^2+k^2+
 \frac{\tau^2}{12}(\mu^2+k^2)^2-
 \frac{h^2}{12}k^4+O(h^4).
 \label{eq:reference-expansion}
\end{equation}
The common limiting speed is thus obtained by refining both space and time at fixed physical mass. The rounding residual in Lemma~\ref{lem:quantization-residual} measures the difference between this reference recurrence and the integer update. The next section controls its accumulated effect over a fixed physical time interval.

\section{Controlled approximation of semilinear waves}
\label{sec:controlled-wave-limit}

The approximation depends on the spatial spacing $h$, the time step $\tau$, and the field spacing $\delta$. The modulus $N$ determines the available field range. We first estimate the error of the unrestricted integer recurrence and then choose $N$ large enough to prevent modular wrapping.

The update and its unrestricted integer lift are defined in
Section~\ref{sec:quantized-model}. The target is
\eqref{eq:intro-pde}; for use below, set
\begin{equation*}
 (D_hw)_j=(w_{j+1}-w_j)/h,\qquad
 (\Delta_hw)_j=(w_{j+1}-2w_j+w_{j-1})/h^2.
\end{equation*}
We use $\rho^n$ for the rounding residual in~\eqref{eq:physical-recurrence}. Its bound is uniform over all integer configurations. The norm $\|\cdot\|_h$ below is the norm $\|\cdot\|_{2,h}$ defined in Lemma~\ref{lem:quantization-residual}.

\subsection{A discrete energy estimate}
\label{subsec:discrete-wave-energy}

Use the inner product and norms
\begin{equation*}
 \langle z,w\rangle_h=h\sum_j z_j\cdot w_j,\qquad
 \|z\|_h^2=\langle z,z\rangle_h,\qquad
 \|z\|_{1,h}^2=\|z\|_h^2+\|D_hz\|_h^2.
\end{equation*}
Summation by parts gives
$\langle-\Delta_hz,w\rangle_h=\langle D_hz,D_hw\rangle_h$.
Let
\begin{equation*}
 A_h=I-\Delta_h,\qquad
 \kappa_h=1-\frac{\tau^2}{4}\left(1+\frac4{h^2}\right).
\end{equation*}
We impose the strict condition
\begin{equation}
 \kappa_h\geq\kappa_0>0.
 \label{eq:strict-energy-cfl}
\end{equation}
For fixed $\nu<1$, this holds on all sufficiently fine grids.
The experiments use $\nu=1/2$.

\begin{lemma}[Uniform forced-wave estimate]
\label{lem:uniform-wave-energy}
Suppose
\begin{equation*}
 \frac{e^{n+1}-2e^n+e^{n-1}}{\tau^2}+A_he^n=g^n.
\end{equation*}
For
\begin{equation*}
 v^{n+1/2}=\frac{e^{n+1}-e^n}{\tau},
 \qquad
 E^{n+1/2}
 =
 \frac12\|v^{n+1/2}\|_h^2
 +\frac12\langle A_he^{n+1},e^n\rangle_h,
\end{equation*}
condition \eqref{eq:strict-energy-cfl} implies
\begin{equation}
 \begin{split}
 E^{n+1/2}
 ={}&
 \frac12\left\langle
 \left(I-\frac{\tau^2}{4}A_h\right)v^{n+1/2},
 v^{n+1/2}\right\rangle_h\\
 &+\frac12\left\|
 A_h^{1/2}\frac{e^{n+1}+e^n}{2}\right\|_h^2,
 \end{split}
 \label{eq:positive-modified-energy}
\end{equation}
\begin{equation}
 \|v^{n+1/2}\|_h+\|e^n\|_{1,h}+\|e^{n+1}\|_{1,h}
 \leq C_{\kappa_0}\sqrt{E^{n+1/2}},
 \label{eq:energy-endpoint-control}
\end{equation}
\begin{equation}
 \left|\sqrt{E^{n+1/2}}-\sqrt{E^{n-1/2}}\right|
 \leq \frac{\tau}{\sqrt{2\kappa_0}}\|g^n\|_h.
 \label{eq:sqrt-energy-bound}
\end{equation}
The constants are independent of $h$, $\tau$, and the number of sites.
\end{lemma}

\begin{proof}
The eigenvalues of $-\Delta_h$ belong to $[0,4/h^2]$.
Writing $e^{n+1}$ and $e^n$ as their average plus or minus
$\tau v^{n+1/2}/2$ proves \eqref{eq:positive-modified-energy}.
Its first term is bounded below by $\kappa_0\|v^{n+1/2}\|_h^2/2$.
For either endpoint,
\begin{equation*}
 \|e^n\|_{1,h}
 \leq
 \left\|A_h^{1/2}\frac{e^{n+1}+e^n}{2}\right\|_h
 +\frac{\tau}{2}\|A_h^{1/2}v^{n+1/2}\|_h.
\end{equation*}
Since $\tau^2\|A_h\|/4\leq1-\kappa_0$, this gives
\eqref{eq:energy-endpoint-control}; the same calculation applies to
$e^{n+1}$.

Take the inner product of the recurrence with
$(e^{n+1}-e^{n-1})/(2\tau)$.  Symmetry of $A_h$ gives the exact identity
\begin{equation*}
 E^{n+1/2}-E^{n-1/2}
 =
 \tau\left\langle g^n,
 \frac{v^{n+1/2}+v^{n-1/2}}2\right\rangle_h.
\end{equation*}
The velocity bound yields
\begin{equation*}
 |E^{n+1/2}-E^{n-1/2}|
 \leq
 \frac{\tau}{\sqrt{2\kappa_0}}\|g^n\|_h
 \left(\sqrt{E^{n+1/2}}+\sqrt{E^{n-1/2}}\right).
\end{equation*}
Cancel the sum of the square roots.  If it vanishes, the claimed inequality
holds directly.
\end{proof}

The identity term in $A_h$ controls the spatially constant mode. In the error equation below, the corresponding $e^n$ term appears on the right-hand side, leaving the original recurrence unchanged.

\subsection{Finite-time convergence and a sufficient modulus}
\label{subsec:finite-time-convergence}

Set $\mathbb T_\ell=\mathbb R/\ell\mathbb Z$ and $x_j=jh$. Write $U_h^n=(U(x_j,n\tau))_j$, and let
$u_h^0=(U(x_j,0))_j$ and $v_h^0=(U_t(x_j,0))_j$.
The following initialization needs only the initial field and velocity:
\begin{align}
 Q^0&=\round(u_h^0/\delta),\qquad W^0=\delta Q^0,
 \nonumber\\
 Q^{-1}
 &=
 \round\!\left(
 \frac{W^0-\tau v_h^0+
 \frac{\tau^2}{2}(\Delta_hW^0-f(W^0))}{\delta}
 \right),
 \qquad P^0=Q^0-Q^{-1}.
 \label{eq:second-order-integer-initialization}
\end{align}
When the smooth initial acceleration is available, the alternative
\nopagebreak
\begin{align}
 Q^0&=\round(u_h^0/\delta),\qquad P^0=Q^0-Q^{-1},\nonumber\\*
 Q^{-1}&=\round\!\left(
 \frac{u_h^0-\tau v_h^0+\tfrac12\tau^2(U_{xx,h}^0-f(u_h^0))}{\delta}
 \right)
 \label{eq:analytic-taylor-initialization}
\end{align}
has the same order of accuracy, where $U_{xx,h}^0=(U_{xx}(x_j,0))_j$. The experiments use this formula with analytic or spectrally evaluated spatial derivatives.
An initialization using exact samples at $-\tau$ also satisfies the
estimates below. The scaled memory variable $\delta P^0/\tau$ approximates the velocity at $-\tau/2$. Even for zero velocity at time zero, the initial acceleration contributes to this staggered value.

\begin{theorem}[A finite-alphabet approximation theorem]
\label{thm:finite-alphabet-wave-convergence}
Fix $\ell>0$, $T>0$, $\tau_0>0$, $s\in\mathbb N$, and $0<\nu<1$.
Let $f:\mathbb R^s\to\mathbb R^s$ be $C^1$, and let
$U\in C^5(\mathbb T_\ell\times[-\tau_0,T];\mathbb R^s)$
solve \eqref{eq:intro-pde}.
Set
\begin{equation*}
 M=\sup_{x,t}|U(x,t)|,\qquad
 \varepsilon_h=h^2+\tau^2+\frac{\delta}{\tau^2},
 \qquad \tau=\nu h.
\end{equation*}
Assume $\tau\leq\min(\tau_0,1)$ and
\eqref{eq:strict-energy-cfl}.
There are constants $C_T$ and $\varepsilon_*>0$, depending on the
stated solution, $f,\ell,s,\nu,\kappa_0,T$, but independent of
$h,\delta,N$, with the following property.
If $\varepsilon_h\leq\varepsilon_*$, initialize by
\eqref{eq:second-order-integer-initialization} or
\eqref{eq:analytic-taylor-initialization} and choose an odd integer
\begin{equation}
 N>\frac{4(M+1)}{\delta}.
 \label{eq:sufficient-modulus}
\end{equation}
Put $n_T=\lfloor T/\tau\rfloor$. Then the finite-state trajectory has
an integer lift without wrap at stored field levels $-1\leq n\leq n_T$
and memory levels $0\leq n\leq n_T$, with transitions performed for
$0\leq n<n_T$, and
\begin{equation}
 \max_{0\leq n\leq\lfloor T/\tau\rfloor}
 \left(
 \|W^n-U_h^n\|_{1,h}
 +
 \left\|\frac{\delta P^n}{\tau}
 -
 \bigl(U_t(x_j,(n-\tfrac12)\tau)\bigr)_j\right\|_h
 \right)
 \leq C_T\varepsilon_h.
 \label{eq:main-convergence-estimate}
\end{equation}
In particular, $\delta=h^4/64$ and $\tau=h/2$ give a second-order
bound $C_T h^2$.
\end{theorem}

\begin{proof}
\emph{Globally Lipschitz forces.}
First suppose $f$ is globally Lipschitz, with constant $L_f$.
Taylor expansion of the smooth solution gives
\begin{equation*}
 \frac{U_h^{n+1}-2U_h^n+U_h^{n-1}}{\tau^2}
 =
 \Delta_hU_h^n-f(U_h^n)+\xi^n,
 \qquad
 \|\xi^n\|_h\leq C(\tau^2+h^2).
\end{equation*}
The Taylor remainders are uniform on the stated compact spacetime interval.
Subtract this identity from \eqref{eq:physical-recurrence}.
For $e^n=W^n-U_h^n$, we obtain
\begin{equation*}
 \frac{e^{n+1}-2e^n+e^{n-1}}{\tau^2}+A_he^n
 =
 e^n-\bigl(f(W^n)-f(U_h^n)\bigr)+\rho^n-\xi^n.
\end{equation*}
The deterministic rounding residual satisfies
\begin{equation*}
 \|\rho^n\|_h\leq
 \frac{3\sqrt{s\ell}}2\,\frac{\delta}{\tau^2}.
\end{equation*}
Let $\mathcal E_n=\sqrt{E^{n-1/2}}$, with $E$ defined for this error.
The endpoint estimate in Lemma~\ref{lem:uniform-wave-energy} gives
$\|e^n\|_h\leq C_{\kappa_0}\mathcal E_n$, and hence
\begin{equation*}
 \mathcal E_{n+1}
 \leq
 (1+C_0\tau)\mathcal E_n
 +C_1\tau\left(h^2+\tau^2+\frac{\delta}{\tau^2}\right).
\end{equation*}
Here $C_0,C_1$ are independent of the grid spacing.
Discrete Gronwall yields
\begin{equation}
 \max_{0\leq n\tau\leq T}\mathcal E_n
 \leq
 \exp(C_0(T+\tau))
 \left(\mathcal E_0+
 C_1(T+\tau)\left(h^2+\tau^2+\frac{\delta}{\tau^2}\right)\right).
 \label{eq:gronwall-energy}
\end{equation}
The exponential factor depends on the physical observation time and is uniform in the spatial mesh.

\emph{Initialization.}
Componentwise rounding gives
$\|e^0\|_{1,h}\leq C\delta(1+h^{-1})$.
Let $e^{-1}=W^{-1}-U_h^{-1}$.
Comparing \eqref{eq:second-order-integer-initialization} with the backward
Taylor expansion of $U$ gives
\begin{equation*}
 \left\|\frac{e^0-e^{-1}}{\tau}\right\|_h
 \leq
 C\left(
 \tau^2+\tau h^2+\frac{\delta}{\tau}
 +\frac{\tau\delta}{h^2}+\tau\delta
 \right).
\end{equation*}
The term $\tau\delta/h^2$ comes from applying $\Delta_h$ to the
rounded initial field.  Since $\tau=\nu h$, it is bounded by a constant
times $\delta/\tau$.
For \eqref{eq:analytic-taylor-initialization}, direct Taylor expansion
instead gives $C(\tau^2+\delta/\tau)$ for the same velocity error,
while the initial field estimate is unchanged.
The positive-energy identity and
$\tau\|A_h^{1/2}\|\leq2$ imply
\begin{equation*}
 \mathcal E_0
 \leq
 C\left(
 \|e^0\|_{1,h}
 +\left\|\frac{e^0-e^{-1}}{\tau}\right\|_h
 \right)
 \leq C\left(h^2+\tau^2+\frac{\delta}{\tau^2}\right).
\end{equation*}
Insert this into \eqref{eq:gronwall-energy} and use the endpoint estimate.
The backward difference of the exact solution differs from
$U_t(\,\cdot\,,(n-\tfrac12)\tau)$ by $O(\tau^2)$.
This proves \eqref{eq:main-convergence-estimate} for the unrestricted
integer lift with globally Lipschitz $f$.

\emph{Locally Lipschitz forces.}
Choose a globally Lipschitz $C^1$ truncation $f_*$ that agrees with
$f$ on the ball $\{|z|\leq M+2\}$.
The exact solution also solves the equation with $f_*$.
The preceding argument applies to the quantized update using $f_*$.
In one spatial dimension there is a constant $C_\ell$, independent of
the grid, such that
\begin{equation*}
 \max_j|z_j|\leq C_\ell\|z\|_{1,h}.
\end{equation*}
Indeed, the discrete mean $\bar z=J^{-1}\sum_jz_j$ satisfies $|\bar z|\leq\ell^{-1/2}\|z\|_h$. Summing nearest-neighbor differences along a path gives $|z_j-z_k|\leq\ell^{1/2}\|D_hz\|_h$. Averaging this inequality over $k$ yields
\begin{equation*}
 \max_j|z_j|\leq\ell^{-1/2}\|z\|_h+\ell^{1/2}\|D_hz\|_h.
\end{equation*}
Thus
\begin{equation*}
 \max_{0\leq n\leq n_T,j}|W_j^n-U(x_j,n\tau)|
 \leq C_\ell C_T\varepsilon_h.
\end{equation*}
Choose $\varepsilon_*$ small enough that this is at most $1/2$;
the corresponding initialization bound controls $W^{-1}$ as well.
The truncated trajectory stays in $\{|z|\leq M+1\}$.
Consequently $f_*(W_j^n)=f(W_j^n)$ at every update.
The sampled field $u_h^0$ and its rounded value $W^0$ lie in the same ball, so the prescribed initializations coincide for $f$ and $f_*$. The truncated and original integer trajectories agree by induction.
This transfers the error estimate to the original nonlinear force.

\emph{Finite modulus.}
The preceding bounds give
\begin{equation*}
 |Q_j^{a,n}|\leq\frac{M+1}{\delta},
 \qquad
 |P_j^{a,n}|
 =
 |Q_j^{a,n}-Q_j^{a,n-1}|
 \leq\frac{2(M+1)}{\delta}.
\end{equation*}
Condition \eqref{eq:sufficient-modulus} places every integer field and
memory register strictly between $-N/2$ and $N/2$, hence within its
centered representative interval.
The reductions modulo $N$ of the integer states satisfy
\eqref{eq:new-update}.  At every step the centered field
representatives equal the integer fields, so the finite-state update
evaluates the same force.  Induction proves agreement of the modular and integer trajectories through the stated final time.
\end{proof}

\begin{remark}[Meaning of the limit]
\label{rem:meaning-quantized-limit}
The theorem approximates a specified smooth solution on any fixed physical
time interval.  Both the local alphabet and the observation step count
change with the mesh.  With $\delta$ proportional to $h^4$, a sufficient
modulus is of order $h^{-4}$; the register length in bits is therefore
of order $\log(1/h)$.
The error constant depends on the fixed physical time interval. Rounding generally makes the finite-state rule nonlinear even when $f$ is linear.  Analytic normal-mode formulas describe the unrounded reference recurrence.
\end{remark}

\begin{corollary}[An energy-based register range]
\label{cor:energy-register-range}
Suppose $f=\nabla V$, $V\geq0$, and the exact solution is smooth.
Let
\begin{equation*}
 \mathcal H_0
 =
 \int_{\mathbb T_\ell}
 \left(
 \frac12|U_t(x,0)|^2+
 \frac12|U_x(x,0)|^2+V(U(x,0))
 \right)\,dx.
\end{equation*}
Then the constant $M$ in \eqref{eq:sufficient-modulus} can be replaced by
\begin{equation*}
 M_T^{\mathrm{en}}
 =
 \ell^{-1/2}
 \left(\|U(\,\cdot\,,0)\|_{L^2}
 +(T+\tau_0)\sqrt{2\mathcal H_0}\right)
 +\ell^{1/2}\sqrt{2\mathcal H_0}.
\end{equation*}
\end{corollary}

\begin{proof}
Multiplication of \eqref{eq:intro-pde} by $U_t$ and
periodic integration by parts show that $\mathcal H(t)=\mathcal H_0$.
Thus
\begin{equation*}
 \|U_x(t)\|_{L^2}\leq\sqrt{2\mathcal H_0},
 \quad
 \|U_t(t)\|_{L^2}\leq\sqrt{2\mathcal H_0},
 \quad
 \|U(t)\|_{L^2}
 \leq\|U(0)\|_{L^2}+|t|\sqrt{2\mathcal H_0}.
\end{equation*}
For a periodic vector-valued function,
$\|z\|_\infty\leq\ell^{-1/2}\|z\|_{L^2}
+\ell^{1/2}\|z_x\|_{L^2}$.
Combining the estimates proves the claim.
\end{proof}

\subsection{An exact nonlinear traveling wave}
\label{subsec:exact-cn-wave}

\begin{proposition}[Quartic traveling waves]
\label{prop:quartic-cn-wave}
Let $\alpha>0$, $A>0$, and let $K(m)$ denote the complete elliptic
integral of the first kind with parameter $m$. Write $\operatorname{cn}$, $\operatorname{sn}$, and $\operatorname{dn}$ for the Jacobi elliptic functions in the same parameter convention.
Set
\begin{equation*}
 \kappa=\frac{4K(1/2)}{\ell},\qquad
 \Omega=\sqrt{\kappa^2+\alpha A^2}.
\end{equation*}
Then
\begin{equation}
 u(x,t)=A\,\operatorname{cn}\left(\kappa x-\Omega t\,\middle|\,\frac12\right)
 \label{eq:quartic-cn-wave}
\end{equation}
is a smooth $\ell$-periodic solution of
\begin{equation*}
 u_{tt}-u_{xx}+\alpha u^3=0.
\end{equation*}
Its initial velocity is
\begin{equation*}
 u_t(x,0)=A\Omega\,
 \operatorname{sn}\left(\kappa x\,\middle|\,\frac12\right)
 \operatorname{dn}\left(\kappa x\,\middle|\,\frac12\right).
\end{equation*}
\end{proposition}

\begin{proof}
In parameter notation, the Jacobi identity is
\begin{equation*}
 \frac{d^2}{dz^2}\operatorname{cn}(z|m)
 =
 (2m-1)\operatorname{cn}(z|m)
 -2m\operatorname{cn}(z|m)^3.
\end{equation*}
At $m=1/2$ it reduces to $\operatorname{cn}''=-\operatorname{cn}^3$.
Substitution of \eqref{eq:quartic-cn-wave} therefore gives
\begin{equation*}
 u_{tt}-u_{xx}+\alpha u^3
 =
 A\bigl(-\Omega^2+\kappa^2+\alpha A^2\bigr)
 \operatorname{cn}\left(\kappa x-\Omega t\,\middle|\,\tfrac12\right)^3
 =0.
\end{equation*}
The real period is $4K(1/2)$.
Finally, $\operatorname{cn}'=-\operatorname{sn}\operatorname{dn}$
gives the stated initial velocity.
\end{proof}

On the covering space, \eqref{eq:quartic-cn-wave} is the Lorentz transform
of the homogeneous quartic oscillator
$A\operatorname{cn}(A\sqrt{\alpha}\,t|1/2)$, with boost speed
$v=\kappa/\Omega<1$.
Its phase velocity $\Omega/\kappa>1$ is not a signal velocity.
Theorem~\ref{thm:finite-alphabet-wave-convergence} gives a convergence test
against this exact solution, with the nonlinear force evaluated at every time step.
The elliptic identities and periods used here are tabulated in the NIST
Digital Library of Mathematical Functions, equation 22.13.14 and Tables 22.13.1 and 22.4.1
\cite{NISTDLMF}. We use the parameter convention $m=k^2$.

\section{Coupled scalar fields and symmetry breaking}
\label{sec:coupled}

\subsection{A common-speed nonlinear limit}

For two scalar fields $(\phi,\chi)$ take $f=\nabla V$, where
\begin{equation}
 \begin{split}
 V(\phi,\chi)={}&\frac{m_\phi^2}{2}\phi^2+
 \frac{m_\chi^2}{2}\chi^2+
 \frac{\eta}{2}(\phi-\chi)^2\\
 &+\frac{\lambda_\phi}{4}\phi^4+
 \frac{\lambda_\chi}{4}\chi^4+
 \frac g2\phi^2\chi^2.
 \end{split}
 \label{eq:coupled-potential}
\end{equation}
The forces are
\begin{align}
 f_\phi&=m_\phi^2\phi+\eta(\phi-\chi)
                    +\lambda_\phi\phi^3+g\phi\chi^2,\\
 f_\chi&=m_\chi^2\chi-\eta(\phi-\chi)
                    +\lambda_\chi\chi^3+g\chi\phi^2.
 \label{eq:coupled-forces}
\end{align}
Nonnegative mass squares, $\eta\ge0$, $\lambda_\phi,\lambda_\chi\ge0$, and $g\ge0$ make $V$ nonnegative. For the continuum field $U=(\phi,\chi)$, the energy
\begin{equation}
 H[U]=\int_0^\ell\left[
       \tfrac12\abs{\partial_tU}^2+
       \tfrac12\abs{\partial_xU}^2+V(U)\right]\dd x
 \label{eq:continuum-energy}
\end{equation}
is conserved along smooth periodic solutions. Section~\ref{sec:numerics} evaluates its discrete analogue along rounded trajectories.

On the covering spacetime the limiting action is
\begin{equation}
 S=\int\dd t\,\dd x\,
 \left[\tfrac12\partial_\mu\phi\,\partial^\mu\phi+
       \tfrac12\partial_\mu\chi\,\partial^\mu\chi-V(\phi,\chi)\right],
 \label{eq:lorentz-action}
\end{equation}
with Minkowski metric $\diag(1,-1)$. Both fields transform as scalars and have the same characteristic cone. The approximation theorem establishes this nonlinear classical system as the limit of the finite-state family. Lorentz transformations act on the continuum fields and their Cauchy data on the covering spacetime.

\subsection{The symmetric invariant manifold}

Set $m_\phi^2=m_\chi^2=0$ and $\lambda_\phi=\lambda_\chi=\lambda$, and use
\begin{equation*}
 u=(\phi+\chi)/\sqrt2,\qquad
 v=(\phi-\chi)/\sqrt2
\end{equation*}
as real coordinates on the lifted fields. The potential is
\begin{equation}
 V(u,v)=\eta v^2+
       \frac{\lambda+g}{8}(u^4+v^4)+
       \frac{3\lambda-g}{4}u^2v^2.
 \label{eq:symmetric-uv-potential}
\end{equation}

\begin{proposition}[Symmetric reduction]
\label{prop:symmetric-reduction}
For equal field and memory registers of the two species at the initial step, the rounded automaton preserves those equalities at every step, for every odd modulus. With equal initial integer representatives, the lifted recurrence has $v=0$. A smooth continuum solution with $v(0)=v_t(0)=0$ also has $v=0$, and its symmetric field satisfies
\begin{equation}
 u_{tt}-u_{xx}+\frac{\lambda+g}{2}u^3=0.
 \label{eq:reduced-quartic-wave}
\end{equation}
The zero-background continuum linearization has $\omega_u^2=k^2$ and $\omega_v^2=k^2+2\eta$.
\end{proposition}
\begin{proof}
The force is equivariant under exchange of the two species. Equal registers therefore receive equal integer kicks and drifts. The same exchange maps $v$ to $-v$ in the continuum equation. Smooth solutions with the same Cauchy data are unique by a wave energy estimate, since the polynomial force is Lipschitz on bounded sets. Exchange-symmetric initial data therefore give $v=0$. Differentiating \eqref{eq:symmetric-uv-potential} at $v=0$ gives \eqref{eq:reduced-quartic-wave}. The two Hessian eigenvalues at the origin are $0$ and $2\eta$.
\end{proof}

The massless $u$ branch results from setting both bare masses to zero; exchange symmetry by itself also permits a quadratic $u^2$ term.

\subsection{Antisymmetric response to unequal quartic couplings}

Keep zero bare masses and set
\begin{equation*}
 \lambda_\phi=\lambda+\varepsilon,\qquad
 \lambda_\chi=\lambda-\varepsilon,
 \qquad 0\le\abs{\varepsilon}<\lambda.
\end{equation*}
The potential acquires
\begin{equation*}
 V_\varepsilon-V_0=
       \frac\varepsilon2(u^3v+uv^3).
\end{equation*}
The exact continuum equations become
\begin{align}
 u_{tt}-u_{xx}
 &+\frac{\lambda+g}{2}u^3+
   \frac{3\lambda-g}{2}uv^2+
   \frac\varepsilon2(3u^2v+v^3)=0,\label{eq:asym-u}\\
 v_{tt}-v_{xx}+2\eta v
 &+\frac{\lambda+g}{2}v^3+
   \frac{3\lambda-g}{2}u^2v+
   \frac\varepsilon2(u^3+3uv^2)=0.\label{eq:asym-v}
\end{align}
For initial data satisfying $v=v_t=0$, the antisymmetric acceleration is $-\varepsilon u^3/2$. It is nonzero wherever $\varepsilon u\ne0$.

For a time-dependent periodic field, write
\begin{equation*}
 \|w\|_{X_T}
 =\sup_{0\leq t\leq T}
 \left(\|w(t)\|_{H^1(\mathbb T_\ell)}
       +\|w_t(t)\|_{L^2(\mathbb T_\ell)}\right).
\end{equation*}

\begin{proposition}[First variation in the quartic coupling]
\label{prop:first-variation}
Fix $0<\varepsilon_0<\lambda$ and smooth initial data independent of $\varepsilon$, with $v(0)=v_t(0)=0$.
Suppose the solutions $(u_\varepsilon,v_\varepsilon)$ exist smoothly on $[0,T]$ for every $|\varepsilon|\leq\varepsilon_0$ and satisfy
\begin{equation*}
 \sup_{|\varepsilon|\leq\varepsilon_0}
 \left(\|u_\varepsilon\|_{X_T}+\|v_\varepsilon\|_{X_T}\right)\leq B.
\end{equation*}
Let $u_0$ be the symmetric solution and let $v_1$ solve
\begin{equation}
 (v_1)_{tt}-(v_1)_{xx}+
 \left(2\eta+\frac{3\lambda-g}{2}u_0^2\right)v_1
 =-\frac12u_0^3,
 \qquad v_1(0)=(v_1)_t(0)=0.
 \label{eq:first-order-generated-field}
\end{equation}
Then
\begin{equation}
 \|u_\varepsilon-u_0\|_{X_T}\leq C|\varepsilon|^2,
 \qquad
 \|v_\varepsilon-\varepsilon v_1\|_{X_T}\leq C|\varepsilon|^3,
 \label{eq:continuum-first-variation-bound}
\end{equation}
where $C$ depends on $B,T,\ell,\eta,\lambda,g,\varepsilon_0$ and the common initial data, but is independent of $\varepsilon$.
\end{proposition}

\begin{proof}
For $w_{tt}-w_{xx}+a(x,t)w=F$ with bounded $a$, multiplying by $w_t$ and adding the derivative of $\|w\|_{L^2}^2/2$ gives, by Gronwall's inequality,
\begin{equation*}
 \|w\|_{X_T}\leq C
 \left(\|w(0)\|_{H^1}+\|w_t(0)\|_{L^2}
             +\int_0^T\|F(t)\|_{L^2}\,dt\right).
\end{equation*}
The same estimate holds for a bounded matrix coefficient. The embedding $H^1(\mathbb T_\ell)\subset L^\infty(\mathbb T_\ell)$ bounds all polynomial coefficients uniformly under the stated hypothesis. Apply this estimate to the difference between the equations at $\varepsilon$ and at zero. The explicit parameter-dependent force is $O(|\varepsilon|)$ in $L^2$, so
$\|u_\varepsilon-u_0\|_{X_T}+\|v_\varepsilon\|_{X_T}\leq C|\varepsilon|$.

In the $u$ equation, factor $u_\varepsilon^3-u_0^3$ into a bounded coefficient times $u_\varepsilon-u_0$. The remaining terms $u_\varepsilon v_\varepsilon^2$ and $\varepsilon(3u_\varepsilon^2v_\varepsilon+v_\varepsilon^3)/2$ are $O(|\varepsilon|^2)$ in $L^2$. The energy estimate gives the first bound in~\eqref{eq:continuum-first-variation-bound}. Next set $r=v_\varepsilon-\varepsilon v_1$ and move
$[2\eta+(3\lambda-g)u_0^2/2]r$ to the left. The remaining sources are constant multiples of
\begin{equation*}
 (u_\varepsilon^2-u_0^2)v_\varepsilon,\qquad
 v_\varepsilon^3,\qquad
 \varepsilon(u_\varepsilon^3-u_0^3),\qquad
 \varepsilon u_\varepsilon v_\varepsilon^2.
\end{equation*}
Each has $L^2$ norm at most $C|\varepsilon|^3$. The initial data for $r$ vanish, so the energy estimate proves the second bound.
\end{proof}

Species exchange and uniqueness also give $u_{-\varepsilon}=u_\varepsilon$ and $v_{-\varepsilon}=-v_\varepsilon$.
To compare with the rounded dynamics, suppose in addition that the continuum fields
$U_\varepsilon=(\phi_\varepsilon,\chi_\varepsilon)$ extend as solutions to $\mathbb T_\ell\times[-\tau_0,T]$ with a common $C^5$ bound.
Write $S_hw=(w(x_j))_j$ and
$v_{h,\varepsilon}^n=(W_\phi^n-W_\chi^n)/\sqrt2$.
Use either initialization in Theorem~\ref{thm:finite-alphabet-wave-convergence} and its Courant and smallness conditions, with a common modulus satisfying
\eqref{eq:sufficient-modulus} for $M=\sup_{\varepsilon,x,t}|U_\varepsilon(x,t)|$.
Then
\begin{equation}
 \max_{0\leq n\leq n_T}
 \|v_{h,\varepsilon}^n-\varepsilon S_hv_1(n\tau)\|_{1,h}
 \leq C\left(|\varepsilon|^3+\varepsilon_h\right),
 \qquad \varepsilon_h=h^2+\tau^2+\delta/\tau^2.
 \label{eq:discrete-first-variation-bound}
\end{equation}
Indeed, the approximation constants are uniform for this family of polynomial forces and uniformly smooth solutions. Also
$\|S_hw\|_{1,h}\leq C_\ell\|w\|_{H^1}$ uniformly in $h$.
The discrete $L^2$ bound follows from the continuous Sobolev embedding; for the difference term, Jensen's inequality gives
$h|(D_hS_hw)_j|^2\leq\int_{x_j}^{x_{j+1}}|w_x|^2\,dx$.
The triangle inequality and~\eqref{eq:continuum-first-variation-bound} prove~\eqref{eq:discrete-first-variation-bound}.
Resolving the leading linear scaling requires $\varepsilon_h\ll|\varepsilon|$; resolving the cubic remainder requires $\varepsilon_h\ll|\varepsilon|^3$.
Rounding can make the trajectory discontinuous in $\varepsilon$ at fixed resolution.

For two initial packets $L$ and $R$, let $v_{L+R}^n,v_L^n,v_R^n$ be the antisymmetric integer-lift fields for the two-packet and one-packet initial data, evolved on the same grid. The quantity
\begin{equation}
 C_v(n\tau)=\norm{v_{L+R}^n-v_L^n-v_R^n}_{2,h}
 \label{eq:collision-v-excess}
\end{equation}
measures the failure of nonlinear superposition. The single-packet responses are subtracted before the norm is taken. We report field amplitudes, the discretized continuum energy, and the quadratic energy diagnostic for the antisymmetric branch separately.

\section{Numerical experiments}
\label{sec:numerics}

\subsection{Arithmetic and error measurements}

The simulations implement \eqref{eq:rounded-force}--\eqref{eq:new-update} with exact rational coefficients and Python arbitrary-precision integers. The main experiments use
\begin{equation*}
 \nu=\frac12,\qquad \delta=\frac{h^4}{64},\qquad N=2^{61}-1.
\end{equation*}
The small-system inversion test and the support test use smaller moduli. Smooth initial fields are sampled in floating-point arithmetic and rounded once according to \eqref{eq:analytic-taylor-initialization}. Thereafter no force decision uses floating point. Every field and memory register at every time level is saved. An independently reduced modular trajectory is computed alongside the unrestricted integer trajectory, and equality is checked at every level. All main trajectories satisfy the no-wrap bound.

Errors are measured on the sampled grid with the norms of Section~\ref{subsec:discrete-wave-energy}. Unless a final-time error is specified, a reported error is the maximum over all stored times. In the exact-solution test the time-difference error is
\begin{equation*}
 \max_{1\le n\le n_T}\left\|
 \frac{W^n-W^{n-1}}{\tau}
 -\frac{U_h^n-U_h^{n-1}}{\tau}\right\|_h.
\end{equation*}
The observed order for an error $e_h$ is $\log_2(e_h/e_{h/2})$. The reported errors are deterministic discretization errors and differences from the reference solver.

For interior times define the continuum-type energy diagnostic
\begin{equation}
 H_h^n=h\sum_j\left[
 \frac12\left|\frac{W_j^{n+1}-W_j^{n-1}}{2\tau}\right|^2
 +\frac12|D_hW_j^n|^2+V(W_j^n)\right].
 \label{eq:numerical-energy-diagnostic}
\end{equation}
Its reported relative variation is $\max_{1\le n<n_T}|H_h^n-H_h^1|/H_h^1$. The diagnostic measures the variation of the discretized continuum energy at times $\tau,\ldots,(n_T-1)\tau$. For the antisymmetric field we additionally report
\begin{equation*}
 E_{v,\mathrm{quad}}^n=\frac h2\sum_j
 \left[|d_t^0v_j^n|^2+|D_hv_j^n|^2+2\eta|v_j^n|^2\right],
\end{equation*}
where $d_t^0v^n=(v^{n+1}-v^{n-1})/(2\tau)$. This is the quadratic energy of the free reference branch. The quartic cross terms remain in the total diagnostic \eqref{eq:numerical-energy-diagnostic}.

The ultraviolet set is $|kh|>\pi/2$. With the unnormalized discrete Fourier transform, the absolute ultraviolet power diagnostic is
\begin{equation}
 A_{\rm UV}^n=
 \frac{(h/J)\sum_{a,\,|kh|>\pi/2}|\widehat W_a^n(k)|^2}
      {h\sum_{a,j}|W_{a,j}^0|^2}.
 \label{eq:absolute-uv-power}
\end{equation}
We also retain the ultraviolet fraction within $v$. The absolute normalization avoids a singular or misleading ratio when $v$ initially vanishes.

\subsection{Convergence to an exact nonlinear traveling wave}
\label{subsec:numerical-cn}

For the solution in Proposition~\ref{prop:quartic-cn-wave}, take $\ell=8$, $A=1$, $\alpha=1$, and $0\le t\le4$. Four meshes, $J=32,64,128,256$, give the errors in Table~\ref{tab:cn-convergence}. The $H_h^1$ errors decrease from $1.42098\times10^{-2}$ to $2.16879\times10^{-4}$, with consecutive observed orders $1.977$, $2.015$, and $2.042$. The time-difference error decreases from $1.57499\times10^{-2}$ to $2.43307\times10^{-4}$.

\begin{table}[t]
\caption{Maximum errors against the exact nonlinear traveling wave over $0\le t\le4$. The last column is the observed $H_h^1$ order.}
\label{tab:cn-convergence}
\centering
\input{tables/convergence_table}
\end{table}

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{convergence.pdf}
\caption{Convergence to the exact nonlinear traveling wave with $\tau=h/2$ and $\delta=h^4/64$. The left panel shows maximum displacement and time-difference errors against the analytic solution. The right panel shows final displacement errors divided by $h^2$.}
\label{fig:cn-convergence}
\end{figure}

An unrounded real-valued recurrence is initialized from the same smooth data and evaluated separately. Its maximum $L_h^2$ difference from the integer trajectory decreases from $8.60269\times10^{-4}$ to $1.15657\times10^{-6}$ over the four meshes. This comparison separates the net effect of initial quantization and force rounding from the ordinary spatial and temporal discretization error. The relative energy variation falls from $2.52000\times10^{-4}$ to $1.11421\times10^{-6}$.

\subsection{Comparison with an independent coupled-field solver}

A separate test uses $\ell=8$, $T=2$, and
\begin{equation*}
 (m_\phi^2,m_\chi^2,\eta,\lambda_\phi,\lambda_\chi,g)
 =\left(\tfrac14,\tfrac12,\tfrac14,\tfrac3{10},\tfrac1{10},\tfrac1{10}\right).
\end{equation*}
Writing $k_0=2\pi/8$, the initial fields are
$\phi_0=0.4\cos(k_0x)+0.1\cos(2k_0x)$,
$\chi_0=0.3\sin(k_0x)$, with velocities
$\dot\phi_0=0.15\sin(k_0x)$ and $\dot\chi_0=-0.1\cos(k_0x)$.
The integer run uses $J=128$. The reference evaluates spatial derivatives by Fourier collocation and integrates the resulting first-order system with DOP853.

The maximum integer-to-reference $L_h^2$ field error is $1.04605\times10^{-4}$ and the final error is $7.45965\times10^{-5}$. For reference solutions on $256$ and $512$ spatial sites, the maximum difference over time on the comparison grid is $1.17331\times10^{-14}$ in the same field norm. Tightening the relative and absolute tolerances from $(2\times10^{-12},2\times10^{-13})$ to $(3\times10^{-13},3\times10^{-14})$, while halving the maximum time step, changes the reference by $1.12274\times10^{-14}$. The reference integrations save both fields and velocities.

\subsection{Light and massive branches under refinement}
\label{subsec:numerical-spectra}

Set both bare masses and all quartic couplings to zero, choose $\eta=1/2$, and take $\ell=32$, $T=16$. The continuum branches have squared masses $0$ and $1$. Each branch is initialized separately with amplitude $0.15$ per Fourier component, modes $m=1,\ldots,4$ for $u$ and $m=0,\ldots,4$ for $v$, and the corresponding continuum traveling-mode velocities. The grids are $J=128,256,512$. At the coarsest grid the largest linear lattice parameter is $4\nu^2+\tau^2=1.015625$, well inside the upper stability boundary $4$.

We project each saved trajectory onto its complex Fourier mode. A recurrence fit supplies an initial frequency estimate. An independent temporal fit then minimizes the full-trajectory residual of $a\cos(\omega t)+b\sin(\omega t)$, with complex amplitudes $a,b$ fitted at each trial frequency. The search interval is $[0.9,1.1]$ times the recurrence estimate and is determined from the data. Repeating the temporal fit on the first and second halves of the observation interval measures time-window sensitivity. The recurrence residual remains a separate implementation diagnostic. The lowest light mode covers only about half a period over $T=16$.

\begin{table}[t]
\caption{Dispersion parameters obtained by temporal sinusoid fitting. For each branch, the frequencies are fitted to $\omega^2=m_{\rm fit}^2+c_{\rm fit}^2k^2$ over the stated modes. The last column is the maximum relative difference from the continuum frequency among those modes.}
\label{tab:speeds}
\centering
\input{tables/speed_table}
\end{table}

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{spectroscopy.pdf}
\caption{Light and massive branches of the quantized linear model. The frequencies are obtained by independent temporal fits to saved trajectories. Grid refinement brings both fitted speeds toward the common continuum value. The analytic lattice relation gives the unrounded frequencies for comparison.}
\label{fig:spectroscopy}
\end{figure}

At $J=512$ the fitted speeds are $c_u=0.99993124$ and $c_v=1.00001936$. The maximum relative differences from the unrounded frequencies are $1.21970\times10^{-5}$ and $4.88019\times10^{-6}$, respectively. The temporal-fit residuals are at most $1.48639\times10^{-6}$ for $u$ and $1.04757\times10^{-5}$ for $v$; the largest relative half-window frequency differences are $8.40973\times10^{-7}$ and $2.10827\times10^{-6}$. The fitted light-branch intercept can include contributions from lattice dispersion, rounding, and frequency fitting. The common continuum characteristic speed is one, as follows from \eqref{eq:reference-expansion} and the principal part of the limiting field equations.

\subsection{Continuous interaction and generated antisymmetric fields}
\label{subsec:numerical-collision}

Take $\ell=32$, $T=12$, $\eta=1/4$, $\lambda=1/5$, $g=1/10$, and $\lambda_\phi=\lambda+\varepsilon$, $\lambda_\chi=\lambda-\varepsilon$. The initial symmetric field consists of two packets. With
\begin{equation*}
 G(z)=0.6\exp[-z^2/(2\cdot1.6^2)]\cos(\pi z/4),
\end{equation*}
write $G_L,G_R$ for its periodic-image sums centered at $10$ and $22$. The implementation sums the nearest three images; the omitted Gaussian tails are below double-precision resolution on this domain. We prescribe
\begin{equation*}
 u(0)=G_L+G_R,\quad u_t(0)=-G_L'+G_R',\quad
 v(0)=v_t(0)=0.
\end{equation*}
The initial $Q$ registers of the two species are exactly equal. The initial memory includes the second-order acceleration term required by \eqref{eq:analytic-taylor-initialization}; it represents the velocity near $-\tau/2$.

The principal case has $\varepsilon=0.1$, $J=512$. Matched left-only and right-only runs use the same model, grid, and initialization rule. Additional runs use $\varepsilon=0.05$, $0.025$, $0$, and $-0.1$, a linear control, and grids $J=256$ and $1024$ for the principal parameter. Single-packet controls are computed at all three meshes for $\varepsilon=0.1$, and at $J=512$ for each of the three positive asymmetries. The nonlinear forces act at every time step.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{collision.pdf}
\caption{Continuously interacting packets for $\varepsilon=0.1$, $J=512$. The antisymmetric field starts from zero. The two one-packet controls distinguish single-packet generation from the nonlinear-superposition difference.}
\label{fig:collision}
\end{figure}

The symmetric nonlinear run and the linear control retain $v=0$ exactly. Changing $\varepsilon$ from $0.1$ to $-0.1$ exactly interchanges the two species' entire saved $Q,P$ trajectories. In the principal asymmetric run, $\max_t\|v(t)\|_h=0.0669694$, compared with $0.0238355$ in either one-packet run. The maximum collision excess \eqref{eq:collision-v-excess} is $0.0440175$. The response includes both single-packet generation and a nonzero superposition difference.

\begin{table}[t]
\caption{Antisymmetric response from symmetric initial data at $J=512$. Maxima are over the saved interval; the quadratic energy is evaluated at interior times. The zero-asymmetry response vanishes exactly.}
\label{tab:epsilon}
\centering
\input{tables/epsilon_table}
\end{table}

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{epsilon_scaling.pdf}
\caption{Parameter scaling in the integer runs at $J=512$. Fits over $\varepsilon=0.025,0.05,0.1$ give an amplitude exponent of $1.00294$ and a quadratic-energy exponent of $2.00254$.}
\label{fig:epsilon}
\end{figure}

The fitted exponents over $\varepsilon=0.025,0.05,0.1$ are $1.00294$ for $\max\|v\|_h$, $2.00254$ for the maximum quadratic branch energy, and $0.99728$ for the maximum collision excess. Table~\ref{tab:collision-refinement} compares complete fields at common grid points and common times. The maximum adjacent-grid difference decreases from $0.00735513$ to $0.00180432$, giving order $2.02729$. The corresponding collision-excess difference in $v$ decreases from $0.00205667$ to $0.000149297$.

\begin{table}[t]
\caption{Mesh comparison for $\varepsilon=0.1$. The fine trajectory is restricted to the coarse space-time grid. The last column gives the norm of the difference between the nonlinear-superposition fields on the two grids.}
\label{tab:collision-refinement}
\centering
\input{tables/collision_refinement_table}
\end{table}

The maximum absolute ultraviolet power \eqref{eq:absolute-uv-power} is $1.66090\times10^{-8}$, $1.40688\times10^{-10}$, and $1.19164\times10^{-12}$ on the three meshes. The relative variation of the total-energy diagnostic is $6.69687\times10^{-4}$, $2.17089\times10^{-4}$, and $4.32243\times10^{-5}$. At $J=1024$, the maximum ultraviolet fraction within $v$ after $t=1$ is $7.17703\times10^{-10}$.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{stability.pdf}
\caption{Energy and ultraviolet diagnostics for the principal collision under refinement. The energy variation is measured relative to the first interior energy sample. The ultraviolet power is normalized by the initial total field power and includes both fields.}
\label{fig:stability}
\end{figure}

\subsection{Independent check of the first variation}

The first-variation equation \eqref{eq:first-order-generated-field} is integrated with the Fourier--DOP853 reference method, together with its symmetric background $u_0$. Independent full continuum solutions use $\varepsilon=0.1,0.05,0.025$ and the same physical initial data. Table~\ref{tab:source-response} gives the maximum over time of the norm of $v_\varepsilon/\varepsilon-v_1$. The fitted exponent for this error is $2.00605$, consistent with the cubic remainder for $v_\varepsilon$ in Section~\ref{sec:coupled}. Spatial refinement of the full reference changes its fields by $1.56100\times10^{-13}$; halving the source solver's maximum time step and tightening tolerances changes the background and first variation by $2.12148\times10^{-14}$.

\begin{table}[t]
\caption{First-variation check. The middle column compares solutions of the independently integrated full continuum and first-variation equations. The last column compares the integer $v$ field at $J=512$ with the full continuum solution.}
\label{tab:source-response}
\centering
\input{tables/source_response_table}
\end{table}

At fixed $J=512$, the absolute integer-to-continuum $v$ errors are approximately $3.02\times10^{-4}$, $3.31\times10^{-4}$, and $7.56\times10^{-4}$ as the asymmetry decreases. The error relative to the response increases as the response decreases at fixed resolution. The combined estimate in Section~\ref{sec:coupled} specifies the joint refinement of $h$ and $\delta$ needed to resolve smaller values of $\varepsilon$. Figure~\ref{fig:source-response} in Appendix~\ref{app:implementation} displays the continuum response comparison.

\subsection{Exact arithmetic, inversion, and support checks}

The arithmetic tests cover $1932$ rational rounding cases including ties, oddness of edge forces, and polynomial numerators with $91$ digits. A two-field system with modulus $257$ is evolved for $37$ steps and then inverted for $37$ steps, with exact recovery even though modular wrapping is allowed. In a separate $129$-site support test, the initial field register is $512$ at one site and zero elsewhere, with zero initial memory. Over $40$ steps no nonzero register lies outside the radius-$n$ classical cone, and the inverse recovers the initial configuration exactly.

An independent validation script evaluates selected force components directly with rational arithmetic, replays saved trajectories forward and backward, checks the ultraviolet Parseval normalization, checks the space-time restriction used in mesh comparisons, and tests the frequency fit against different initial estimates. A separate serialized-output validator checks all $44$ raw data packages and all $9864$ saved transitions in $30$ integer trajectories, with exact agreement in both forward and inverse directions. The operator support statement used in the analysis is proved in Section~\ref{sec:reversibleca}.
\FloatBarrier

\section{Discussion}
\label{sec:discussion}

The finite-state dynamics is exactly invertible, has a causal support bound for the full local matrix algebra, and approximates smooth nonlinear waves with a quantified error. The edge decomposition determines operator support. The strict Courant condition controls the spatial difference operator, and the deterministic rounding bound controls force quantization. Together these estimates yield an error of order $h^2+\tau^2+\delta/\tau^2$ on a fixed physical time interval.

\subsection{Continuum scaling and register size}

For $\tau=\nu h$ and $\delta\asymp h^4$, a sufficient modulus grows as $h^{-4}$. Each register therefore requires $O(\log(1/h))$ bits, and the number of states per site grows polynomially in $h^{-1}$. The finite-time register bound ensures that the modular trajectory realizes the unrestricted integer recurrence throughout the approximation interval. The free nondispersive family in Appendix~\ref{app:auxiliary} gives an explicit case in which the register bound holds for all time.

The present nonlinear estimate applies to smooth solutions in one spatial dimension. The proof uses a mesh-independent estimate of the supremum norm in terms of the $H^1$ norm; this step requires a stronger estimate in higher dimensions. Longer observation times require control of the growth of the error constant and of the corresponding register range. Solutions with lower regularity require consistency estimates adapted to their regularity.

The numerical results agree with the second-order estimate for an exact quartic traveling wave and for coupled fields. The antisymmetric response follows the continuum source equation when its amplitude is resolved by the mesh and field spacing. The condition $h^2+\tau^2+\delta/\tau^2\ll|\varepsilon|$ expresses this resolution requirement, while the cubic remainder requires the stronger comparison with $|\varepsilon|^3$.

\subsection{Permutation dynamics and quantum observables}

The permutation representation supplies a finite-dimensional unitary dynamics and a well-defined algebra of off-diagonal operators. The fields used in the approximation theorem are diagonal configuration observables. Conjugation by a permutation preserves the diagonal algebra, so these particular field observables commute at every pair of times, independently of the chosen state. A continuum quantum field construction would therefore require a noncommuting observable map and convergence of the associated correlation functions.

Such a construction also requires a physical state, an energy generator with an appropriate spectral condition, and a scaling prescription for local operators. On a finite system, logarithms of a permutation unitary have branch ambiguities, and locality of the update alone does not ensure locality of a logarithm. The exact operator support bounds proved here constrain propagation in such a construction.


\appendix
\section{An all-time bounded family for the free wave}
\label{app:auxiliary}

At $\nu=1$ and $f=0$, every edge argument is an integer, so rounding leaves the force unchanged. The following explicit solution treats this endpoint of the Courant range; the convergence theorem assumes $0<\nu<1$.

\begin{proposition}[A global no-wrap family at the nondispersive point]
\label{prop:global-nowrap-free-family}
For the unrounded massless integer rule
\begin{equation*}
 Q_x^{n+1}-2Q_x^n+Q_x^{n-1}
 =
 Q_{x+1}^n-2Q_x^n+Q_{x-1}^n,
\end{equation*}
let $a,b:\mathbb Z/J\mathbb Z\to\mathbb Z$ and set
\begin{equation*}
 Q_x^n=a(x-n)+b(x+n),\qquad
 P_x^n=Q_x^n-Q_x^{n-1}.
\end{equation*}
If $B=\|a\|_\infty+\|b\|_\infty$, then
$|Q_x^n|\leq B$ and $|P_x^n|\leq2B$ for all integer $n$.
Every odd $N>4B$ gives a finite-state trajectory with no wrap for
all time.
\end{proposition}

\begin{proof}
Substitution verifies the recurrence separately for each traveling
component.  The displayed bounds follow from the triangle inequality.
Translation preserves the norms of $a$ and $b$, so the same bounds
hold at every positive and negative time.
\end{proof}

\section{Reproduction and supplementary numerical checks}
\label{app:implementation}

The accompanying sources include the \texttt{cell\_auto} Python package, a rational-parameter YAML configuration, and programs that generate the trajectories, numerical summaries, tables, and figures. After preparing the source directory as described in the README, run the following commands from the project root:
\begin{verbatim}
python -m pip install -r experiment/requirements-lock.txt
python experiment/reproduce.py
python experiment/verify_outputs.py
python experiment/verify_results.py
python make_manuscript_tables.py
latexmk -pdf manuscript.tex
\end{verbatim}
The calculations were run with Python 3.12 and independently reproduced with Python 3.13. Each run records its software versions in \path{experiment/results/environment.json}. The README gives the source layout, software requirements, and reproduction commands. A simulation run stores integer arrays $Q,P$ and physical fields at all time levels, together with the fitted frequencies and numerical summaries. Reference trajectories store fields and physical velocities, and each control is a separate named run. Fixed-width integer storage is used only after an explicit range check; the update itself uses arbitrary-precision integers.

\begin{figure}[htbp]
\centering
\includegraphics[width=\linewidth]{source_response.pdf}
\caption{Continuum first-variation test. The left panel compares $\|v_\varepsilon/\varepsilon\|_h$ with $\|v_1\|_h$. The right panel plots $\max_t\|v_\varepsilon/\varepsilon-v_1\|_h$ against $\varepsilon$, with an observed exponent close to two. These curves come from continuum reference integrations; the separate integer errors are reported in Table~\ref{tab:source-response}.}
\label{fig:source-response}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width=\linewidth]{causal_support.pdf}
\caption{Exact register support in the compact-perturbation experiment. Dashed lines show the one-cell-per-step support bound, corresponding to speed $2$ in the physical units used here. The continuum cone has unit speed in the limiting equation. The full operator support bound is proved in Theorem~\ref{thm:pairwisecone}.}
\label{fig:causal-support}
\end{figure}


\FloatBarrier


\backmatter
\section*{Statements and Declarations}
\noindent\textbf{Funding.} This work was supported by Northeastern University, Grant Number 10011345052601.

\noindent\textbf{Competing interests.} The author declares no competing interests.

\noindent\textbf{Ethics approval and consent.} Not applicable.

\noindent\textbf{Data and code availability.} The accompanying source files contain exact rational configurations and programs that generate the integer trajectories, continuum reference solutions, numerical summaries, tables, and figures. The README specifies the reproduction commands and software requirements.

\noindent\textbf{Author contribution.} Zhenyu Sun is responsible for the conception, mathematical analysis, numerical methodology, interpretation of the results, and preparation of the manuscript.

\noindent\textbf{Use of AI-assisted tools.} ChatGPT was used for proofreading and language editing.

\begingroup
\interlinepenalty=10000
\setlength{\bibsep}{6pt}
\bibliography{references}
\endgroup
\end{document}