Second-order Sobolev regularity of three-dimensional \texorpdfstring{$p$}{p}-harmonic functions

For every finite $p\geq5$, a weak $p$-harmonic function on an arbitrary open subset of $\R^3$ belongs to $W^{2,2}_{\loc}$. The conclusion concerns the full unweighted distributional Hessian and includes neighborhoods of critical points without assumptions on the critical set. We prove an estimate for smooth regularized solutions in which the Hessian energy is divided by the actual gradient norm. Gaussian curvature of the level surfaces supplies the component that is lost in the usual pointwise coercivity estimate. The equation excludes closed level components in a ball, so Gauss--Bonnet and the total curvature of the boundary curves control the integrated correction. A radial cutoff argument makes the estimate uniform in the regularization parameter, and smooth approximation identifies every second distributional derivative.