\section{A coupled Cartesian scheme and its exact energy law}
\label{sec:coupled-energy}

\subsection{Known geometry at each linear solve}
Choose $0=t_0<\cdots<t_N=T$ and set $\tau_n=t_n-t_{n-1}$.
Given $q_0$ and a compatible exact initial state $X_0$, initialize
by the joint mass projection $X_h^0=P_M^0X_0$ at the pose $q_0$.
For $n\ge1$, first update the pose using the previous computed body
velocity:
\begin{equation}\label{eq:pose-update}
 c_n=c_{n-1}+\tau_n V_h^{n-1},\qquad
 \theta_n=\theta_{n-1}+\tau_n\omega_h^{n-1}.
\end{equation}
The angle is represented as a real-valued continuous angle.
The piecewise affine pose reconstruction has derivative
$\dot q_h=(V_h^{n-1},\omega_h^{n-1})$ on the $n$th interval.
The coupled velocity solve below therefore uses a known pose, while
its new body velocity determines the following pose update.

At $q_n$ choose any admitted Cartesian enrichment, kernel width,
volume-cell partition, and surface-panel partition. Write
\[
 \mathbb K_h^n=\mathbb K_h(q_n),\qquad
 P_M^n:\mathbb H\longrightarrow\mathbb K_h^n
\]
for the current joint space of \eqref{eq:joint-discrete-space} and
its orthogonal projector in \eqref{eq:mass-inner-product}.
All the indicated parameters may change at every time step. The
projection is in the combined fluid and rigid-body energy metric.
The same physical bulk mass integrals are used in every configuration.

Let $f^n$, $F_e^n$, and $T_e^n$ be approximations of the external
data on the interval. The primary choice is their exact time average
as ambient data; the fluid integral is taken on the \emph{current}
exterior region:
\begin{equation}\label{eq:discrete-load}
 \mathcal L^n(Y)=\int_{F(q_n)}f^n(x)\cdot v(x)\,dx
                         +F_e^n\cdot W+T_e^n\eta.
\end{equation}
The scheme finds $X_h^n=(u_h^n,V_h^n,\omega_h^n)\in\mathbb K_h^n$
such that
\begin{equation}\label{eq:coupled-scheme}
 \left(\frac{X_h^n-P_M^nX_h^{n-1}}{\tau_n},Y_h\right)_M
       +a(u_h^n,v_h)
       +\mathcal G(\omega_h^{n-1};X_h^n,Y_h)
       =\mathcal L^n(Y_h),\qquad Y_h\in\mathbb K_h^n.
\end{equation}
The previous angular velocity is a known coefficient of a skew
operator, whereas the current translational velocity is unknown in
that operator. Both fluid and rigid velocities are solved together.

For a current test, orthogonality gives
$(P_M^nX_h^{n-1},Y_h)_M=(X_h^{n-1},Y_h)_M$.
Thus \eqref{eq:coupled-scheme} can be assembled directly using the
cross-space mass pairing with the full previous state, without
performing a separate projection solve.
For each admissible current geometry the linear system has a unique
joint state: its bilinear form evaluated on $(Y_h,Y_h)$ equals
$\tau_n^{-1}\norm{Y_h}_M^2+a(v_h,v_h)>0$ for $Y_h\ne0$.
The skew term does not change this coercivity. This statement is
independent of the rank of the multiplier constraint matrix.

The geometric algorithm is considered until a configuration leaves
the collision-free region used by the analysis. Stability below holds
on every admissible prefix. The convergence argument closes the
geometry bound against the assumed collision-free exact branch for
sufficiently small discretization parameters; collision avoidance is
not an extra conclusion of the energy identity alone.

