We report constraints on the three-dimensional orbital architecture for all four planets known to orbit the nearby M dwarf Gliese 876 based solely on Doppler measurements and demanding long-term orbital stability. Our data set incorporates publicly available radial velocities taken with the ELODIE and CORALIE spectrographs, High Accuracy Radial velocity Planet Searcher (HARPS), and Keck HIgh Resolution Echelle Spectrometer (HIRES) as well as previously unpublished HIRES velocities. We first quantitatively assess the validity of the planets thought to orbit GJ 876 by computing the Bayes factors for a variety of different coplanar models using an importance sampling algorithm. We find that a four-planet model is preferred over a three-planet model. Next, we apply a Newtonian Markov chain Monte Carlo algorithm to perform a Bayesian analysis of the planet masses and orbits using an N-body model in three-dimensional space. Based on the radial velocities alone, we find that a 99 per cent credible interval provides upper limits on the mutual inclinations for the three resonant planets (Fcbbecbbe < 7.87°, suggesting that the amount of planet-planet scattering in the system has been low. We investigate the distribution of the respective resonant arguments of each planet pair and find that at least one argument for each planet pair and the Laplace argument librate. The libration amplitudes in our three-dimensional orbital model support the idea of the outer three planets having undergone significant past disc migration.
