The study of gapped quantum many-body systems in three spatial dimensions has uncovered the existence of quantum states hosting quasiparticles that are confined, not by energetics but by the structure of local operators, to move along lower dimensional submanifolds. These so-called fracton phases are beyond the usual topological quantum field theory description, and thus require new theoretical frameworks to describe them. Here we consider coupling fracton models to topological quantum field theories in (3+1) dimensions by starting with two copies of a known fracton model and gauging the $mathbb{Z}_2$ symmetry that exchanges the two copies. This yields a class of exactly solvable lattice models that we study in detail for the case of the X-cube model and Haahs cubic code. The resulting phases host finite-energy non-Abelian immobile quasiparticles with robust degeneracies that depend on their relative positions. The phases also host non-Abelian string excitations with robust degeneracies that depend on the string geometry. Applying the construction to Haahs cubic code in particular provides an exactly solvable model with finite energy yet immobile non-Abelian quasiparticles that can only be created at the corners of operators with fractal support.