\subsection{Cancellation of the coupled reactions}
Opposite fluid and structural work is the standard multiplier
energy mechanism \citep{boffi2015,boffigastaldi2017}. Here it is
applied to both regularized volume and surface moments in the same
positive fluid--body mass space.
Let $b_h^n(Y_h,\Lambda_h)$ be the sum of the current volume and
surface moment pairings of
$\mathcal J_{\eps_n}v_h-\mathcal R_{q_n}(W_h,\eta_h)$ with the two
multiplier components. The equivalent mixed formulation on
$X_{h,\mathcal A^n}\times\R^3$ is
\begin{align}
 \left(\frac{X_h^n-X_h^{n-1}}{\tau_n},Y_h\right)_M
  +a(u_h^n,v_h)+\mathcal G(\omega_h^{n-1};X_h^n,Y_h)
  +b_h^n(Y_h,\Lambda_h^n)&=\mathcal L^n(Y_h),
       \label{eq:coupled-mixed}\\
 b_h^n(X_h^n,\Upsilon_h)&=0
       \qquad\text{for all }\Upsilon_h.\nonumber
\end{align}
The minus rigid velocity in each moment is essential: the same
reaction acts with opposite fluid and body work. Testing by $X_h^n$
annuls the entire reaction term exactly, including the interior
moments. A multiplier exists because the residual of the constrained
solution vanishes on the kernel of the finite-dimensional constraint
map and hence belongs to the range of its transpose. It may be
nonunique. No estimate of this multiplier is used in the energy or
velocity-convergence arguments.

\subsection{Exact identity and bounds independent of mesh history}
\begin{theorem}[Coupled discrete energy]
\label{thm:coupled-energy}
Every step of \eqref{eq:coupled-scheme} satisfies
\begin{equation}\label{eq:coupled-energy}
 \begin{split}
 \norm{X_h^n}_M^2-\norm{X_h^{n-1}}_M^2
  +\norm{X_h^n-X_h^{n-1}}_M^2
  +2\tau_n a(u_h^n,u_h^n)
       =2\tau_n\mathcal L^n(X_h^n).
 \end{split}
\end{equation}
Moreover,
\begin{equation}\label{eq:coupled-loss-splitting}
 \begin{split}
 \norm{X_h^n-X_h^{n-1}}_M^2
  =\norm{X_h^n-P_M^nX_h^{n-1}}_M^2
       +\norm{(I-P_M^n)X_h^{n-1}}_M^2.
 \end{split}
