\documentclass[pdflatex,sn-mathphys-num]{sn-jnl}
\usepackage{amsmath,amssymb,amsthm,mathtools}
\usepackage{booktabs,array,xcolor}
\usepackage{graphicx,microtype}
\usepackage{enumitem}
\usepackage{etoolbox}
\usepackage{needspace}

% Author manuscript layout; the bundled Springer class is unchanged.
\geometry{reset,a4paper,left=23mm,right=23mm,top=23mm,bottom=24mm,
  bindingoffset=0mm,headheight=14pt,headsep=6mm,footskip=12mm}
\makeatletter
\renewcommand\normalsize{%
  \@setfontsize\normalsize{11pt}{13.5pt}%
  \abovedisplayskip=9pt plus2pt minus2pt
  \belowdisplayskip=\abovedisplayskip
  \abovedisplayshortskip=3pt plus2pt
  \belowdisplayshortskip=5pt plus2pt minus2pt
  \let\@listi\@listI}
\renewcommand\small{\@setfontsize\small{10pt}{12pt}}
\renewcommand\footnotesize{\@setfontsize\footnotesize{9pt}{11pt}}
\renewcommand\Titlefont{\normalfont\fontsize{18pt}{22pt}\bfseries\titraggedcenter}
\renewcommand\addressfont{\normalfont\fontsize{10pt}{12pt}\titraggedcenter}
\renewcommand\abstractheadfont{\normalfont\normalsize\bfseries\titraggedcenter}
\renewcommand\abstractfont{\normalfont\normalsize\leftskip=24pt\rightskip=24pt\parfillskip=0pt plus1fil}
\renewcommand\keywordfont{\normalfont\small\leftskip=0pt\rightskip=0pt\parfillskip=0pt plus1fil}
\renewcommand\figurecaptionfont{\normalfont\small}
\renewcommand\tablecaptionfont{\normalfont\small}
\def\ps@manuscript{%
  \let\@oddhead\@empty\let\@evenhead\@empty
  \def\@oddfoot{\hfil\normalfont\normalsize\thepage\hfil}%
  \let\@evenfoot\@oddfoot}
