Uniqueness and convergence for non-convex multi-matrix models with commutator-squared interaction

We study D ≥ 3 Hermitian matrices with potential V = ½ Σ X_i² + (g/4) Σ X_i⁴ − (h/4) Σ_{i<j} [X_i, X_j]². This is the bosonic Yang–Mills matrix integral with a mass term and a quartic confining term, and for h > 0 it is not convex. We prove that the planar Schwinger–Dyson equations have exactly one tracial solution, and that the empirical distribution of the matrix model converges to it almost surely. This holds in an explicit region of the (g, h) plane for every D. At g = 1 the region contains: D 3 4 5 6 0 ≤ h ≤ 0.491 0.294 0.208 0.160 Without any computer assistance, the same argument gives 0 ≤ h ≤ 0.243 for D = 3. The proof extends the method of the two-matrix paper (doi:10.5281/zenodo.24129786). Two estimates are new: Spectator lemma. Along the flow that generates the moments of (X_i + iX_j)*(X_i + iX_j), the commutators with the other matrices can only grow. Hence every solution obeys the two-matrix bound on ‖X_i + iX_j‖, for every pair. Hessian bound. Exact sum-of-squares certificates for D = 3, 4, 5, 6 bound the Hessian of the quartic part below by −c_D max‖X_i + iX_j‖². Here c_D = 0.545, 0.845, 1.15, 1.46, against 0.765, 1.265, 1.765, 2.265 from the two-matrix estimate. Averaging over subsets extends the bound to all D. We also explain, with numerical evidence, why the route through the Dirac-type operator Σ X_i ⊗ σ_i gives less at D = 3. Contents: the paper (PDF and a single LaTeX file); the four certificates (exact rationals, JSON); their checker, verify_dmatrix.py (Python 3 standard library only; 0.2–7 s per certificate); 17 tests per certificate and the exact range computation; the semidefinite programs and exact rounding that produced the certificates; the scripts behind the numerical remarks. Running ./check_all.sh verifies everything in about 1.5 minutes and ends with "ALL CHECKS PASSED". Licensed CC BY 4.0. Ilia Marghia; made in collaboration with Claude (Anthropic 2026).