\end{equation}
Let $\ell_f^n\in U'$ denote the fluid part of \eqref{eq:discrete-load},
and put
\[
 A_N=\sum_{n=1}^N\tau_n
        \left(\frac{|F_e^n|^2}{m_e}
                         +\frac{|T_e^n|^2}{I_e}\right)^{1/2},
 \qquad
 B_N=\norm{X_h^0}_M^2+
                     \mu^{-1}\sum_{n=1}^N\tau_n\norm{\ell_f^n}_{U'}^2.
\]
Then
\begin{equation}\label{eq:coupled-stability}
 \begin{split}
 \max_{0\le n\le N}\norm{X_h^n}_M^2
   +\sum_{n=1}^N\norm{X_h^n-X_h^{n-1}}_M^2
   +\mu\sum_{n=1}^N\tau_n\norm{u_h^n}_U^2
       \le C(B_N+A_N^2),
 \end{split}
\end{equation}
with an absolute constant $C$. No time-step restriction, smallness
of the moment widths, bound on the number of mesh changes, or
multiplier inf--sup condition is required for this estimate.
\end{theorem}
\begin{proof}
Test \eqref{eq:coupled-scheme} with $X_h^n$. The projection may be
removed inside this pairing. The skew term is zero, and the identity
$2(X-Y,X)_M=\norm X_M^2-\norm Y_M^2+\norm{X-Y}_M^2$
proves \eqref{eq:coupled-energy}. The vector
$(I-P_M^n)X_h^{n-1}$ is orthogonal to the current space, which contains
$X_h^n-P_M^nX_h^{n-1}$. Pythagoras proves
\eqref{eq:coupled-loss-splitting}.

Write $\beta_n=(|F_e^n|^2/m_e+|T_e^n|^2/I_e)^{1/2}$.
Cauchy--Schwarz in the rigid mass norm gives
\[
 |F_e^n\cdot V_h^n+T_e^n\omega_h^n|
       \le\beta_n\norm{X_h^n}_M.
\]
Use $a(u_h^n,u_h^n)=\mu\norm{u_h^n}_U^2$ and apply Young's
inequality only to the fluid work. Summing the energy identity through
step $k$ gives
\begin{equation}\label{eq:prefix-energy-bound}
 \norm{X_h^k}_M^2
  +\sum_{n=1}^k\norm{X_h^n-X_h^{n-1}}_M^2
  +\mu\sum_{n=1}^k\tau_n\norm{u_h^n}_U^2
       \le B_N+2\sum_{n=1}^k\tau_n\beta_n\norm{X_h^n}_M.
\end{equation}
Let $Z=\max_{0\le n\le N}\norm{X_h^n}_M$. If the maximum is
attained at zero, $Z^2\le B_N$. Otherwise choose a maximizing index
in \eqref{eq:prefix-energy-bound} to obtain $Z^2\le B_N+2A_NZ$.
Hence
\[
 Z\le A_N+\sqrt{A_N^2+B_N}\le2A_N+\sqrt{B_N}.
\]
Use this bound in \eqref{eq:prefix-energy-bound} with $k=N$ to
control the two accumulated nonnegative terms. This proves
\eqref{eq:coupled-stability}. The argument is an algebraic maximum
estimate and does not introduce a time-step condition through a
discrete Gronwall denominator.
\end{proof}

For exact ambient time averages,
$\norm{\ell_f^n}_{U'}\le C\norm{f^n}_{L^2(\Omega)}$ and Jensen's
inequality bounds its squared time sum by
$C\norm f_{L^2(0,T;L^2(\Omega))}^2$.
Similarly $A_N$ is bounded by the $L^1$-in-time rigid dual-mass norm
of the external body load. In particular, the prescribed data and
positive excess inertias give a uniform bound for
$\max_n(|V_h^n|+|\omega_h^n|)$, and hence for the speed of the pose
reconstruction in \eqref{eq:pose-update}.

\begin{corollary}[Physical constraint defects]
\label{cor:coupled-defects}
Suppose the data bounds in Theorem~\ref{thm:coupled-energy} are
uniform over a sequence of admissible computations. Then
\begin{align}
 \left(\sum_n\tau_n
  \norm{u_h^n-\mathcal R_{q_n}\xi_h^n}_{L^2(B(q_n))}^2\right)^{1/2}
      &\le C(\bar\eps+\bar\kappa),\label{eq:integrated-volume-defect}\\
 \left(\sum_n\tau_n
  \norm{\gamma_{q_n}u_h^n-\mathcal R_{q_n}\xi_h^n}
                                  _{L^2(\Gamma(q_n))}^2\right)^{1/2}
      &\le C(\sqrt{\bar\eps}+\sqrt{\bar\ell}).
       \label{eq:integrated-surface-defect}
\end{align}
Here bars denote the maxima of the indicated widths over the run.
The surrogate-to-physical kinetic-energy discrepancy in
\eqref{eq:energy-discrepancy} tends to zero in the weighted $L^1$
time sum as $\bar\eps+\bar\kappa\to0$.
\end{corollary}
\begin{proof}
The stability bound controls $u_h$ in the weighted $L^2(U)$ norm.
It also controls $\xi_h$ uniformly by strict positivity of $m_e,I_e$.
Thus $\sum_n\tau_n\norm{X_h^n}_{\mathbb V}^2$ is bounded.
Apply Lemma~\ref{lem:volume-surface-defect} and
Proposition~\ref{prop:surrogate-physical-energy}.
\end{proof}

\begin{corollary}[Defects relative to the actual interpolated motion]
\label{cor:path-defects}
Under the uniform bounds of Corollary~\ref{cor:coupled-defects}, suppose
also that the piecewise affine pose path remains in the fixed
admissible configuration family. For $t\in(t_{n-1},t_n]$, define
\[
 U_h(t)=u_h^n,\qquad
 b_h^{\rm path}(t,x)=V_h^{n-1}
          +\omega_h^{n-1}J(x-c_h(t)).
\]
The latter is the material velocity of the reconstructed body, because
$\dot q_h=(V_h^{n-1},\omega_h^{n-1})$. Then
\begin{align}
 \left(\int_0^T
   \norm{U_h-b_h^{\rm path}}_{L^2(B(q_h(t)))}^2\,dt\right)^{1/2}
    &\le C(\bar\eps+\bar\kappa+\sqrt{\tau_{\max}}),
       \label{eq:path-volume-defect}\\
 \left(\int_0^T
   \norm{\gamma_{q_h(t)}U_h-b_h^{\rm path}}
                          _{L^2(\Gamma(q_h(t)))}^2\,dt\right)^{1/2}
    &\le C(\sqrt{\bar\eps}+\sqrt{\bar\ell}
                               +\sqrt{\tau_{\max}}).
       \label{eq:path-surface-defect}
\end{align}
The constants may depend on the uniform data bound and the fixed
admissible family, but not on the mesh history. The nodal defects in
\eqref{eq:integrated-volume-defect}--\eqref{eq:integrated-surface-defect}
use the newly solved velocity $\xi_h^n$; the present bounds instead
use the actual velocity of the affine pose path.
\end{corollary}
\begin{proof}
Let $w_n=u_h^n-L_{h,q_n}\xi_h^n$. The uniform mass bound controls
all body velocities and therefore the speed of $q_h$. If $R_{t,n}$
is the exact rigid map from $B(q_n)$ to $B(q_h(t))$, its displacement
and rotational change on bounded sets are at most $C\tau_n$.
Zero extension of $w_n$ and integration along this rigid map give
\[
 \|Ew_n\circ R_{t,n}-Ew_n\|_{L^2(\R^2)}
       \le C\tau_n\|w_n\|_U.
\]
The corresponding $H^1$ norm is bounded by $C\|w_n\|_U$.
A uniform multiplicative trace estimate therefore bounds the trace
difference on $\Gamma(q_n)$ by
$C\sqrt{\tau_n}\|w_n\|_U$, exactly as in the moving-trace argument
of Part~II \citep{partii2026}. The rigid map preserves volume and
arclength. Combining these shift estimates with the current volume
and surface defects bounds $w_n$ on the interpolated body by
$C(\eps_n+\kappa_n+\tau_n)\|X_h^n\|_{\mathbb V}$, and its
trace by
$C(\sqrt{\eps_n}+\sqrt{\ell_n}+\sqrt{\tau_n})
\|X_h^n\|_{\mathbb V}$.

The common lift reproduces its rigid affine field on every body in
the admissible family, including the intermediate body. At the same
physical point its difference from the path velocity obeys
\[
 |L_{h,q_n}\xi_h^n-b_h^{\rm path}(t)|
       \le C|\xi_h^n-\xi_h^{n-1}|+C\tau_n
 \quad\text{on }B(q_h(t))\cup\Gamma(q_h(t)).
\]
Indeed, the first term accounts for the change of rigid velocity and
the second for $|c_n-c_h(t)|\le C\tau_n$; all coefficients and body
diameters are uniformly bounded. The additional lag contribution is
controlled by the full joint jumps in the energy identity:
\[
 \sum_n\tau_n|\xi_h^n-\xi_h^{n-1}|^2
   \le C\tau_{\max}\sum_n
                  \|X_h^n-X_h^{n-1}\|_M^2\le C\tau_{\max}.
\]
Integrate over each interval, sum, and use
$\sum_n\tau_n\|X_h^n\|_{\mathbb V}^2\le C$.
The terms of order $\tau_{\max}$ are bounded by a constant depending
on $T$ times $\sqrt{\tau_{\max}}$. This proves both assertions.
\end{proof}

The discrete dissipation is $2\mu\int_\Omega|D(u_h)|^2$.
Interior strain can be nonzero at a finite resolution, as can the
physical surface slip. The energy law alone does not prove that this
interior strain vanishes; that conclusion uses the strong velocity
convergence established subsequently. No physical torque is inferred
from a raw discrete surface multiplier.