  A long-standing and formidable challenge faced by all conservative schemes
for relativistic magnetohydrodynamics (RMHD) is the recovery of primitive
variables from conservative ones. This process involves solving highly
nonlinear equations subject to physical constraints. An ideal solver should be
"robust, accurate, and fast -- it is at the heart of all conservative RMHD
schemes," as emphasized in [S.C. Noble et al., ApJ, 641:626-637, 2006]. Despite
over three decades of research, seeking efficient solvers that can provably
guarantee stability and convergence remains an open problem.
  This paper presents the first theoretical analysis for designing a robust,
physical-constraint-preserving (PCP), and provably (quadratically) convergent
Newton-Raphson (NR) method for primitive variable recovery in RMHD. Our key
innovation is a unified approach for the initial guess, devised based on
sophisticated analysis. It ensures that the NR iteration consistently converges
and adheres to physical constraints. Given the extreme nonlinearity and
complexity of the iterative function, the theoretical analysis is highly
nontrivial and technical. We discover a pivotal inequality for delineating the
convexity and concavity of the iterative function and establish theories to
guarantee the PCP property and convergence. We also develop theories to
determine a computable initial guess within a theoretical "safe" interval.
Intriguingly, we find that the unique positive root of a cubic polynomial
always falls within this interval. Our PCP NR method is versatile and can be
seamlessly integrated into any RMHD scheme that requires the recovery of
primitive variables, potentially leading to a broad impact in this field. As an
application, we incorporate it into a discontinuous Galerkin method, resulting
in fully PCP schemes. Several numerical experiments demonstrate the efficiency
and robustness of the PCP NR method.
