Lipschitz regularity of harmonic quasiconformal maps between Lyapunov domains in $\mathbb{R}^n$

We prove that every sense-preserving harmonic $K$--quasiconformal homeomorphism $f\colon D\to\Omega$ between Lyapunov domains (equivalently, bounded $C^{1,\alpha}$ domains) in $\mathbb{R}^n$, $\alpha\in(0,1]$, is globally Lipschitz on $\overline D$. The argument is based on a boundary iteration scheme: an initial H\"older modulus for the boundary trace (coming from quasiconformality) is improved via the $C^{1,\alpha}$ graph representation of $\partial\Omega$, yielding higher H\"older regularity for the normal component. This boundary gain is converted into a near-boundary gradient bound for harmonic functions through a basepoint boundary H\"older-to-gradient estimate obtained by flattening the boundary and using local harmonic-measure bounds. Quasiconformality then propagates the resulting control from one component to the full differential, and iteration gives boundedness of $|Df|$ up to the boundary. Along the way we briefly survey several standard tools from the theory of quasiconformal harmonic mappings (QCH), including boundary H\"older continuity, distortion of derivatives, and boundary-to-interior propagation principles that enter the iteration.