Coupling-Robust Accuracy in Multiphysics Physics Informed Neural Networks via Kronecker-Preconditioned Optimization

Physics-informed neural networks (PINNs) for coupled multiphysics systems suffer systematic accuracy degradation as inter-equation coupling strengthens. We provide a theoretical explanation through neural tangent kernel (NTK) analysis: for linearly coupled systems, the standard NTK's spectral radius grows as Ω(γ^2) with coupling strength γ, shrinking the stable learning rate, while block-diagonal Gauss--Newton (GN) preconditioning yields a preconditioned NTK K_P = JH^{+}J^\top whose spectral radius is bounded by S (number of networks), independent of γ. Adam's diagonal preconditioning destroys this projector structure -- inflating λ_{\max} far above S for any coupling type -- and its residual-dynamics kernel grows as Θ(γ), placing its stable learning rate strictly between gradient descent and GN. For one-way coupling the limitation is class-wide: no diagonal preconditioner, fixed or adaptive, halves the driving residual in fewer than Ω(γ) iterations (Ω(γ^2) if fixed), whereas block-diagonal GN requires O(1). We verify Ω(γ^2) growth across linearly coupled benchmarks and confirm λ_{\max}(K_P) = S in all three 1D systems, including nonlinearly coupled NP+P. Combining the Kronecker-preconditioned optimizer SOAP with inverse-gradient-norm loss balancing (SOAP+GradNorm) yields coupling-robust accuracy: across 222 experiments spanning three 1D systems and a 2D electroosmotic flow benchmark, SOAP+GradNorm maintains final-epoch L_2 accuracy across coupling strengths, with \leq 2.3\times degradation in nonlinear NP+P while Adam+GradNorm fails (L_2 > 0.1). SOAP+GradNorm further scales to a 2D, 6-PDE electroosmotic flow at EDL-resolved conditions down to \varepsilon = 0.01 -- a regime all prior PINN electrokinetics studies have avoided -- where Adam+GradNorm fails entirely (L_2 > 0.3).