\let\ps@headings\ps@manuscript
\let\ps@titlepage\ps@manuscript
\let\ps@plain\ps@manuscript
\patchcmd{\@maketitle}{\vskip21pt}{\vskip0pt}{}{}
\patchcmd{\@maketitle}{\vskip20pt}{\vskip6pt}{}{}
\patchcmd{\@maketitle}{\vskip36pt}{\vskip18pt}{}{}
\patchcmd{\@maketitle}{Corresponding author(s). E-mail(s):}{Corresponding author:}{}{}
\patchcmd{\email}{;\ }{}{}{}
\AtBeginDocument{\renewenvironment{proof}[1][\proofname]{%
  \par\pushQED{\qed}\normalfont\normalsize
  \topsep=7pt plus2pt\trivlist
  \item[\hskip\labelsep\itshape #1\@addpunct{.}]\ignorespaces
}{\popQED\endtrivlist\@endpefalse}}
\makeatother
\normalsize
\newtheoremstyle{paperplain}{9pt plus2pt minus1pt}{9pt plus2pt minus1pt}
  {\normalfont\itshape}{0pt}{\bfseries}{.}{.5em}{}
\newtheoremstyle{paperdefinition}{9pt plus2pt minus1pt}{9pt plus2pt minus1pt}
  {\normalfont}{0pt}{\bfseries}{.}{.5em}{}
\newtheoremstyle{paperremark}{9pt plus2pt minus1pt}{9pt plus2pt minus1pt}
  {\normalfont}{0pt}{\itshape}{.}{.5em}{}

\hypersetup{
  pdftitle={Gibbs Approximation from Coarse-Graining of Finite Systems},
  pdfauthor={Sun Zhenyu and Chen Chen}
}

\numberwithin{equation}{section}

\theoremstyle{paperplain}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{paperdefinition}
\newtheorem{definition}[theorem]{Definition}
\theoremstyle{paperremark}
\newtheorem{remark}[theorem]{Remark}
\BeforeBeginEnvironment{theorem}{\Needspace{11\baselineskip}}
\BeforeBeginEnvironment{proposition}{\Needspace{6\baselineskip}}
\BeforeBeginEnvironment{lemma}{\Needspace{6\baselineskip}}
\BeforeBeginEnvironment{corollary}{\Needspace{6\baselineskip}}

\setlist{topsep=4pt,itemsep=2pt,parsep=0pt}
\setlength{\emergencystretch}{2em}
\setlength{\parindent}{1.2em}
\setlength{\bibsep}{4pt}
\renewcommand{\bibfont}{\normalfont\small}
\AtBeginDocument{\raggedbottom}

\begin{document}

\title[Gibbs Approximation from Coarse-Graining]{Gibbs Approximation from Coarse-Graining of Finite Systems}
\author*[1]{Sun Zhenyu}\email{sunzhenyu@neuq.edu.cn}
\author[2]{Chen Chen}
\affil*[1]{\orgdiv{School of Mathematics and Statistics}, \orgname{Northeastern University at Qinhuangdao}, \orgaddress{\city{Qinhuangdao}, \postcode{066004}, \country{China}}}
\affil[2]{\orgdiv{Sanli Honors College, School of Physics and Electronics Engineering}, \orgname{Shanxi University}, \orgaddress{\city{Taiyuan}, \state{Shanxi}, \postcode{030006}, \country{China}}}

\abstract{Thermodynamic integrability in finite systems admits concave probability-response potentials beyond the canonical ensemble. We investigate how replacing an unresolved block of energy levels by its effective free energy constrains this class and determines the Gibbs distribution, both exactly and quantitatively. For a response family defined across all finite spectra, exact single-block replacement implies conditional independence within each block and selects the Gibbs family with its canonical free-energy normalization. On a bounded energy window, conditional log-odds discrepancies determine a common reference inverse-temperature scale and yield a dimension-independent total-variation bound. A discrete calibration protocol incorporates grid rounding, interpolation under stated regularity bounds, and finite-sample concentration guarantees. We bound the conditional log-odds discrepancies in terms of coarse-grained block-probability residuals through finite block suppression, a bounded outer-control displacement, and a local internal curvature bound. The estimate separates auxiliary system size, measurement precision, and probability conditioning. Smooth local perturbations show that finitely many observations do not uniquely determine the global response. Applications to finite-capacity reservoirs and fixed-coupling mean-force models distinguish work-form closure from temperature-energy scaling. For two subsystems sharing small, bounded, symmetric log-temperature fluctuations, the leading mutual information is half the squared log-temperature variance times the product of their reference canonical heat capacities.}

\keywords{Gibbs distribution, coarse-graining, thermodynamic integrability, probability ratios, finite reservoirs, temperature fluctuations}
\maketitle

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

The canonical Gibbs distribution describes the equilibrium populations of a system weakly coupled to a thermal reservoir. Thermodynamic integrability constrains their response to changes in the energy levels by requiring the existence of free energy and entropy. For finite systems, these integrability conditions also admit noncanonical probability responses. We ask whether consistency under coarse-graining selects the Gibbs distribution, and what quantitative conclusion remains when system sizes, energy ranges, and measurement precision are finite.

Zhai and Sun~\cite{ZhaiSun2024} formulated the entropy-integrability equations (EIE) for finite diagonal systems and constructed canonical and noncanonical solutions, including temperature-scale mixtures and correlations induced by a shared fluctuating temperature. We distinguish the closure condition for the reversible work differential from the scaling relation between temperature and energy. Finite-dimensional convex duality~\cite{Rockafellar1970,HiriartUrrutyLemarechal2001} associates each spectrally stable response potential with a concave thermodynamic entropy, without prescribing the Shannon form.

Several complementary approaches in statistical physics, information theory, and decision theory characterize exponential weights through consistency or independence principles. Classical axiomatizations of Shannon entropy rely on recursive grouping identities~\cite{Shannon1948a,Faddeev1956,BaezFritzLeinster2011}, which connect to equilibrium statistical mechanics via Jaynes's maximum-entropy variational principle~\cite{Jaynes1957a}. Directly at the level of state probabilities, Cerreia-Vioglio et al.~\cite{CerreiaVioglioEtAl2021} introduced axiomatic tests for the Boltzmann distribution. In the language of stochastic choice, where utility corresponds to negative energy, Sandomirskiy and Tamuz~\cite{SandomirskiyTamuz2024} characterized the multinomial logit rule and its quantitative stability under continuity, decomposability, and neutrality, and established a related characterization of the Boltzmann law via convolutional independence of energy distributions~\cite{SandomirskiyTamuz2025}. A later revision by Sandomirskiy, Sung, Tamuz, and Wincelberg, \emph{A Separability Foundation for Random Coefficients Logit} (\href{https://arxiv.org/abs/2312.04827v4}{arXiv:2312.04827v4}, 2026), characterized random-coefficients mixtures of linear logit models under marginal consistency, monotonicity, and continuity. Fu et al.~\cite{FuEtAl2026} derived exponential weights from invariance under coarse-graining scale transformations and uniform energy shifts.

The coarse-graining operation studied in the present work replaces a selected block of energy levels by the value of its own response potential. Operationally, this prescription specifies both the aggregate probability of the unresolved block and the reference zero of its effective free energy. When imposed across a family of response potentials defined for all finite system dimensions, this single-block replacement identity already implies conditional independence within the resolved block. The canonical Gibbs distribution, the Shannon entropy formula, and the Shannon grouping law then emerge as direct consequences, without requiring separate assumptions on conditional factorization or on the density of attainable probability vectors.

Our main quantitative results address what remains of this characterization when experiments are confined to a fixed energy window and a finite collection of system dimensions. On any bounded energy interval, three-state probability-ratio cycles control the additivity defect of binary log-odds within the same window. This yields a single endpoint-calibrated Gibbs scale, a dimension-independent bound in total-variation distance, and a discrete calibration protocol requiring only finitely many spectral preparations. To connect these log-odds bounds to errors in coarse-grained block probabilities, we introduce a finite suppression protocol that transfers control from an auxiliary configuration, in which the target block carries small probability mass, to the configuration of interest. Because probability residuals control only the outer gradient of the potential defect, recovering the internal conditional law involves a trade-off between the outer control displacement and a local bound on the second internal derivative, without requiring derivatives of the empirical measurement error.

Section~\ref{sec:response} introduces response potentials and entropy duality. The exact characterization is proved in Section~\ref{sec:exact}; finite-window estimates and their connection to coarse-graining errors follow in Sections~\ref{sec:finite-window} and~\ref{sec:coarse-error}. Section~\ref{sec:physical-models} examines reservoir models and shared temperature fluctuations. The appendices discuss finite-observation non-identifiability, sparse transition graphs, and numerical implementation.


\section{Response potentials and entropy}
\label{sec:response}

Let $\Delta_N=\{p\in\mathbb{R}^N:p_i\geq0,\ \sum_i p_i=1\}$ denote the $(N-1)$-dimensional probability simplex and $\Delta_N^\circ$ its relative interior. We consider diagonal state populations $P(\beta,E)$ parameterized by independently variable energy levels $E=(E_1,\ldots,E_N)\in\mathbb{R}^N$ and inverse temperature $\beta>0$. The entropy-integrability equations (EIE) read~\cite{ZhaiSun2024}
\begin{align}
\partial_{E_j}P_i&=\partial_{E_i}P_j,
\label{eq:closure}\\
\beta\partial_\beta P_i&=\sum_jE_j\partial_{E_j}P_i.
\label{eq:descent}
\end{align}
Equation~\eqref{eq:closure} expresses the closure of the reversible work differential $\sum_i P_i\,\mathrm{d}E_i$ at fixed $\beta$, whereas~\eqref{eq:descent} requires invariance under the reciprocal rescaling $(\beta,E)\mapsto(\lambda\beta,E/\lambda)$, which preserves the dimensionless spectrum $x=\beta E$. When a physical apparatus gives access only to a lower-dimensional submanifold of control parameters, \eqref{eq:closure}--\eqref{eq:descent}, if satisfied on an ambient spectral domain, are sufficient for integrability on the accessible submanifold; integrability of the restricted response alone need not imply these ambient equations.

\begin{proposition}[Local potential and physical normalization]
\label{prop:normal}
Any $C^1$ solution of~\eqref{eq:closure}--\eqref{eq:descent} on a sufficiently small convex domain of dimensionless spectra admits the local representation $P(\beta,E)=\nabla\Psi(\beta E)$. Probability normalization $\sum_i P_i=1$ is equivalent to translation covariance along the uniform direction $\mathbf{1}=(1,\ldots,1)^\top$,
\begin{equation}
\Psi(x+t\mathbf{1})=\Psi(x)+t,
\label{eq:translation}
\end{equation}
whenever the segment connecting $x$ and $x+t\mathbf{1}$ lies within the domain. The reversible work form is generated by the free energy $F(\beta,E)=\beta^{-1}\Psi(\beta E)+B(\beta)$. Requiring the thermodynamic internal energy $U=\partial_\beta(\beta F)$ to coincide with the spectral expectation $E\cdot P$ fixes $B(\beta)=C/\beta$ for a constant $C$. The associated thermodynamic entropy takes the Legendre form
\begin{equation}
S_{\mathrm{th}}=x\cdot\nabla\Psi(x)-\Psi(x)+C_0,
\qquad \mathrm{d}S_{\mathrm{th}}=\beta\sum_i E_i\,\mathrm{d}P_i.
\label{eq:clausius}
\end{equation}
\end{proposition}
\begin{proof}
The characteristic vector field $\beta\partial_\beta-E\cdot\nabla_E$ of~\eqref{eq:descent} has coordinate invariants $x_i=\beta E_i$, so $P(\beta,E)=Q(x)$ locally. By~\eqref{eq:closure}, the differential $1$-form $Q(x)\cdot\mathrm{d}x$ is closed, and the Poincar\'e lemma on a convex domain yields $Q=\nabla\Psi$ for a scalar potential $\Psi$. Integrating the normalization identity $\mathbf{1}\cdot\nabla\Psi(x)=1$ along $\mathbf{1}$ gives~\eqref{eq:translation}. Integrating $\partial_{E_i}F=P_i$ with respect to $E$ produces $F(\beta,E)=\beta^{-1}\Psi(\beta E)+B(\beta)$. Differentiating $\beta F$ with respect to $\beta$ gives $\partial_\beta(\beta F)=E\cdot P+\frac{\mathrm{d}}{\mathrm{d}\beta}(\beta B(\beta))$, so imposing $U=E\cdot P$ requires $(\beta B)'=0$. Finally, $\mathrm{d}\Psi=P\cdot\mathrm{d}x$ implies $\mathrm{d}(x\cdot P-\Psi)=x\cdot\mathrm{d}P$. Combined with the first-law identity $\mathrm{d}U=P\cdot\mathrm{d}E+E\cdot\mathrm{d}P$, this recovers the Clausius relation~\eqref{eq:clausius} for the heat differential $\delta Q=E\cdot\mathrm{d}P$.
\end{proof}

In our global analysis, an $N$-state response potential is a finite $C^1$ function $\Psi_N:\mathbb{R}^N\to\mathbb{R}$ satisfying~\eqref{eq:translation} with positive gradient $\nabla\Psi_N\in\Delta_N^\circ$. The potential is symmetric if it is invariant under coordinate permutations, and spectrally stable if it is concave on $\mathbb{R}^N$ (or equivalently on the zero-sum hyperplane $H_N=\mathbf{1}^\perp$). When $\Psi_N$ is $C^2$, translation covariance and concavity imply
\begin{equation}
\operatorname{Hess}\Psi_N\,\mathbf{1}=0,
\qquad v^\top\operatorname{Hess}\Psi_N\,v\leq0\quad(\forall v\in\mathbb{R}^N),
\qquad C=-x^\top\operatorname{Hess}\Psi_N(x)\,x\geq0,
\label{eq:stability}
\end{equation}
so that energy levels and populations respond monotonically in opposition ($\partial_{x_i}P_i\le 0$) and the heat capacity $C$ is nonnegative.

\begin{proposition}[Concave entropy duality]
\label{prop:duality}
Let $\Psi:\mathbb{R}^N\to\mathbb{R}$ be a finite, concave, translation-covariant potential whose supergradients lie in $\Delta_N$. Then the conjugate function
\begin{equation}
s_\Psi(p)=\inf_{x\in\mathbb{R}^N}\{x\cdot p-\Psi(x)\},
\qquad
\Psi(x)=\min_{p\in\Delta_N}\{x\cdot p-s_\Psi(p)\},
\label{eq:duality}
\end{equation}
defines a proper, concave, upper semicontinuous entropy $s_\Psi:\Delta_N\to\mathbb{R}\cup\{-\infty\}$. Conversely, every proper, concave, upper semicontinuous function $s:\Delta_N\to\mathbb{R}\cup\{-\infty\}$ generates a finite, concave, translation-covariant potential $\Psi$ via the second relation in~\eqref{eq:duality}, and at any differentiability point $x$ of $\Psi$ the unique minimizer is $p=\nabla\Psi(x)$. Moreover, if $s_\Psi$ is finite throughout $\Delta_N^\circ$ and $\Psi$ is $C^1$ with $\nabla\Psi(x)\in\Delta_N^\circ$, then $\nabla\Psi$ is surjective onto the open simplex:
\begin{equation}
\nabla\Psi(\mathbb{R}^N)=\Delta_N^\circ.
\label{eq:onto}
\end{equation}
\end{proposition}
\begin{proof}
As an infimum of affine functions, $s_\Psi$ is concave and upper semicontinuous. Since $\Psi$ is finite and concave on $\mathbb{R}^N$, it admits a supporting supergradient $p\in\Delta_N$ at every point $x$, which gives $s_\Psi(p)=x\cdot p-\Psi(x)>-\infty$; combining this with the definition of $s_\Psi$ establishes the variational recovery of $\Psi(x)$ in~\eqref{eq:duality}. Conversely, any proper upper semicontinuous concave function $s$ on the compact simplex $\Delta_N$ is bounded above. For each $x$, the lower semicontinuous objective $p\mapsto x\cdot p-s(p)$ is bounded below and finite somewhere, so it attains a finite minimum in~\eqref{eq:duality}. This yields a finite, concave, translation-covariant potential $\Psi$ whose supergradients coincide with the minimizers $p$. Differentiability of $\Psi$ at $x$ renders the supergradient unique. Finally, if $s_\Psi$ is finite on $\Delta_N^\circ$, concave duality furnishes a supporting supergradient $x\in\mathbb{R}^N$ at every $p\in\Delta_N^\circ$, meaning $p$ minimizes $x\cdot q-s_\Psi(q)$ and hence $p=\nabla\Psi(x)$.
\end{proof}

Proposition~\ref{prop:duality} is the concave specialization of Fenchel--Moreau duality~\cite{Rockafellar1970,HiriartUrrutyLemarechal2001} adapted to the simplex via translation covariance. Whenever $P=\nabla\Psi(x)$, the dual entropy $s_\Psi(P)=x\cdot P-\Psi(x)$ coincides with the Clausius entropy~\eqref{eq:clausius} up to the integration constant $C_0$. For reference throughout the paper, we denote the canonical Gibbs response at scale parameter $a\in\mathbb{R}$ by
\begin{equation}
\mathsf{G}_{a,N,i}(x)=\frac{e^{-ax_i}}{\sum_j e^{-ax_j}}.
\label{eq:gibbs}
\end{equation}
Spectrally stable Gibbs responses correspond to $a\ge 0$, with $a=0$ giving the uniform distribution. For $a>0$, the canonically normalized potential and its dual entropy are
\begin{equation}
\Psi_{a,N}(x)=-a^{-1}\log\sum_j e^{-ax_j},\qquad
s_{a,N}(p)=-a^{-1}\sum_i p_i\log p_i,
\label{eq:gibbsduals}
\end{equation}
where $a$ sets the inverse-temperature scale relative to $\beta$.


\section{Exact replacement of an unresolved block}
\label{sec:exact}

Suppose now that a symmetric, spectrally stable response potential $\Psi_N$ is specified for every integer $N\geq1$, normalized by $\Psi_1(t)=t$, and write $P_i^N=\partial_i\Psi_N$. Consider an $n$-state outer spectrum $x\in\mathbb{R}^n$ in which the $a$-th state ($1\le a\le n$) is resolved into an internal block of $m$ sublevels with relative spectrum $y\in\mathbb{R}^m$. The resolved $(n+m-1)$-state spectrum and the corresponding $n$-state coarse-grained spectrum are defined by
\begin{align}
S_a(x,y)&=(x_1,\ldots,x_{a-1},x_a+y_1,\ldots,x_a+y_m,x_{a+1},\ldots,x_n),\label{eq:splitfield}\\
T_a(x,y)&=(x_1,\ldots,x_{a-1},x_a+\Psi_m(y),x_{a+1},\ldots,x_n).
\label{eq:coarsefield}
\end{align}
To track the outer states in the resolved spectrum, let $\kappa_a(i)=i$ for $i<a$ and $\kappa_a(i)=i+m-1$ for $i>a$.

\begin{definition}[Single-block free-energy replacement]
\label{def:replacement}
The family $\{\Psi_N\}_{N\ge 1}$ satisfies exact free-energy replacement if, for all $n,m\ge 1$, every $1\le a\le n$, and all spectra $x\in\mathbb{R}^n$ and $y\in\mathbb{R}^m$,
\begin{align}
P_{\kappa_a(i)}^{n+m-1}(S_a(x,y))&=P_i^n(T_a(x,y)),\qquad i\ne a,
\label{eq:replacement-outside}\\
\sum_{b=1}^m P_{a+b-1}^{n+m-1}(S_a(x,y))&=P_a^n(T_a(x,y)).
\label{eq:replacement-mass}
\end{align}
\end{definition}

Although~\eqref{eq:replacement-mass} follows from~\eqref{eq:replacement-outside} by probability conservation, we state both relations explicitly because the aggregate block probability is a primary coarse-grained observable. Note that the value $\Psi_m(y)$ enters directly into the coarse-grained level $x_a+\Psi_m(y)$ in~\eqref{eq:coarsefield}; consequently, shifting $\Psi_m$ by an additive constant alters the replacement protocol even though it leaves the isolated $m$-state gradient $P^m(y)=\nabla\Psi_m(y)$ unchanged.

\begin{theorem}[Exact free-energy replacement]
\label{thm:exact-replacement}
A family of symmetric, spectrally stable response potentials $\{\Psi_N\}_{N\ge 1}$ satisfies exact free-energy replacement if and only if there exists a universal constant $a_*>0$ such that, for every $N\ge 1$,
\begin{equation}
\Psi_N(x)=-a_*^{-1}\log\sum_{i=1}^N e^{-a_*x_i}.
\label{eq:exact-gibbs-potential}
\end{equation}
In particular, conditional independence of internal block populations, the Shannon form of the dual entropy, and the Shannon grouping identity all follow as consequences of Definition~\ref{def:replacement}.
\end{theorem}

\begin{proof}
Fix $n,m,a$ and compare the resolved and coarse-grained potentials
\[
F(x,y)=\Psi_{n+m-1}(S_a(x,y)),\qquad
G(x,y)=\Psi_n(T_a(x,y)).
\]
By the chain rule, the replacement identities~\eqref{eq:replacement-outside}--\eqref{eq:replacement-mass} are equivalent to $\nabla_x(F-G)=0$, which implies $F(x,y)-G(x,y)=K(y)$ for a function $K$ independent of $x$. Differentiating this relation with respect to $y_b$ yields
\begin{equation}
\partial_{y_b}K(y)=
P_{a+b-1}^{n+m-1}(S_a(x,y))-P_a^n(T_a(x,y))P_b^m(y).
\label{eq:exact-K}
\end{equation}
For $n\geq2$, choose an outer index $j\ne a$ and consider an auxiliary experiment in which state $j$ of an $n$-state zero spectrum is split into $k$ degenerate zero-energy sublevels. By permutation symmetry of the resulting $(n+k-1)$-state zero spectrum, every resolved state carries probability $(n+k-1)^{-1}$. Applying~\eqref{eq:replacement-outside} to this auxiliary split shows that the coarse-grained spectrum $t^{(k)}\in\mathbb{R}^n$ with $t_j^{(k)}=\Psi_k(0)$ and $t_i^{(k)}=0$ ($i\ne j$) satisfies
\[
P_a^n(t^{(k)})=\frac{1}{n+k-1}\xrightarrow[k\to\infty]{}0.
\]
Now, for any fixed $y\in\mathbb{R}^m$, choose $x^{(k)}\in\mathbb{R}^n$ by setting $x_i^{(k)}=t_i^{(k)}$ for $i\ne a$ and $x_a^{(k)}=t_a^{(k)}-\Psi_m(y)$, so that $T_a(x^{(k)},y)=t^{(k)}$. Since both $P_{a+b-1}^{n+m-1}(S_a(x^{(k)},y))$ and $P_a^n(t^{(k)})P_b^m(y)$ lie in $[0,P_a^n(t^{(k)})]$, the right-hand side of~\eqref{eq:exact-K} is bounded in absolute value by $(n+k-1)^{-1}$. Because the left-hand side of~\eqref{eq:exact-K} is independent of $x$, letting $k\to\infty$ forces $\partial_{y_b}K(y)=0$, which establishes the conditional factorization
\begin{equation}
P_{a+b-1}^{n+m-1}(S_a(x,y))=P_a^n(T_a(x,y))P_b^m(y).
\label{eq:derived-split-independence}
\end{equation}
(For $n=1$, \eqref{eq:derived-split-independence} holds trivially by translation covariance.)

Specializing~\eqref{eq:derived-split-independence} to binary splits ($m=2$) and invoking permutation symmetry to group any pair of states $(i,j)$ of an $N$-state system ($N\ge 3$) into a two-state block shows that all pairwise log-odds depend only on the energy gap through the universal binary function $h$:
\[
\log\frac{P_i^N(z)}{P_j^N(z)}=h(z_i-z_j),\qquad
h(t)=\log\frac{P_1^2(t,0)}{P_2^2(t,0)}.
\]
Applying the identity $\frac{P_1^3}{P_3^3}=\frac{P_1^3}{P_2^3}\frac{P_2^3}{P_3^3}$ to a three-state spectrum gives Cauchy's functional equation $h(s+t)=h(s)+h(t)$ for all $s,t\in\mathbb{R}$. Since $P^2$ is continuous, $h(t)=-a_*t$ for a single constant $a_*\in\mathbb{R}$ common to all dimensions. Concavity of $\Psi_2$ implies that $P_1^2(t,0)$ is nonincreasing in $t$, so $a_*\ge 0$. If $a_*=0$, every response is uniform; taking $n=m=2$ in~\eqref{eq:replacement-outside} would then equate the resolved outer probability $1/3$ with the coarse-grained probability $1/2$, a contradiction. Hence $a_*>0$.

Integrating $\nabla\Psi_m=\mathsf{G}_{a_*,m}$ yields $\Psi_m(y)=-a_*^{-1}\log\sum_{b=1}^m e^{-a_*y_b}+k_m$ for a dimension-dependent constant $k_m$. Applying~\eqref{eq:replacement-outside} to a two-state zero outer spectrum with one state split into $m$ zero-energy sublevels gives
\[
\frac{1}{m+1}=\frac{1}{1+e^{-a_*\Psi_m(0)}},
\]
so $\Psi_m(0)=-a_*^{-1}\log m$ and thus $k_m=0$ for all $m$. The converse is immediate from the partition-function identity $\sum_{b=1}^m e^{-a_*(x_a+y_b)}=e^{-a_*(x_a+\Psi_m(y))}$.
\end{proof}

\begin{corollary}[Entropy and grouping]
\label{cor:exact-grouping}
Under the hypotheses of Theorem~\ref{thm:exact-replacement}, the dual entropy on the closed simplex $\Delta_N$ is $s_N(p)=-a_*^{-1}\sum_{i=1}^N p_i\log p_i$ (with $0\log 0=0$), and any partitioned distribution $(p_{ab})$ with block sums $q_a=\sum_b p_{ab}$ satisfies the Shannon grouping law
\begin{equation}
s_{\sum_a n_a}((p_{ab}))
=s_k(q)+\sum_{a:q_a>0}q_a s_{n_a}((p_{ab}/q_a)_b).
\label{eq:grouping-consequence}
\end{equation}
\end{corollary}
\begin{proof}
Evaluating the Fenchel conjugate~\eqref{eq:duality} of $\Psi_N$ in~\eqref{eq:exact-gibbs-potential} yields $s_N(p)=-a_*^{-1}\sum_i p_i\log p_i$, and substituting $\log p_{ab}=\log q_a+\log(p_{ab}/q_a)$ for $q_a>0$ directly produces~\eqref{eq:grouping-consequence}.
\end{proof}

Faddeev-type characterizations assume a recursive splitting rule on probability distributions~\cite{Faddeev1956}; Section~4 of Baez, Fritz, and Leinster~\cite{BaezFritzLeinster2011} reviews this formulation. Theorem~\ref{thm:exact-replacement} derives both conditional independence and the grouping law from free-energy replacement in spectral space. The proof also identifies the role of unbounded auxiliary dimensions: splitting an external state into $k\to\infty$ degenerate sublevels provides a sequence of configurations in which the tested block probability $P_a^n(t^{(k)})$ becomes arbitrarily small.

\subsection{Response independence versus potential additivity}
\label{subsec:product}

For two subsystems with independent spectral controls $x\in\mathbb{R}^N$ and $y\in\mathbb{R}^M$, the additive composite spectrum has levels $(x\oplus y)_{ij}=x_i+y_j$. Strict response independence requires the composite distribution to factorize across all spectral pairs:
\begin{equation}
P_{ij}^{NM}(x\oplus y)=P_i^N(x)P_j^M(y).
\label{eq:product-response}
\end{equation}
As a finite-spectrum counterpart to the decomposability theorems of Sandomirskiy and Tamuz~\cite{SandomirskiyTamuz2024,SandomirskiyTamuz2025}, strict response independence fixes the Gibbs family up to free-energy constants.

\begin{proposition}[Strict product response independence]
\label{prop:product-gibbs}
Let $\{P^N\}_{N\ge 1}$ be a family of continuous, strictly positive, permutation-equivariant responses on all finite spectra. If~\eqref{eq:product-response} holds for all $N,M\ge 1$, then $P^N=\mathsf{G}_{a,N}$ for a universal constant $a\in\mathbb{R}$ (with $a\geq0$ whenever $P^N$ is generated by a concave potential). Unlike Definition~\ref{def:replacement}, condition~\eqref{eq:product-response} leaves the additive constants of the generating potentials $\Psi_N$ undetermined.
\end{proposition}
\begin{proof}
Fix $x\in\mathbb{R}^N$ and distinct indices $i\ne j$, set $d=x_i-x_j$, and couple the system to a binary auxiliary spectrum $y=(0,d)$. Because the composite states $(i,1)$ and $(j,2)$ have identical total energy $x_i+0=x_j+d$, permutation symmetry of $P^{2N}$ and~\eqref{eq:product-response} imply
\[
P_i^N(x)P_1^2(0,d)=P_j^N(x)P_2^2(0,d).
\]
Every pairwise probability ratio therefore depends solely on the energy gap via the binary ratio $P_2^2(0,d)/P_1^2(0,d)$, and the three-state cycle identity forces this dependence to be exponential as in the proof of Theorem~\ref{thm:exact-replacement}.
\end{proof}

By contrast, requiring the generating potentials to be additive under the product embedding,
\begin{equation}
\Psi_{NM}(x\oplus y)=\Psi_N(x)+\Psi_M(y)+C_{N,M},
\label{eq:product-potential}
\end{equation}
constrains only the directional derivatives along the subspace $\{x\oplus y\}$. Differentiating~\eqref{eq:product-potential} with respect to $x_i$ and $y_j$ yields marginal consistency, $\sum_j P_{ij}^{NM}(x\oplus y)=P_i^N(x)$ and $\sum_i P_{ij}^{NM}(x\oplus y)=P_j^M(y)$, but does not enforce the joint factorization~\eqref{eq:product-response}. Indeed, consider the two-scale canonical mixture
\begin{equation}
\Psi_N^\gamma(x)=\sum_{r\in\{r_1,r_2\}}
\frac{\gamma_r}{r}\left(-\log\sum_{i=1}^N e^{-rx_i}\right),
\qquad 0<r_1<r_2,\quad \gamma_r>0,\quad\gamma_{r_1}+\gamma_{r_2}=1.
\label{eq:finite-mixture-potential}
\end{equation}
For every $N$, $\Psi_N^\gamma$ is smooth, symmetric, translation covariant, and strictly concave on $H_N$, and satisfies $\min_i x_i-(\log N)\sum_r\gamma_r/r \le \Psi_N^\gamma(x)\le\min_i x_i$, so its dual entropy is finite throughout the closed simplex $\Delta_N$. Because each canonical summand is additive on $x\oplus y$, $\Psi_N^\gamma$ satisfies~\eqref{eq:product-potential} with $C_{N,M}=0$. Yet for two identical nondegenerate binary systems with gap $z>0$, setting $f_r=(1+e^{rz})^{-1}$ gives
\[
P_{11}^4(x\oplus x)-P_1^2(x)P_1^2(x)
=\sum_r\gamma_r f_r^2-\left(\sum_r\gamma_r f_r\right)^2>0,
\]
violating~\eqref{eq:product-response}. Since differentiability of $\Psi_4^\gamma$ makes its supergradient unique, the Fenchel inequality at $z=x\oplus x$ is strict for the product distribution $p\otimes p$ with $p=P^2(x)$, so
\begin{equation}
s_4(p\otimes p)<z\cdot(p\otimes p)-\Psi_4^\gamma(z)=2s_2(p).
\label{eq:strict-mixture-entropy}
\end{equation}
Additivity of the response potential on product spectra is thus strictly weaker than statistical independence of the subsystem responses.


\section{Finite-window Gibbs approximation}
\label{sec:finite-window}

Theorem~\ref{thm:exact-replacement} assumes exact identities over all spectra and arbitrarily large auxiliary dimensions. Measurements instead give probabilities at finitely many spectra within a bounded energy window. We use two- and three-state calibrations of energy diameter at most $R>0$, together with the conditional log-odds of the target system, to bound its total-variation distance from a common Gibbs law. The bound has no explicit dependence on the target dimension.

For $R>0$, denote the domain of spectra with energy diameter at most $R$ by
\[
 \mathcal{X}_N(R)=\Bigl\{x\in\mathbb{R}^N:\max_i x_i-\min_i x_i\le R\Bigr\}.
\]
Let $P^N:\mathcal{X}_N(R)\to\Delta_N^\circ$ be permutation-equivariant and translation-invariant responses. We define the isolated binary log-odds function and its endpoint-calibrated inverse-temperature slope by
\begin{equation}
 h(t)=\log\frac{P_1^2(t,0)}{P_2^2(t,0)},\qquad
 a_R=-\frac{h(R)}{R}.
 \label{eq:window-binary-calibration}
\end{equation}
Permutation symmetry ensures that $h(-t)=-h(t)$, and concavity of a generating binary potential $\Psi_2$ guarantees that $h(t)$ is nonincreasing with $a_R\ge 0$. At any spectrum $x\in\mathcal{X}_N(R)$, we measure the deviation of pairwise conditional probabilities from the isolated binary law via the conditional log-odds discrepancy
\begin{equation}
 d_N(x)=\max_{i\ne j}
 \left|\log\frac{P_i^N(x)}{P_j^N(x)}-h(x_i-x_j)\right|,
 \label{eq:window-conditional-error}
\end{equation}
which vanishes identically for $N=2$.

To control the deviation of $h$ from linearity on $[-R,R]$, we first record an elementary stability lemma for Cauchy's equation on a bounded interval (compare~\cite{Hyers1941,BorosPalesVolkmann2000}) together with a sharp conversion from pairwise log-odds bounds to total-variation distance.

\begin{lemma}[Additivity on a bounded interval]
\label{lem:local-cauchy}
Let $h:[-R,R]\to\mathbb{R}$ be continuous and odd. If
\[
 |h(s+t)-h(s)-h(t)|\le\delta
 \qquad\forall s,t\in[-R,R]\text{ with }s+t\in[-R,R],
\]
then
\begin{equation}
 \left|h(t)-\frac{h(R)}{R}t\right|\le2\delta,
 \qquad |t|\le R.
 \label{eq:local-cauchy-bound}
\end{equation}
\end{lemma}

\begin{proof}
For $t\in[0,R]$, define the local average
\[
 A(t)=\frac{1}{R}\int_{-R}^0\bigl(h(s+t)-h(s)\bigr)\,\mathrm{d}s.
\]
The approximate additivity hypothesis immediately implies $|A(t)-h(t)|\le\delta$. Moreover, continuity of $h$ allows us to differentiate $A(t)$ explicitly:
\[
 A'(t)=\frac{h(t)-h(t-R)}{R}.
\]
Applying the additivity bound to the decomposition $(t-R)+R=t$ yields $|A'(t)-h(R)/R|\le\delta/R$. Integrating from $t=0$ (where $A(0)=0$) gives $|A(t)-h(R)t/R|\le (t/R)\delta\le\delta$, and combining this with $|A(t)-h(t)|\le\delta$ establishes~\eqref{eq:local-cauchy-bound} on $[0,R]$. The bound on $[-R,0]$ follows by oddness.
\end{proof}

\begin{lemma}[Probability distance from pairwise log-odds]
\label{lem:odds-distance}
Let $p,q\in\Delta_N^\circ$. If
\[
 \max_{i,j}\left|\log\frac{p_i}{p_j}-\log\frac{q_i}{q_j}\right|\le B,
\]
then
\begin{equation}
 \|p-q\|_1\le2\tanh(B/4),\qquad
 \operatorname{TV}(p,q):=\tfrac{1}{2}\|p-q\|_1\le\tanh(B/4).
 \label{eq:odds-tv}
\end{equation}
\end{lemma}

\begin{proof}
Set $r_i=p_i/q_i$, $m=\min_i r_i$, and $M=\max_i r_i$, so that $m\le 1\le M$ and $M/m\le e^B$. If $m=M=1$, the bound is trivial. Otherwise, convexity of $r\mapsto|r-1|$ on $[m,M]$ together with $\sum_i q_i=1$ and $\sum_i q_i r_i=1$ yields
\[
 \operatorname{TV}(p,q)
 =\frac{1}{2}\sum_{i=1}^N q_i|r_i-1|
 \le\frac{(M-1)(1-m)}{M-m}.
\]
For a fixed ratio $r=M/m\ge 1$, substituting $M=rm$ shows that the right-hand side is maximized at $m=r^{-1/2}$, where it equals $(\sqrt{r}-1)/(\sqrt{r}+1)=\tanh((\log r)/4)\le\tanh(B/4)$.
\end{proof}

\begin{theorem}[Finite-interval Gibbs estimate]
\label{thm:finite-window-gibbs}
Suppose $h$ is continuous on $[-R,R]$ and the three-state conditional discrepancy satisfies
\[
 d_3(z)\le\varepsilon_3\qquad\forall z\in\mathcal{X}_3(R).
\]
Then for every $N\ge 2$ and every target spectrum $x\in\mathcal{X}_N(R)$ with $d_N(x)\le\varepsilon_N$, setting $B_N=\varepsilon_N+6\varepsilon_3$ yields
\begin{align}
 \left|\log\frac{P_i^N(x)}{P_j^N(x)}+a_R(x_i-x_j)\right|
 &\le B_N,\label{eq:finite-window-odds}\\
 \|P^N(x)-\mathsf{G}_{a_R,N}(x)\|_1
 &\le2\tanh(B_N/4).\label{eq:finite-window-l1}
\end{align}
Here the reference scale $a_R$ is determined entirely by the binary response at gap $R$, with $a_R\ge 0$ when $\Psi_2$ is concave and $a_R>0$ whenever $h(R)<0$.
\end{theorem}

\begin{proof}
For any $s,t\in[-R,R]$ with $s+t\in[-R,R]$, the three-state spectrum $z=(s+t,t,0)$ belongs to $\mathcal{X}_3(R)$. The pairwise log-odds of $P^3(z)$ sum to zero around the cycle $(1\to 2\to 3\to 1)$. Approximating each log-ratio by $h$ therefore gives
\[
 |h(s+t)-h(s)-h(t)|\le3\varepsilon_3.
\]
Lemma~\ref{lem:local-cauchy} then gives $|h(u)+a_Ru|\le6\varepsilon_3$ for all $u\in[-R,R]$. Applying the triangle inequality to~\eqref{eq:window-conditional-error} proves~\eqref{eq:finite-window-odds}, and Lemma~\ref{lem:odds-distance} yields~\eqref{eq:finite-window-l1}.
\end{proof}

Theorem~\ref{thm:finite-window-gibbs} uses only spectra within the prescribed energy window. When $P^N=\nabla\Psi_N$ throughout $\mathcal{X}_N(R)$, integration along line segments from the origin also bounds the free-energy variation.

\begin{corollary}[Free-energy variation]
\label{cor:finite-window-potential}
Suppose $P^N=\nabla\Psi_N$ and the hypotheses of Theorem~\ref{thm:finite-window-gibbs} hold uniformly for all $x\in\mathcal{X}_N(R)$ with error bound $B_N$. Define the centered reference potential
\[
 \Phi_{a,N}(x)=
 \begin{cases}
 -a^{-1}\log\bigl(N^{-1}\sum_{i=1}^N e^{-ax_i}\bigr),&a\ne0,\\
 N^{-1}\sum_{i=1}^N x_i,&a=0.
 \end{cases}
\]
Then for every $x\in\mathcal{X}_N(R)$,
\begin{equation}
 |\Psi_N(x)-\Psi_N(0)-\Phi_{a_R,N}(x)|
 \le\bigl(\max_i x_i-\min_i x_i\bigr)\tanh(B_N/4).
 \label{eq:finite-window-potential}
\end{equation}
\end{corollary}

\begin{proof}
For any $x\in\mathcal{X}_N(R)$, the ray $tx$ ($0\le t\le 1$) remains within $\mathcal{X}_N(R)$. Since any two probability vectors $p,q\in\Delta_N$ satisfy $|(p-q)\cdot x|\le\operatorname{TV}(p,q)(\max_i x_i-\min_i x_i)$, integrating $(\nabla\Psi_N(tx)-\mathsf{G}_{a_R,N}(tx))\cdot x$ over $t\in[0,1]$ and applying~\eqref{eq:finite-window-l1} yields~\eqref{eq:finite-window-potential}.
\end{proof}

\subsection{Discrete calibration on a finite spectral grid}
\label{subsec:finite-calibration}

In an actual experiment or numerical simulation, one prepares only a finite number of spectra rather than the continuum $\mathcal{X}_3(R)$. We therefore formulate a discrete counterpart to Theorem~\ref{thm:finite-window-gibbs} on a uniform calibration grid of mesh size $\Delta=R/M$ ($M\ge 2$).

\begin{lemma}[Additivity on a finite grid]
\label{lem:grid-cauchy}
Let $M\ge2$ and let $h_{-M},\ldots,h_M\in\mathbb{R}$ satisfy $h_{-k}=-h_k$. If
\[
 |h_{j+k}-h_j-h_k|\le\delta
 \qquad\forall j,k\ge1\text{ with }j+k\le M,
\]
then
\begin{equation}
 \left|h_k-\frac{k}{M}h_M\right|\le2\delta\qquad(-M\le k\le M).
 \label{eq:grid-cauchy-bound}
\end{equation}
\end{lemma}

\begin{proof}
By oddness, $|h_{j+k}-h_j-h_k|\le\delta$ extends to all integers $j,k\in[-M,M]$ with $j+k\in[-M,M]$. For $0\le k\le M$, define the discrete running average
\[
 A_k=\frac{1}{M}\sum_{s=-M}^{-1}(h_{s+k}-h_s).
\]
Then $A_0=0$, $|A_k-h_k|\le\delta$, and telescoping the sum gives the increment
\[
 A_k-A_{k-1}=\frac{h_{k-1}-h_{k-1-M}}{M}.
\]
Applying approximate additivity to the decomposition $(k-1-M)+M=k-1$ shows that $|(A_k-A_{k-1})-h_M/M|\le\delta/M$. Summing over $1\le j\le k$ yields $|A_k-(k/M)h_M|\le (k/M)\delta\le\delta$, which combines with $|A_k-h_k|\le\delta$ to give~\eqref{eq:grid-cauchy-bound}.
\end{proof}

The prefactor $2$ in~\eqref{eq:grid-cauchy-bound} is sharp uniformly in $M$ under endpoint calibration: the odd sequence defined by $h_0=h_M=0$ and $h_k=3k/M-1$ for $1\le k<M$ has additivity defect $1$ on the positive grid while attaining $h_{M-1}=2-3/M\to 2$ as $M\to\infty$. Subtracting the linear trend $3k/M$ renders the sequence nonincreasing without altering either its additivity defect or its deviation from the endpoint-calibrated line, showing that monotonicity does not reduce this worst-case constant.

To calibrate the response and test a finite collection of target spectra $\mathcal{T}$, we prepare:
\begin{enumerate}
 \item the $M$ binary spectra $(0,k\Delta)$ for $1\le k\le M$;
 \item the $\binom{M}{2}$ ternary spectra $(0,j\Delta,k\Delta)$ for $1\le j<k\le M$;
 \item each target spectrum $x\in\mathcal{T}$ (with $\max_i x_i-\min_i x_i\le R$).
\end{enumerate}
The total number of distinct spectral preparations is at most
\begin{equation}
 \mathcal{E}=M+\binom{M}{2}+|\mathcal{T}|.
 \label{eq:grid-operation-count}
\end{equation}
Let $h_k=h(k\Delta)$ for $1\le k\le M$, with $h_0=0$ and $h_{-k}=-h_k$, and let $\varepsilon_3^{\mathrm{grid}}$ denote the maximum conditional log-odds discrepancy across the $\binom{M}{2}$ ternary grid preparations. For any $t\in[-R,R]$, let $\kappa(t)\in\{-M,\ldots,M\}$ be the integer nearest to $t/\Delta$ (breaking ties oddly), and define the discrete target discrepancy
\[
 d_{N,M}^{\mathrm{grid}}(x)=
 \max_{i\ne j}
 \left|\log\frac{P_i^N(x)}{P_j^N(x)}-h_{\kappa(x_i-x_j)}\right|,
\]
which compares each pairwise log-odds in the target spectrum against the nearest calibrated binary gap.

\begin{theorem}[Gibbs estimate from finitely many preparations]
\label{thm:finite-grid-gibbs}
Suppose the binary response is generated by a differentiable concave potential ($a_R\ge 0$). For every target spectrum $x\in\mathcal{T}\cap\mathcal{X}_N(R)$, setting
\[
 B_{N,M}(x)=d_{N,M}^{\mathrm{grid}}(x)+6\varepsilon_3^{\mathrm{grid}}+\frac{a_R\Delta}{2}
\]
yields
\begin{equation}
 \|P^N(x)-\mathsf{G}_{a_R,N}(x)\|_1
 \le2\tanh\bigl(B_{N,M}(x)/4\bigr).
 \label{eq:finite-grid-l1}
\end{equation}
If all pairwise energy gaps of $x$ lie on the calibration grid $\Delta\mathbb{Z}$, the rounding term $a_R\Delta/2$ drops out of $B_{N,M}(x)$.
\end{theorem}

\begin{proof}
For any $j,k\ge 1$ with $j+k\le M$, the ternary preparation $(0,j\Delta,(j+k)\Delta)$ yields $|h_{j+k}-h_j-h_k|\le 3\varepsilon_3^{\mathrm{grid}}$, so Lemma~\ref{lem:grid-cauchy} gives $|h_m+a_Rm\Delta|\le 6\varepsilon_3^{\mathrm{grid}}$ for all $-M\le m\le M$. Since $|\kappa(x_i-x_j)\Delta-(x_i-x_j)|\le\Delta/2$ and $a_R\ge 0$, the triangle inequality bounds the pairwise log-odds discrepancy relative to $\mathsf{G}_{a_R,N}(x)$ by $B_{N,M}(x)$, and Lemma~\ref{lem:odds-distance} completes the proof.
\end{proof}

Theorem~\ref{thm:finite-grid-gibbs} bounds the error at the tested spectra $x\in\mathcal{T}$. A Lipschitz bound on the target log-odds and a finite $\ell^\infty$-net extend this estimate to the full window $\mathcal{X}_N(R)$.

\begin{corollary}[Uniform extension from a finite spectral net]
\label{cor:finite-net-extension}
Fix $N\ge 2$ and suppose $\mathcal{T}_N\subset\mathcal{X}_N(R)$ is a finite $\nu$-net up to translations, in the sense that for every $x\in\mathcal{X}_N(R)$ there exist $z\in\mathcal{T}_N$ and $c\in\mathbb{R}$ with $\|x-z-c\mathbf{1}\|_\infty\le\nu$. Assume further that the log-odds functions $\ell_{ij}(x)=\log(P_i^N(x)/P_j^N(x))$ are $L_\ell$-Lipschitz with respect to $\|\cdot\|_\infty$ on $\mathcal{X}_N(R)$ (for instance, $L_\ell\le 2K/m$ whenever $P_i^N\ge m>0$ and $\|\nabla P_i^N\|_1\le K$). If $|\ell_{ij}(z)+a_R(z_i-z_j)|\le B_*$ holds for all $z\in\mathcal{T}_N$, then every $x\in\mathcal{X}_N(R)$ satisfies
\begin{equation}
 \|P^N(x)-\mathsf{G}_{a_R,N}(x)\|_1
 \le2\tanh\left(\frac{B_*+(L_\ell+2|a_R|)\nu}{4}\right).
 \label{eq:finite-net-extension}
\end{equation}
In particular, fixing each coordinate in turn at zero and letting the remaining $N-1$ coordinates range over $\{0,\Delta_0,\ldots,R\}$ with $\Delta_0=R/M_0$ constructs such a net of cardinality at most $N(M_0+1)^{N-1}$ with mesh radius $\nu=\Delta_0/2$.
\end{corollary}

\begin{proof}
By translation invariance, choosing $z\in\mathcal{T}_N$ and $c\in\mathbb{R}$ with $\|x-z-c\mathbf{1}\|_\infty\le\nu$ gives $|\ell_{ij}(x)-\ell_{ij}(z)|\le L_\ell\nu$ and $|(x_i-x_j)-(z_i-z_j)|\le 2\nu$. Combining these inequalities with the net bound $B_*$ and applying Lemma~\ref{lem:odds-distance} yields~\eqref{eq:finite-net-extension}.
\end{proof}

\subsection{Finite-sample concentration guarantees}
\label{subsec:finite-sample}

When the response probabilities of each of the $\mathcal{E}$ preparations are estimated from $n$ independent experimental trials, the discrete protocol admits an explicit finite-sample confidence bound. Let $V$ denote the total number of outcome coordinates summed across all $\mathcal{E}$ preparations. For a prescribed failure probability $\zeta\in(0,1)$, set
\[
 t_n=\sqrt{\frac{\log(2V/\zeta)}{2n}},
\]
and let $\widehat{p}$ denote the empirical frequency vectors, with $\widehat{m}$ the minimum empirical probability observed across all $V$ coordinates. Whenever $\widehat{m}>t_n$, define
\[
 \Lambda_n=\frac{2t_n}{\widehat{m}-t_n}.
\]
From the empirical frequencies, form the plug-in estimators $\widehat{h}_k=\log(\widehat{p}_2/\widehat{p}_1)$ on the binary preparations $(0,k\Delta)$ (with $\widehat{h}_0=0$ and $\widehat{h}_{-k}=-\widehat{h}_k$), the empirical ternary discrepancy $\widehat{\varepsilon}_3^{\mathrm{grid}}$, the empirical target discrepancy $\widehat{d}_{N,M}^{\mathrm{grid}}(x)$, and the empirical slope
\[
 \widehat{a}=-\frac{\widehat{h}_M}{R},\qquad
 \widehat{a}_+=\max\{0,\widehat{a}\}.
\]

\begin{corollary}[Finite-sample guarantee]
\label{cor:grid-sampling}
Under the assumptions of Theorem~\ref{thm:finite-grid-gibbs}, suppose each preparation is sampled $n$ times independently. With probability at least $1-\zeta$, on the event $\widehat{m}>t_n$, every target spectrum $x\in\mathcal{T}\cap\mathcal{X}_N(R)$ satisfies
\begin{align}
 \|P^N(x)-\mathsf{G}_{\widehat{a}_+,N}(x)\|_1
 &\le2\tanh\bigl(\widehat{B}_{N,M}(x)/4\bigr),\label{eq:sample-grid-l1}\\
 \widehat{B}_{N,M}(x)
 &=\widehat{d}_{N,M}^{\mathrm{grid}}(x)
 +6\widehat{\varepsilon}_3^{\mathrm{grid}}
 +\frac{\widehat{a}_+\Delta}{2}+2\Lambda_n,
 \label{eq:sample-grid-error}
\end{align}
where $2\Lambda_n$ may be tightened to $\Lambda_n$ whenever $\widehat{a}\ge 0$, and the rounding term $\widehat{a}_+\Delta/2$ vanishes if all pairwise gaps of $x$ lie on $\Delta\mathbb{Z}$.
\end{corollary}

\begin{proof}
By Hoeffding's inequality~\cite{Hoeffding1963} and a union bound over all $V$ outcome coordinates, the event $\max|\widehat{p}_i-p_i|\le t_n$ holds with probability at least $1-\zeta$. On this event, whenever $\widehat{m}>t_n$, the elementary bound $|\log(\widehat{p}_i/\widehat{p}_j)-\log(p_i/p_j)|\le 2t_n/(\widehat{m}-t_n)=\Lambda_n$ holds for every empirical log-odds. Because the empirical ternary frequencies satisfy the three-state ratio cycle identically, applying Lemma~\ref{lem:grid-cauchy} directly to the sequence $\widehat{h}_k$ yields $|\widehat{h}_k+\widehat{a}k\Delta|\le 6\widehat{\varepsilon}_3^{\mathrm{grid}}$ without incurring an extra sampling penalty at the algebraic step. Passing from the empirical target log-odds to the true target log-odds adds at most $\Lambda_n$. Finally, if $\widehat{a}<0$, concavity of $\Psi_2$ implies $h(R)\le 0$, so $|\widehat{a}_+-\widehat{a}|R=|\widehat{h}_M|\le|\widehat{h}_M-h(R)|\le\Lambda_n$; clipping the slope to $\widehat{a}_+=0$ therefore adds at most one further $\Lambda_n$. Applying Lemma~\ref{lem:odds-distance} completes the proof.
\end{proof}

The condition $\widehat{m}>t_n$ is sufficient for the Hoeffding confidence intervals used here to have positive lower endpoints, allowing a finite uniform log-odds error bound. It is not a necessary threshold for other confidence procedures; rare outcomes nevertheless make log-odds estimation less precise.

\subsection{Numerical illustrations}
\label{subsec:finite-window-numerics}

To illustrate the discrete calibration protocol of Theorem~\ref{thm:finite-grid-gibbs}, we take $R=3$ and $M=24$, corresponding to $24$ binary and $276$ ternary calibration preparations, and evaluate four target spectra of dimensions $3,4,6,8$ (for $\mathcal{E}=304$ preparations in total). Expressed in units of $R$, the target spectra are $(0,0.243,1)$, $(0,0.25,0.5,1)$, $(0,0.04,0.07,0.35,0.877,1)$, and $(j/7)_{j=0}^7$. We compare four response families: the exact Gibbs response at unit scale ($a=1$), equal-weight two-scale mixtures with scales $(0.85,1.15)$ and $(0.5,1.5)$, and a three-scale mixture $(0.4,1,2)$ with weights $(0.25,0.5,0.25)$. The resulting endpoint-calibrated slopes $a_R$ are $1$, $0.970434$, $0.744790$, and $0.804275$, respectively. Across the four target spectra, the maximum true total-variation distances from $\mathsf{G}_{a_R,N}$ are $0$, $0.007390$, $0.064792$, and $0.068035$, while the corresponding upper bounds from Theorem~\ref{thm:finite-grid-gibbs} are $0.029367$, $0.105157$, $0.575179$, and $0.603586$ (Figure~\ref{fig:finite-window}, left). For the exact Gibbs model, the on-grid four-state target $(0,0.25,0.5,1)R$ exhibits zero bound and zero error. The maximum off-grid bound $0.029367$, attained at the six-state target, combines the rounded-target discrepancy $d_{6,24}^{\mathrm{grid}}=0.055$ and the explicit allowance $a_R\Delta/2=0.0625$, giving $\tanh((0.055+0.0625)/4)$ despite zero actual error.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{figures/finite_window_bounds.pdf}
\caption{Finite-preparation Gibbs bounds and sharpness of the grid additivity constant. Left: largest true total-variation distance and the corresponding upper bound from Theorem~\ref{thm:finite-grid-gibbs} across the four target spectra for each response family. Right: worst-case deviation from the endpoint-calibrated line over odd sequences with unit additivity defect, computed by linear programming; at $M=32$ the maximum is $1.90625$, approaching the sharp asymptotic constant $2$ of Lemma~\ref{lem:grid-cauchy}.}
\label{fig:finite-window}
\end{figure}

Figure~\ref{fig:finite-sample} illustrates the finite-sample bound of Corollary~\ref{cor:grid-sampling} for the weak two-scale mixture $(0.85,1.15)$ with $R=2$ and $M=8$ ($\mathcal{E}=40$ preparations and $V=121$ outcome coordinates), using multinomial samples of size $n\in\{10^4,10^5,10^6,10^7\}$ per preparation at confidence level $1-\zeta=0.95$. As $n$ increases, the sampling penalty $\Lambda_n$ decays as $n^{-1/2}$ and the empirical upper bound converges toward the deterministic finite-grid bound.

\begin{figure}[t]
\centering
\includegraphics[width=0.82\linewidth]{figures/finite_sample_bounds.pdf}
\caption{Finite-sample Gibbs approximation bounds for the equal mixture of scales $0.85$ and $1.15$ ($R=2, M=8$). Curves display maxima across the four target spectra of the true total-variation distance to $\mathsf{G}_{\widehat{a}_+,N}$ and the $95\%$ confidence bound of Corollary~\ref{cor:grid-sampling} as functions of the number of repetitions $n$ per preparation.}
\label{fig:finite-sample}
\end{figure}


\section{Coarse-graining error bounds}
\label{sec:coarse-error}

The stability theorems of Section~\ref{sec:finite-window} take conditional log-odds discrepancies $d_N(x)$ as their input. We now establish how $d_N(x)$ can itself be bounded from coarse-grained block-probability residuals using only a finite family of symmetric, spectrally stable $C^2$ response potentials.

Consider a binary split ($m=2$) of outer state $a$ in an $n$-state system, producing a resolved system of dimension $N=n+1$. The potential defect between the resolved and coarse-grained descriptions is
\begin{equation}
K(x,y)=\Psi_{n+1}(S_a(x,y))-\Psi_n(T_a(x,y)),\qquad x\in\mathbb{R}^n,\ y\in\mathbb{R}^2,
\label{eq:approx-K}
\end{equation}
and the observable coarse-graining probability residuals are
\begin{equation}
r_i(x,y)=
\begin{cases}
P_{\kappa_a(i)}^{n+1}(S_a(x,y))-P_i^n(T_a(x,y)),&i\ne a,\\[4pt]
\displaystyle\sum_{b=1}^2P_{a+b-1}^{n+1}(S_a(x,y))-P_a^n(T_a(x,y)),&i=a.
\end{cases}
\label{eq:probability-residuals}
\end{equation}
By the chain rule, $r_i(x,y)=\partial_{x_i}K(x,y)$. Thus, experimental bounds on the probability residuals $r_i$ directly control the outer gradient $\nabla_x K$, whereas the conditional distribution inside the split block is governed by the internal derivative $\partial_{y_1}K$. The following theorem bridges these two derivatives by combining an outer control displacement to a block-suppressing configuration $x^*$ with a finite-difference estimate in the internal coordinate $y_1$.

\begin{theorem}[Finite-control conditional estimate]
\label{thm:finite-control}
Fix a target configuration $(x,y)\in\mathbb{R}^n\times\mathbb{R}^2$, a suppressing outer field $x^*\in\mathbb{R}^n$ at $\ell^1$ distance $D=\|x-x^*\|_1$, and an internal half-width $h_0>0$. Assume that:
\begin{enumerate}[label=\textup{(\roman*)}]
\item the probability residuals satisfy $|r_i|\le e$ for all $1\le i\le n$ throughout the two-dimensional control rectangle
\begin{equation}
\mathcal{C}=\bigl\{(x+s(x^*-x),y+t e_1):0\le s\le1,\ |t|\le h_0\bigr\};
\label{eq:control-rectangle}
\end{equation}
\item along the two boundary lines $z\in\{x,x^*\}$ and $|t|\le h_0$, the internal second derivative is bounded by
\begin{equation}
\left|\frac{\mathrm{d}^2}{\mathrm{d}t^2}K(z,y+t e_1)\right|\le L;
\label{eq:internal-curvature}
\end{equation}
\item at the suppressing field $x^*$, the coarse-grained block probability satisfies $u^*=P_a^n(T_a(x^*,y))\le\eta$;
\item at the target configuration $(x,y)$, the resolved block mass $M$ and the internal conditional probabilities $q,v$ satisfy the conditioning bounds
\begin{equation}
\begin{aligned}
M&=\sum_{b=1}^2P_{a+b-1}^{n+1}(S_a(x,y))\ge\mu>0,\\
q&=\frac{P_a^{n+1}(S_a(x,y))}{M},\qquad
v=P_1^2(y),\qquad q,v\in[\rho,1-\rho],
\end{aligned}
\label{eq:conditioning}
\end{equation}
for some $\rho\in(0,1/2]$.
\end{enumerate}
Then, for every step size $h\in(0,h_0]$,
\begin{equation}
\left|\log\frac{P_a^{n+1}(S_a(x,y))}{P_{a+1}^{n+1}(S_a(x,y))}
       -\log\frac{P_1^2(y)}{P_2^2(y)}\right|
\le
\frac{\eta+2e+De/h+Lh}{\mu\rho(1-\rho)}.
\label{eq:finite-control-odds}
\end{equation}
In particular, whenever $De>0$, $L>0$, and $\sqrt{De/L}\le h_0$, optimizing over $h$ at $h=\sqrt{De/L}$ yields the numerator
\begin{equation}
\eta+2e+2\sqrt{DLe}.
\label{eq:square-root-error}
\end{equation}
\end{theorem}

\begin{proof}
Define the outer difference $E(t)=K(x,y+t e_1)-K(x^*,y+t e_1)$ for $|t|\le h_0$. Integrating $\nabla_x K=r$ along the segment connecting $x^*$ to $x$ gives $|E(t)|\le D e$ for all $|t|\le h_0$, while assumption~(ii) gives $|E''(t)|\le 2L$. Taylor's theorem with integral remainder around $t=0$ then bounds the central difference:
\[
\left|E'(0)-\frac{E(h)-E(-h)}{2h}\right|
\le\frac{1}{2h}\int_{-h}^h|E'(0)-E'(t)|\,\mathrm{d}t\le Lh,
\qquad
|E'(0)|\le \frac{De}{h}+Lh.
\]
At the suppressing configuration $x^*$, let $p^*=P_a^{n+1}(S_a(x^*,y))$. Because the resolved block mass at $x^*$ differs from $u^*$ by at most $e$, we have $0\le p^*\le u^*+e$. Differentiating $K$ with respect to $y_1$ at $(x^*,y)$ gives $\partial_{y_1}K(x^*,y)=p^*-u^*v$, whence $|\partial_{y_1}K(x^*,y)|\le \eta+e$. Combining this with $|E'(0)|=|\partial_{y_1}K(x,y)-\partial_{y_1}K(x^*,y)|$ yields
\[
|\partial_{y_1}K(x,y)|\le\eta+e+\frac{De}{h}+Lh.
\]
At the target configuration $(x,y)$, writing $u=P_a^n(T_a(x,y))$, we have $\partial_{y_1}K(x,y)=Mq-uv$ and $|M-u|\le e$. Consequently,
\[
|q-v|=\frac{|Mq-uv+(u-M)v|}{M}\le\frac{\eta+2e+De/h+Lh}{M}
\le\frac{\eta+2e+De/h+Lh}{\mu}.
\]
Since the derivative of the logit map $w\mapsto\log(w/(1-w))$ on $[\rho,1-\rho]$ is bounded by $[\rho(1-\rho)]^{-1}$, the mean-value theorem gives~\eqref{eq:finite-control-odds}.
\end{proof}

The square-root scaling $\sqrt{DLe}$ in~\eqref{eq:square-root-error} reflects the ill-posedness of numerical differentiation when only zero-order probability residuals $r_i$ are measured. By the chain rule, $\frac{\mathrm{d}^2}{\mathrm{d}t^2}\Psi_n(T_a(z,y+t e_1))=\partial_{aa}\Psi_n(T_a)(P_1^2)^2+P_a^n(T_a)\partial_{11}\Psi_2$, so if the diagonal Hessian entries of $\Psi_{n+1},\Psi_n,\Psi_2$ are bounded in magnitude by $H_{n+1},H_n,H_2$, the internal curvature constant may be chosen as
\begin{equation}
L=H_{n+1}+H_n+H_2.
\label{eq:curvature-from-potentials}
\end{equation}

To realize the suppression bound $u^*\le\eta$ constructively using a finite auxiliary system, we employ a finite-dimensional version of the degenerate refinement used in the proof of Theorem~\ref{thm:exact-replacement}.

\begin{proposition}[Suppression via finite degenerate refinement]
\label{prop:finite-suppression}
Fix $n\ge2$, an outer index $j\ne a$, and an integer $k\ge2$. If replacing $k$ degenerate zero-energy states in an $(n+k-1)$-state zero spectrum incurs an outer probability residual of at most $e_0$, then the coarse-grained spectrum $t^*\in\mathbb{R}^n$ with $t_j^*=\Psi_k(0)$ and $t_i^*=0$ ($i\ne j$) satisfies
\begin{equation}
P_a^n(t^*)\le\frac{1}{n+k-1}+e_0.
\label{eq:finite-suppression-mass}
\end{equation}
Consequently, for any internal binary spectrum $y\in\mathbb{R}^2$, setting $x_i^*=t_i^*$ ($i\ne a$) and $x_a^*=t_a^*-\Psi_2(y)$ ensures $T_a(x^*,y)=t^*$, so that condition~\textup{(iii)} of Theorem~\ref{thm:finite-control} holds with $\eta=(n+k-1)^{-1}+e_0$.
\end{proposition}
\begin{proof}
By permutation symmetry, every state of the $(n+k-1)$-state zero spectrum has probability $(n+k-1)^{-1}$, and the replacement residual bounds $|P_a^n(t^*)-(n+k-1)^{-1}|\le e_0$.
\end{proof}

Proposition~\ref{prop:finite-suppression} makes explicit the finite resource trade-off underlying block suppression: reducing the residual mass $(n+k-1)^{-1}$ requires preparing an auxiliary system of $n+k-1$ states and setting the $j$-th outer level to $\Psi_k(0)$. This control requirement enters through $D=\|x-x^*\|_1$ and may increase as $k$ grows. Combining Theorem~\ref{thm:finite-control} with the finite-window theorems of Section~\ref{sec:finite-window} yields a complete chain from coarse-graining residuals to Gibbs approximation.

\begin{corollary}[Gibbs approximation from coarse-graining residuals]
\label{cor:coarse-to-gibbs}
Fix an energy window $R>0$ and a finite set of dimensions containing $2$ and $3$. Suppose every pair of levels in each tested spectrum $x\in\mathcal{X}_N(R)$ ($N\ge 3$) can be represented as a binary split satisfying Theorem~\ref{thm:finite-control} with dimension-dependent bounds $\eta_N,e_N,D_N,L_N,h_N,\mu_N,\rho_N$, and set
\begin{equation}
\theta_N=
\inf_{0<h\le h_N}
\frac{\eta_N+2e_N+D_Ne_N/h+L_Nh}
{\mu_N\rho_N(1-\rho_N)},\qquad \theta_2=0.
\label{eq:coarse-theta}
\end{equation}
If these bounds hold across $\mathcal{X}_3(R)$, then every tested $x\in\mathcal{X}_N(R)$ satisfies
\begin{equation}
\left\|P^N(x)-\mathsf{G}_{a_R,N}(x)\right\|_1
\le2\tanh\left(\frac{\theta_N+6\theta_3}{4}\right),
\qquad a_R=-\frac{1}{R}\log\frac{P_1^2(R,0)}{P_2^2(R,0)}\ge0.
\label{eq:coarse-gibbs-bound}
\end{equation}
Similarly, if the three-state tests are performed only on the $\binom{M}{2}$ ternary grid preparations of Theorem~\ref{thm:finite-grid-gibbs}, then~\eqref{eq:coarse-gibbs-bound} holds for every tested target spectrum whose pairwise gaps lie on $\Delta\mathbb{Z}$.
\end{corollary}
\begin{proof}
By translation invariance, the isolated binary log-odds of a split pair with internal spectrum $y$ equals $h(y_1-y_2)$, so Theorem~\ref{thm:finite-control} implies $d_N(x)\le\theta_N$. Applying Theorem~\ref{thm:finite-window-gibbs} in the continuum case, or Theorem~\ref{thm:finite-grid-gibbs} in the discrete grid case, yields~\eqref{eq:coarse-gibbs-bound}.
\end{proof}

When Proposition~\ref{prop:finite-suppression} is used with $e_0\le e_N$, the numerator of~\eqref{eq:coarse-theta} becomes
\begin{equation}
\frac{1}{N+k-2}+3e_N+\frac{D_Ne_N}{h}+L_Nh,
\label{eq:resource-bound}
\end{equation}
For each fixed $k$, taking $e_N\to0$ with the control and conditioning bounds fixed, and then $h\to0$, leaves only the suppression contribution. Letting $k\to\infty$ then removes this contribution, provided the target conditioning bounds remain positive. This recovers exact conditional factorization; extending the hypotheses to all finite spectra gives the characterization of Theorem~\ref{thm:exact-replacement}.

\subsection{Discretizing the control rectangle}
\label{subsec:rectangle-grid}

If the probability residuals $r_i$ in condition~(i) of Theorem~\ref{thm:finite-control} are evaluated on a discrete grid over the parameter rectangle $(s,t)\in[0,1]\times[-h_0,h_0]$ with mesh spacings $(\Delta_s,\Delta_t)$ (including the boundary lines) and satisfy the Lipschitz modulus
\begin{equation}
|r_i(s,t)-r_i(s',t')|
\le A_s|s-s'|+A_t|t-t'|,
\label{eq:rectangle-modulus}
\end{equation}
then an empirical maximum residual $\widehat{e}$ measured within tolerance $\sigma_r$ extends to the continuum rectangle with bound
\begin{equation}
e=\widehat{e}+\sigma_r+\tfrac{1}{2}A_s\Delta_s+\tfrac{1}{2}A_t\Delta_t.
\label{eq:rectangle-grid-error}
\end{equation}
When each residual is obtained as the difference of two directly measured block probabilities with individual error at most $\sigma$, we have $\sigma_r\le 2\sigma$. The present finite-difference and interpolation estimates use the a priori bounds~\eqref{eq:internal-curvature} and~\eqref{eq:rectangle-modulus}, respectively. Appendix~\ref{app:finite-observations} shows that finite observations alone do not uniquely determine the global response; it does not assert that these particular regularity assumptions are necessary for every quantitative interpolation method.

\subsection{Analytic bounds for two-scale mixtures}
\label{subsec:coarse-numerics}

To illustrate the progression from coarse-graining residuals to conditional log-odds bounds and finally to finite-grid total-variation bounds, we analyze the symmetric two-scale canonical mixture with scales $r_\pm=1\pm\delta$ ($0<\delta<1$):
\begin{equation}
  \Psi_N^{(\delta)}(z)
  =\frac{1}{2}\left[
    \frac{\Psi_N^{\mathrm{can}}(r_-z)}{r_-}
    +\frac{\Psi_N^{\mathrm{can}}(r_+z)}{r_+}\right],
  \qquad \Psi_N^{\mathrm{can}}(z)=-\log\sum_{i=1}^N e^{-z_i}.
  \label{eq:two-scale-potential}
\end{equation}
Here every coarse-grained spectrum $T_a(z,y)$ retains the exact potential value $\Psi_2^{(\delta)}(y)$ defined in~\eqref{eq:two-scale-potential}.

\begin{proposition}[Analytic probability-error envelope]
\label{prop:two-scale-envelope}
For the two-scale family~\eqref{eq:two-scale-potential}, the response deviation from the unit-scale Gibbs distribution satisfies
\begin{equation}
  \|P^{N,(\delta)}(z)-\mathsf{G}_{1,N}(z)\|_1
     \le\frac{\delta^2}{4}\operatorname{osc}(z)^2,
  \qquad \operatorname{osc}(z)=\max_i z_i-\min_i z_i.
  \label{eq:two-scale-response-envelope}
\end{equation}
On any binary control rectangle~\eqref{eq:control-rectangle}, let $R_f$ denote the maximum oscillation of the resolved spectrum $S_a$ across its four vertices, and set
\begin{equation}
  R_c=\max_{z\in\{x,x^*\}}\operatorname{osc}(T_a(z,y))+h_0,\qquad
  R_y=\operatorname{osc}(y)+h_0,\qquad
  b_\Psi=\delta^2\left[\frac{R_y^2}{8r_-}+\frac{\log2}{r_-^3}\right].
  \label{eq:two-scale-potential-envelope}
\end{equation}
Then every coarse-graining probability residual on $\mathcal{C}$ is bounded uniformly by
\begin{equation}
  e=\frac{\delta^2}{8}(R_f^2+R_c^2)+\frac{b_\Psi}{4},
  \label{eq:two-scale-rectangle-envelope}
\end{equation}
and the internal curvature constant in Theorem~\ref{thm:finite-control} may be taken as $L=3/4$.
\end{proposition}

\begin{proof}
For the canonical family $p_i(r)=e^{-rz_i}/\sum_j e^{-rz_j}$, direct differentiation gives $p_i''(r)=p_i(r)[(z_i-\langle z\rangle_r)^2-\operatorname{Var}_r(z)]$, so $\sum_i|p_i''(r)|\le 2\operatorname{Var}_r(z)\le\operatorname{osc}(z)^2/2$ by Popoviciu's inequality. Taylor expansion of the symmetric average $\frac{1}{2}[p(1-\delta)+p(1+\delta)]$ around $r=1$ yields~\eqref{eq:two-scale-response-envelope}, and half of this $\ell^1$ bound controls the probability discrepancy of any single event.

Similarly, for a binary spectrum $w$, differentiating $\phi_w(r)=r^{-1}\Psi_2^{\mathrm{can}}(rw)$ twice gives $\phi_w''(r)=-r^{-1}\operatorname{Var}_r(w)-2r^{-3}H(p(r))$, where $H(p)=-\sum_{i=1}^2 p_i\log p_i\le\log 2$. Since $\operatorname{Var}_r(w)\le\operatorname{osc}(w)^2/4$, symmetric Taylor expansion on $[r_-,r_+]$ bounds $|\Psi_2^{(\delta)}(w)-\Psi_2^{\mathrm{can}}(w)|\le b_\Psi$ whenever $\operatorname{osc}(w)\le R_y$.

By convexity of $\operatorname{osc}(\cdot)$, the resolved oscillation on $\mathcal{C}$ is maximized at the four vertices, giving $\operatorname{osc}(S_a)\le R_f$, while $|\partial_{y_1}\Psi_2^{(\delta)}|\le 1$ implies $\operatorname{osc}(T_a)\le R_c$ on $\mathcal{C}$. Decomposing each residual $r_i$ by inserting the unit-scale canonical resolved and coarse-grained probabilities (which cancel identically when evaluated at the canonical block potential $\Psi_2^{\mathrm{can}}$) bounds $|r_i|$ by the sum of the resolved mixture error $\delta^2 R_f^2/8$, the coarse mixture error $\delta^2 R_c^2/8$, and the shift in canonical coarse probability induced by $|\Psi_2^{(\delta)}-\Psi_2^{\mathrm{can}}|\le b_\Psi$. Because $|\partial_j \mathsf{G}_{1,n,i}|\le 1/4$, the last contribution is at most $b_\Psi/4$, establishing~\eqref{eq:two-scale-rectangle-envelope}. Finally, each diagonal Hessian entry of $\Psi_m^{(\delta)}$ is bounded in magnitude by $(r_-+r_+)/8=1/4$, so~\eqref{eq:curvature-from-potentials} yields $L=3/4$.
\end{proof}

We apply Proposition~\ref{prop:two-scale-envelope} first to the target spectrum $(0,1/2,1)$, represented as a binary split of its first two levels with $x=(0,1)$ and $y=(0,1/2)$, using the suppressing field $x^*=(12,1)$ ($D=12$) and internal half-width $h_0=1/2$. The step size is chosen as $h=\min\{h_0,\sqrt{De/L}\}$. To obtain a full finite-grid total-variation bound via Corollary~\ref{cor:coarse-to-gibbs}, we take $R=1$ and $M=4$, giving the six ternary grid spectra $\mathcal{T}_3=\{(0,j/4,k/4):1\le j<k\le 4\}$. Assigning each of the three pairs in every ternary spectrum its own control rectangle with $D=12$, $h_0=1/2$, and $L=3/4$, and using the analytical lower bounds
\[
  \mu_3=\frac{2}{2+e^{(1+\delta)R}},\qquad
  \rho_3=\frac{1}{1+e^{(1+\delta)R}},
\]
yields the uniform grid bound $\max_{z\in\mathcal{T}_3}\operatorname{TV}(P^{3,(\delta)}(z),\mathsf{G}_{a_R,3}(z))\le\tanh(7\theta_3/4)$ shown in Figure~\ref{fig:coarse-control}. The true conditional and total-variation discrepancies are $O(\delta^2)$ as $\delta\to0$. The square-root contribution $\sqrt{DLe}$ in Theorem~\ref{thm:finite-control} instead decreases linearly with $\delta$ when it dominates the bound. For smaller $\delta$, the fixed suppression error $\eta$ determines the limiting bound.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{figures/coarse_control_bounds.pdf}
\caption{Coarse-graining error propagation for the symmetric two-scale mixture $r_\pm=1\pm\delta$ with control displacement $D=12$ and internal half-width $h_0=1/2$. (a)~True conditional log-odds discrepancy for the binary split of $(0,1/2,1)$ and the theoretical upper bound of Theorem~\ref{thm:finite-control}. (b)~Analytic continuum residual envelope $e$ from Proposition~\ref{prop:two-scale-envelope} compared against the maximum residual sampled on a $41\times41$ grid over $\mathcal{C}$. (c)~Maximum total-variation distance across the six ternary grid spectra ($R=1, M=4$) and the resulting bound $\tanh(7\theta_3/4)$ from Corollary~\ref{cor:coarse-to-gibbs}.}
\label{fig:coarse-control}
\end{figure}


\section{Reservoir models and temperature fluctuations}
\label{sec:physical-models}

We now examine how concrete physical reservoir models and temperature-fluctuation mechanisms behave under the two integrability conditions~\eqref{eq:closure} and~\eqref{eq:descent}. Setting $k_B=1$ and $x=\beta E$, we quantify violations of work-form closure on $x$-dependent responses via the antisymmetric defect
\[
  \mathcal{R}_{ij}(x)
  =\partial_{x_j}P_{N,i}(x)-\partial_{x_i}P_{N,j}(x),
\]
and violations of temperature-energy scaling via the Euler scaling operator
\[
  D=\beta\partial_\beta-\sum_{i=1}^N E_i\partial_{E_i}.
\]
The conditions $\mathcal{R}=0$ and $DP_N=0$ probe distinct physical properties: the former tests whether state populations are generated by a reversible work potential, whereas the latter tests whether $E_i$ and $\beta^{-1}$ are the sole energy scales governing the reduced state.

\subsection{Microcanonical reservoir with constant heat capacity}
\label{subsec:constant-capacity}

Consider a weakly coupled finite thermal reservoir whose density of states $\omega_B(u)$ defines the entropy $S_B(u)=\log[\omega_B(u)/\omega_B(u_0)]=C\log(u/u_0)$ for $u>0$. In this convention, the temperature is $T_B=(S_B')^{-1}=u/C$ and the heat capacity is $C>0$. For the same density of states, the phase-space-volume entropy used by Campisi et al.~\cite{CampisiTalknerHanggi2009Bath} gives heat capacity $C+1$, explaining the corresponding unit shift in the power-law exponent. Calibrating the total microcanonical energy $a=C/\beta$ so that the baseline reservoir inverse temperature is $S_B'(a)=\beta$ independent of the system spectrum $E$ produces the weights
\begin{equation}
  P_{N,i}^{(C)}(x)
  =\frac{(1-x_i/C)^C}{\sum_j(1-x_j/C)^C},
  \qquad x_i<C.
  \label{eq:capacity-response}
\end{equation}

\Needspace{15\baselineskip}
\begin{proposition}[Closure defect for a fixed-energy-calibrated reservoir]
\label{prop:capacity-defect}
For $N\ge2$, the finite-reservoir response~\eqref{eq:capacity-response} satisfies $DP_N^{(C)}=0$, but its closure defect for $i\ne j$ is
\begin{equation}
  \mathcal{R}_{ij}^{(C)}(x)
  =\frac{P_{N,i}^{(C)}(x)P_{N,j}^{(C)}(x)(x_j-x_i)}
  {C(1-x_i/C)(1-x_j/C)}.
  \label{eq:capacity-defect}
\end{equation}
Consequently, $P_N^{(C)}$ is not generated by a response potential on any open domain containing nondegenerate spectra. As $C\to\infty$, writing $p_i(x)=e^{-x_i}/\sum_j e^{-x_j}$, one has uniformly on compact subsets of $\mathbb{R}^N$:
\begin{align}
  P_{N,i}^{(C)}(x)
  &=p_i(x)\left[
     1-\frac{x_i^2-\sum_jp_j(x)x_j^2}{2C}
    \right]+O(C^{-2}),
    \label{eq:capacity-asymptotic}\\
  \mathcal{R}_{ij}^{(C)}(x)
  &=\frac{p_i(x)p_j(x)(x_j-x_i)}{C}+O(C^{-2}).
    \label{eq:capacity-defect-asymptotic}
\end{align}
\end{proposition}

\begin{proof}
Because $P_N^{(C)}$ depends on $(\beta,E)$ solely through $x=\beta E$, $DP_N^{(C)}=0$ holds identically. Writing $\phi_C(t)=C\log(1-t/C)$ with $\phi_C'(t)=-(1-t/C)^{-1}$, differentiation gives $\partial_{x_j}P_{N,i}^{(C)}=P_{N,i}^{(C)}P_{N,j}^{(C)}(1-x_j/C)^{-1}$ for $i\ne j$, and antisymmetrizing yields~\eqref{eq:capacity-defect}. Expanding $\phi_C(t)=-t-t^2/(2C)+O(C^{-2})$ uniformly on compact sets yields~\eqref{eq:capacity-asymptotic}--\eqref{eq:capacity-defect-asymptotic}.
\end{proof}

Physically, the failure of closure in~\eqref{eq:capacity-defect} stems from holding the total energy $a=C/\beta$ fixed while shifting the system spectrum: a uniform shift $E_i\mapsto E_i+d$ reduces the residual reservoir energy $a-E_i$ unless $a$ is shifted simultaneously by $d$, violating translation invariance. Notice also that for any reservoir with positive heat capacity ($S_B''(a)=-\beta_B^2/C_B<0$), the log-weight $\log w_B(E)=S_B(a-E)$ is strictly concave in $E$, whereas for any nondegenerate positive superposition of Boltzmann factors $B(E)=\int_0^\infty e^{-bE}\,\nu(\mathrm{d}b)$,
\[
  \frac{\mathrm{d}^2}{\mathrm{d}E^2}\log B(E)
  =\operatorname{Var}_{\nu_E}(b)\ge0,
  \qquad
  \nu_E(\mathrm{d}b)=\frac{e^{-bE}\nu(\mathrm{d}b)}{B(E)}.
\]
The opposite signs of the second derivatives preclude representing finite-bath microcanonical weights as positive temperature mixtures (a distinction also noted in explicit finite Ising-chain calculations~\cite{FariasDavis2023}).

\subsection{Spectrum-dependent reservoir calibration}
\label{subsec:spectrum-dependent-reservoir}

If instead the total energy calibration is adjusted with the arithmetic mean of the system energy levels so as to preserve translation invariance, a finite microcanonical bath can generate an exact noncanonical EIE potential.

\begin{proposition}[Spectrum-calibrated non-Gibbs potential]
\label{prop:spectrum-dependent-reservoir}
For a reservoir with entropy $S_B(u)=\log(u/u_0)$ ($C=1$), calibrating the total energy as
\[
  a(\beta,E)=\overline{E}+\beta^{-1},
  \qquad
  \overline{E}=\frac{1}{N}\sum_{i=1}^N E_i,
\]
yields on the domain $\max_i(x_i-\overline{x})<1$ the affine response
\begin{equation}
  P_{N,i}(x)=\frac{1-(x_i-\overline{x})}{N},
  \qquad
  \overline{x}=\frac{1}{N}\sum_{i=1}^N x_i,
  \label{eq:affine-reservoir-response}
\end{equation}
which satisfies both EIE conditions~\eqref{eq:closure}--\eqref{eq:descent} and is generated by the following quadratic potential, strictly concave on $H_N$:
\begin{equation}
  \Psi_N(x)=\overline{x}
      -\frac{1}{2N}\sum_{i=1}^N(x_i-\overline{x})^2,
  \qquad
  \nabla^2\Psi_N=-\frac{1}{N}\left(I-\frac{1}{N}\mathbf{1}\mathbf{1}^\top\right).
  \label{eq:affine-reservoir-potential}
\end{equation}
\end{proposition}

\begin{proof}
Since $e^{S_B(a-E_i)}\propto a-E_i=\beta^{-1}[1-(x_i-\overline{x})]$ and $\sum_{i=1}^N[1-(x_i-\overline{x})]=N$, normalizing gives~\eqref{eq:affine-reservoir-response}, which is the gradient of~\eqref{eq:affine-reservoir-potential}.
\end{proof}

Because $a(\beta,E+d\mathbf{1})=a(\beta,E)+d$, this calibration keeps the reservoir excitations $a-E_i$ invariant under uniform shifts of the energy zero. The quadratic entropy $s_N(p)=-\frac{N}{2}\sum_{i=1}^N(p_i-1/N)^2$ on $\Delta_N$ generates a global concave extension of the local potential~\eqref{eq:affine-reservoir-potential} through Fenchel duality. Its gradient is the Euclidean projection of $N^{-1}(\mathbf{1}-x)$ onto $\Delta_N$. It agrees with~\eqref{eq:affine-reservoir-response} on the stated positivity domain; outside that domain, some populations can vanish, so the extension need not belong to the strictly positive global response class used in Section~\ref{sec:exact}.

\subsection{Hamiltonians of mean force at fixed coupling}
\label{subsec:fixed-coupling}

For a system coupled to a canonical bath at inverse temperature $\beta$ with total Hamiltonian $H_{\mathrm{tot}}=H_S+H_B+V$, the Hamiltonian of mean force $H^*$ is defined by~\cite{Jarzynski2004,TalknerHanggi2020}
\[
  e^{-\beta H^*}
  =\frac{\operatorname{Tr}_B e^{-\beta H_{\mathrm{tot}}}}{Z_B},
  \qquad Z_B=\operatorname{Tr}_B e^{-\beta H_B}.
\]
When the bare system levels $E_i$ are varied while holding the bath Hamiltonian $H_B$ and interaction $V$ fixed, the bare populations automatically satisfy the closure condition~\eqref{eq:closure}, even for noncommuting quantum interactions.

\begin{proposition}[Closure at fixed system--bath coupling]
\label{prop:mean-force-closedness}
Let $H_{\mathrm{tot}}(E)=\sum_{i=1}^N E_i\Pi_i+H_B+V$, where the orthogonal projectors $\Pi_i$ sum to the identity and, together with the bath Hamiltonian $H_B$ and coupling $V$, are independent of $E$. Assume the partition functions are finite and the trace derivatives below are justified, as is automatic in finite dimension. Then the mean-force free energy $F^*(\beta,E)=-\beta^{-1}\log(Z_{\mathrm{tot}}/Z_B)$ generates the bare populations $P_{N,i}=\operatorname{Tr}(\rho_{\mathrm{tot}}\Pi_i)$:
\begin{equation}
  \partial_{E_i}F^*=P_{N,i},
  \qquad
  \partial_{E_j}P_{N,i}=\partial_{E_i}P_{N,j}.
  \label{eq:mean-force-closedness}
\end{equation}
\end{proposition}

\begin{proof}
By Duhamel's formula for the derivative of the matrix exponential and cyclicity of the trace, $\partial_{E_i}Z_{\mathrm{tot}}=-\beta\operatorname{Tr}(\Pi_i e^{-\beta H_{\mathrm{tot}}})$. Since $Z_B$ does not depend on $E$, differentiating $F^*$ gives $\partial_{E_i}F^*=P_{N,i}$, and symmetry of mixed partial derivatives yields~\eqref{eq:mean-force-closedness}.
\end{proof}

Despite satisfying~\eqref{eq:mean-force-closedness}, fixed-coupling models generally violate the scaling equation~\eqref{eq:descent} because $H_B$ and $V$ introduce fixed energy scales that are not rescaled with $E$. As a concrete solvable example, consider a classical two-state bath $b\in\{0,1\}$ with joint energies $H(i,b)=E_i+(\Delta+g_i)b$, where $\Delta>0$ and $h_i=\Delta+g_i>0$. Writing $n_i=(1+e^{\beta h_i})^{-1}$ for the conditional bath excitation probability, we obtain
\begin{align}
  P_{N,i}(\beta,E)
  &=\frac{e^{-\beta E_i}(1+e^{-\beta h_i})}
          {\sum_j e^{-\beta E_j}(1+e^{-\beta h_j})},
    \label{eq:binary-bath-response}\\
  DP_{N,i}
  &=-\beta P_{N,i}
       \left(h_i n_i-\sum_{j=1}^N P_{N,j}h_jn_j\right),
    \label{eq:binary-bath-scaling}\\
  D\log\frac{P_{N,i}}{P_{N,j}}
  &=-\beta(h_i n_i-h_j n_j).
    \label{eq:binary-bath-log-odds}
\end{align}
At a fixed temperature, $DP_N=0$ precisely when all quantities $h_i n_i$ coincide. Distinct state-dependent couplings therefore produce a scaling defect generically, although accidental cancellations can occur at particular temperatures. Furthermore, the thermodynamic internal energy conjugate to $F^*$,
\begin{equation}
  U^*=\partial_\beta(\beta F^*)
     =\langle H^*+\beta\partial_\beta H^*\rangle
     =\sum_{i=1}^N P_{N,i}(E_i+h_in_i)-\frac{\Delta}{1+e^{\beta\Delta}},
  \label{eq:mean-force-energy}
\end{equation}
differs from the bare spectral expectation $U_{\mathrm{bare}}=\sum_i P_{N,i}E_i$ by the shift in bath and interaction energy. Figure~\ref{fig:physical-compatibility} contrasts the closure defect of the constant-capacity microcanonical reservoir against the scaling defect and energy renormalization of this fixed-coupling model.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{figures/physical_compatibility.pdf}
\caption{Complementary integrability defects for finite-capacity and fixed-coupling reservoirs on the spectrum $E=(0,0.7,1.8)$. (a)~Frobenius norm $\|\mathcal{R}^{(C)}\|_F$ of the closure defect and (b)~population distance $\|P^{(C)}-\mathsf{G}_{1,N}\|_1$ at $\beta=1$ for the constant-capacity reservoir~\eqref{eq:capacity-response}, compared with the leading $C^{-1}$ asymptotics~\eqref{eq:capacity-asymptotic}--\eqref{eq:capacity-defect-asymptotic} (dashed lines). (c)~Euclidean norm $\|DP_N\|_2$ of the scaling defect at $\beta=1$, $\Delta=1$, and $g=\lambda(0,0.4,1.2)$, and (d)~difference between the mean-force internal energy $U^*$ and bare energy $U_{\mathrm{bare}}$ as $\beta$ varies at $\lambda=1$ for the fixed-coupling model~\eqref{eq:binary-bath-response}, whose work form is closed.}
\label{fig:physical-compatibility}
\end{figure}

\subsection{Shared temperature fluctuations and mutual information}
\label{subsec:latent-temperature}

Suppose now that a system experiences slow fluctuations in its effective inverse temperature, described by a probability measure $\gamma$ on $(0,\infty)$ independent of $(\beta,E)$~\cite{BeckCohen2003,ZhaiSun2024}:
\[
  p_{N,i}(x)=\frac{e^{-x_i}}{\sum_j e^{-x_j}},
  \qquad
  P_{N,i}^{\gamma}(x)=\int_{(0,\infty)} p_{N,i}(rx)\,\gamma(\mathrm{d}r).
\]
By subtracting the zero-field value before integrating, the generating potential can be defined for an arbitrary probability measure $\gamma$ without requiring $\int r^{-1}\,\gamma(\mathrm{d}r)<\infty$.

\begin{proposition}[Centered mixture potential]
\label{prop:centered-mixture}
For every probability measure $\gamma$ on $(0,\infty)$, the integral
\begin{equation}
  \widetilde{\Psi}_N^\gamma(x)
  =\int_{(0,\infty)}\frac{\Psi_N^{\mathrm{can}}(rx)-\Psi_N^{\mathrm{can}}(0)}{r}\,\gamma(\mathrm{d}r),
  \qquad
  \Psi_N^{\mathrm{can}}(x)=-\log\sum_{i=1}^N e^{-x_i},
  \label{eq:centered-mixture-potential}
\end{equation}
defines a finite, concave, $C^1$, translation-covariant response potential with gradient $\nabla\widetilde{\Psi}_N^\gamma=P_N^\gamma$. If $\int_{(0,\infty)} r\,\gamma(\mathrm{d}r)<\infty$, then $\widetilde{\Psi}_N^\gamma$ is $C^2$ with negative semidefinite Hessian
\[
  \nabla^2\widetilde{\Psi}_N^\gamma(x)
  =\int_{(0,\infty)} r\bigl(p_N(rx)p_N(rx)^\top-\operatorname{diag}p_N(rx)\bigr)\,\gamma(\mathrm{d}r)\preceq 0.
\]
\end{proposition}

\begin{proof}
Since $\min_i x_i\le r^{-1}[\Psi_N^{\mathrm{can}}(rx)-\Psi_N^{\mathrm{can}}(0)]\le\max_i x_i$ and $p_{N,i}(rx)\in[0,1]$ uniformly in $r>0$, dominated convergence justifies differentiation under the integral sign. The Hessian integrand is bounded in operator norm by $2r$, so a finite first moment guarantees twice continuous differentiability.
\end{proof}

When two noninteracting subsystems $A$ and $B$ with spectra $x\in\mathbb{R}^N$ and $y\in\mathbb{R}^M$ are exposed to the same fluctuating temperature scale $R\sim\gamma$, their joint distribution
\[
  P_{ij}^{AB}=\int_{(0,\infty)} p_{N,i}(rx)p_{M,j}(ry)\,\gamma(\mathrm{d}r)
\]
becomes statistically correlated: for two identical binary systems with gap $z\ne 0$, $P_{11}^{AB}-P_1^A P_1^B=\operatorname{Var}_\gamma((1+e^{Rz})^{-1})$, which vanishes if and only if $\gamma$ is a Dirac mass. For weak fluctuations around unit scale, the resulting mutual information is governed universally by the product of the canonical heat capacities of the two subsystems.

\begin{proposition}[Mutual information from shared temperature fluctuations]
\label{prop:mutual-information-capacity}
Let $R=e^{\tau\xi}$, where $\xi$ is a bounded, symmetrically distributed random variable with $\mathbb{E}[\xi^2]=1$. For fixed spectra $x\in\mathbb{R}^N$ and $y\in\mathbb{R}^M$, let $p_i=p_{N,i}(x)$ and $q_j=p_{M,j}(y)$, and denote their reference canonical heat capacities by
\[
  C_A=\sum_{i=1}^N p_i\left(x_i-\sum_k p_kx_k\right)^2,
  \qquad
  C_B=\sum_{j=1}^M q_j\left(y_j-\sum_k q_ky_k\right)^2.
\]
Then, as $\tau\to 0$, the mutual information between $A$ and $B$ expands as
\begin{equation}
  I(A:B)=\frac{\tau^4}{2}C_A C_B+O(\tau^6),
  \label{eq:mutual-information-capacity}
\end{equation}
with remainder uniform for $(x,y)$ in compact sets.
\end{proposition}

\begin{proof}
Set $L_x=x\cdot\nabla_x$ and $L_y=y\cdot\nabla_y$. Because $\xi$ is symmetric, all odd moments vanish and Taylor expansion around $\tau=0$ yields the covariance matrix
\[
  \delta_{ij}:=P_{ij}^{AB}-P_i^A P_j^B
     =\tau^2(L_xp_i)(L_yq_j)+O(\tau^4).
\]
Since $\sum_i\delta_{ij}=\sum_j\delta_{ij}=0$, expanding the relative entropy $I(A:B)=\sum_{i,j}P_{ij}^{AB}\log(P_{ij}^{AB}/(P_i^A P_j^B))$ to second order in $\delta_{ij}$ and using $P_i^A=p_i+O(\tau^2)$, $P_j^B=q_j+O(\tau^2)$ gives
\[
  I(A:B)
   =\frac{1}{2}\sum_{i,j}\frac{\delta_{ij}^2}{P_i^AP_j^B}+O(\tau^6)
   =\frac{\tau^4}{2}
     \left(\sum_{i=1}^N\frac{(L_xp_i)^2}{p_i}\right)
     \left(\sum_{j=1}^M\frac{(L_yq_j)^2}{q_j}\right)+O(\tau^6).
\]
Substituting $L_xp_i=-p_i(x_i-\sum_k p_kx_k)$ identifies the two factors as $C_A$ and $C_B$.
\end{proof}

If the symmetry assumption on $\xi$ is relaxed to $\mathbb{E}[\xi]=0$ and $\mathbb{E}[\xi^2]=1$, with boundedness retained, the leading term $\frac{1}{2}\tau^4 C_A C_B$ remains unchanged with remainder $O(\tau^5)$. Figure~\ref{fig:common-temperature-information} illustrates the leading law and the predicted sixth-order remainder for the symmetric two-point distribution $\Pr(\xi=\pm 1)=1/2$.

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{figures/common_temperature_information.pdf}
\caption{Subsystem mutual information generated by a shared binary temperature fluctuation $R=e^{\pm\tau}$ with equal weights, for spectra $x=(0,0.8,1.9)$ and $y=(0,1.1)$. (a)~Exact mutual information $I(A:B)$ and the leading heat-capacity product law $\frac{1}{2}\tau^4 C_A C_B$ from Proposition~\ref{prop:mutual-information-capacity}. (b)~Absolute relative deviation $|I(A:B)/(\frac{1}{2}\tau^4 C_A C_B)-1|$, illustrating the $O(\tau^2)$ relative error and hence the $O(\tau^6)$ absolute remainder in~\eqref{eq:mutual-information-capacity}.}
\label{fig:common-temperature-information}
\end{figure}

\subsection{Boundary entropy and low-temperature asymptotics}
\label{subsec:low-temperature-entropy}

Spectral stability and concave duality do not by themselves imply finite entropy at pure ground states ($\beta\to\infty$). For example, the Cauchy response on a two-level system with gap $\Delta>0$ and $z=\beta\Delta$,
\[
  f(z)=\frac{1}{2}-\frac{1}{\pi}\arctan z,
\]
has Clausius entropy $S(z)=zf(z)-\int_0^z f(t)\,\mathrm{d}t=-\frac{1}{2\pi}\log(1+z^2)$ and heat capacity $C(z)=\frac{z^2}{\pi(1+z^2)}\to\frac{1}{\pi}$ as $\beta\to\infty$. Equivalently, its dual entropy on $\Delta_2$,
\[
  s_2(p)=\frac{1}{\pi}\log\sin(\pi p),
\]
diverges to $-\infty$ at the pure-state boundaries $p\in\{0,1\}$. By contrast, finite mixtures of positive canonical scales~\eqref{eq:finite-mixture-potential} possess bounded dual entropies throughout the closed simplex $\Delta_N$ and vanishing heat capacity as $\beta\to\infty$, while still violating the Shannon grouping law.


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

Our analysis establishes a direct quantitative link between operational coarse-graining and the canonical Gibbs distribution in finite systems. At the exact level, requiring a block of energy levels to be replaceable by its own response potential across all finite spectra forces conditional independence within the block and uniquely selects the Gibbs family with its canonical partition-function normalization. At the quantitative level, restricting observations to a bounded energy window of diameter $R$ shows that two- and three-state preparations suffice to calibrate a common inverse-temperature scale $a_R$, while pairwise conditional log-odds discrepancies control the total-variation distance to the Gibbs law uniformly in the target dimension $N$.

Theorem~\ref{thm:finite-control} bounds conditional log-odds discrepancies from coarse-grained block-probability residuals and a finite block-suppression experiment. The estimate separates the residual block mass $\eta$, the probability error $e$, the outer control displacement $D$, and the internal curvature bound $L$. Its denominator records the conditioning of the target probabilities through $\mu$ and $\rho$. Measurement and interpolation tolerances contribute to $e$. Increasing the auxiliary degeneracy $k$ reduces the suppression term to $(N+k-2)^{-1}$, at the cost of preparing more states and controlling the level shift $\Psi_k(0)$.

Work-form closure and temperature-energy scaling distinguish the reservoir models in Section~\ref{sec:physical-models}. A finite-capacity microcanonical bath with fixed total-energy calibration satisfies scaling but violates closure at order $C^{-1}$. A fixed-coupling mean-force model has a closed work form with respect to the bare levels, but its bath and interaction energies generally break scaling. The smooth temperature-scale mixtures considered here satisfy both integrability equations. A shared fluctuating temperature nevertheless correlates noninteracting subsystems: for weak bounded symmetric log-temperature fluctuations, the leading mutual information is $\tfrac12\tau^4 C_A C_B$. Integrability of each subsystem therefore remains compatible with a measurable failure of joint statistical independence.


\appendix

\section{Non-identifiability from finitely many observations}
\label{app:finite-observations}

The following construction shows that exact agreement with the canonical Gibbs family at finitely many spectra, including finitely many derivatives and coarse-graining replacement tests, does not uniquely determine a smooth stable response potential globally. This non-uniqueness is compatible with quantitative bounds between observations obtained from concavity or additional regularity assumptions, such as those used in Corollary~\ref{cor:finite-net-extension} and Section~\ref{subsec:rectangle-grid}.

\begin{proposition}[Indistinguishability under finite observations]
\label{prop:finite-observations}
Fix any finite collection of potential evaluations and finitely many derivative evaluations at finitely many spectra across finitely many dimensions (including all inner, coarse, and resolved spectra of finitely many replacement experiments), and suppose all recorded values agree with the unit-scale canonical Gibbs family $\Psi_{\mathrm{can},N}(x)=-\log\sum_{i=1}^N e^{-x_i}$. Then for any chosen dimension $N_0\ge 2$, there exists a smooth, symmetric, spectrally stable, strictly positive response family $\{\widetilde{\Psi}_N\}_{N\ge 1}$ that matches every recorded observation identically, has a finite dual entropy on each closed simplex $\Delta_N$, and yet differs from the Gibbs family in dimension $N_0$.
\end{proposition}

\begin{proof}
Project all observed $N_0$-dimensional spectra onto the zero-sum hyperplane $H_{N_0}=\mathbf{1}^\perp$ and take the union of their coordinate-permutation orbits; this forms a finite subset $\mathcal{O}\subset H_{N_0}$. Choose a smooth, nontrivial, permutation-invariant bump function $b:H_{N_0}\to\mathbb{R}$ with compact support disjoint from a neighborhood of $\mathcal{O}$, and define
\[
\widetilde{\Psi}_{N_0}(x)=\Psi_{\mathrm{can},N_0}(x)+\epsilon\, b(x-\overline{x}\mathbf{1}),
\qquad \overline{x}=\frac{1}{N_0}\sum_{i=1}^{N_0}x_i,
\]
while keeping $\widetilde{\Psi}_N=\Psi_{\mathrm{can},N}$ for all $N\ne N_0$. Translation covariance and permutation symmetry hold by construction. Because $\operatorname{supp}(b)$ is compact, $\nabla\Psi_{\mathrm{can},N_0}$ is bounded away from the boundary of $\Delta_{N_0}$ and $\operatorname{Hess}\Psi_{\mathrm{can},N_0}$ is uniformly negative definite on $H_{N_0}$ over $\operatorname{supp}(b)$. For sufficiently small $\epsilon\ne 0$, $\widetilde{\Psi}_{N_0}$ therefore retains a strictly positive gradient and a negative definite Hessian on $H_{N_0}$. Since $b$ vanishes in a neighborhood of $\mathcal{O}$, all recorded functional and derivative values---as well as all coarse-grained spectra constructed from them---coincide with their canonical values. The perturbation has a nonzero gradient somewhere, so the response differs from the unit-scale Gibbs law. Since it agrees with that law outside the compact projected support, it cannot be a Gibbs response at any other constant scale. Finally, since $|\widetilde{\Psi}_{N_0}-\Psi_{\mathrm{can},N_0}|$ is uniformly bounded,~\eqref{eq:duality} ensures that $s_{\widetilde{\Psi}_{N_0}}$ differs from the Shannon entropy by a bounded amount and hence remains finite on $\Delta_{N_0}$.
\end{proof}


\section{Pairwise probability ratios and sparse transition graphs}
\label{app:transition-graphs}

The pairwise log-odds condition used in Section~\ref{sec:exact} requires the ratio $P_i^N(x)/P_j^N(x)$ between every pair of states to depend solely on the energy gap $x_i-x_j$. In a continuous-time Markov jump process with transition rates $W_{i\leftarrow j}>0$ obeying detailed balance $P_j W_{i\leftarrow j}=P_i W_{j\leftarrow i}$, this all-pairs condition holds whenever every pair of states is connected by a transition whose forward-to-backward rate ratio is a universal function of the energy difference, $W_{i\leftarrow j}/W_{j\leftarrow i}=R(x_i-x_j)$.

By contrast, imposing gap-local rate ratios only along the edges of a sparse transition graph does not enforce the Gibbs distribution at a fixed spectrum. Consider, for example, the three-state linear chain $1\leftrightarrow 2\leftrightarrow 3$ with fixed spectrum $x=(0,1,3)$ and edge rate ratios $W_{i\leftarrow j}/W_{j\leftarrow i}=e^{-(x_i-x_j)^3}$ on the two allowed transitions $(1,2)$ and $(2,3)$. The unique stationary distribution satisfies detailed balance with weights
\[
P\propto (1,\,e^{-1},\,e^{-9}),
\]
since $(x_2-x_1)^3=1$ and $(x_3-x_2)^3=8$, so that no single inverse-temperature scale $a$ reproduces both population ratios $P_2/P_1=e^{-1}$ and $P_3/P_1=e^{-9}$. Complete connectivity alone at a fixed spectrum does not exclude non-exponential laws. If detailed balance with the same continuous gap-ratio function is required for every independently variable three-state spectrum, the cycle identity instead forces $R(u+v)=R(u)R(v)$ and hence $R(u)=e^{-au}$. On a chain, the same conclusion follows if the endpoint ratio is independent of the intermediate level for every choice of both endpoint energies. A separate sufficient condition is complete passivity, meaning passivity of every finite tensor power under arbitrary joint cyclic operations~\cite{PuszWoronowicz1978}. For a strictly positive state of a finite system with a nonconstant Hamiltonian, complete passivity selects a Gibbs state at nonnegative inverse temperature. If the Hamiltonian is proportional to the identity, every state is completely passive and no such selection follows.


\section{Numerical implementation and reproducibility}
\label{app:reproduction}

Online Resource~1 implements the finite-window calibration protocol (Section~\ref{sec:finite-window}), the coarse-graining control-rectangle bounds (Section~\ref{sec:coarse-error}), and the reservoir and fluctuating-temperature models (Section~\ref{sec:physical-models}). Spectra, mixture weights, grid resolutions, and sample sizes are specified in the text and in the programs. The multinomial simulations use seed 26052026; the random probability-distance checks use seed 27042026. Fixed seeds make sampling repeatable within a fixed software environment, whose versions are recorded in the supplied summaries. Floating-point results may differ slightly across platforms. The programs compare analytical derivatives with finite differences, test the scaling expression~\eqref{eq:binary-bath-scaling}, measure the decay of the mutual-information remainder, compute the finite-grid linear-programming optima, and check equality cases of Lemma~\ref{lem:odds-distance}. These numerical checks accompany the proofs and do not establish the continuum regularity assumptions used in the error bounds.


\backmatter
\section*{Statements and Declarations}
\subsection*{Funding}
Sun Zhenyu acknowledges support from Northeastern University under Grant No.~10011345052601.

\subsection*{Data and code availability}
Online Resource~1 (\texttt{ESM\_1.zip}) contains the Python source code, generated CSV data tables, figure-generation scripts, and numerical verification summaries supporting this study. Its \texttt{README.md} provides instructions for reproducing all numerical results and figures.

\subsection*{Competing interests}
The authors declare no competing interests.

\subsection*{Use of AI-assisted tools}
ChatGPT (OpenAI) was used for proofreading and language editing, mathematical consistency checks, assistance with local corrections to statements and proofs, and preparation and verification of the LaTeX and computational materials. The authors are responsible for the mathematical arguments, numerical results, and final text.

\bibliography{references}

\end{document}