RUBIES: The Evolution of the Ionization Parameter from 0 < z < 9

High-redshift galaxies have smaller radii, harder ionizing continua, and higher ionizing photon production efficiencies than lower redshift systems, which implies a corresponding evolution in nebular conditions. A key metric to quantify gas properties is the ionization parameter, q, the ratio of the local ionizing photon flux to the local hydrogen density. The ionization parameter is often inferred from observed emission line ratios, e.g., O32=[O III]/[O II]. Prior to JWST, statistical samples of ionization parameter-sensitive emission lines in the rest-frame optical remained inaccessible at high-z. We investigate the dimensionless ionization parameter, U=q/c at 3<z<9, inferred using Cloudy photoionization models from the O32 ratios for 434 galaxies in the RUBIES survey with JWST/NIRSpec PRISM and G395M spectroscopy, constituting the largest high-z population study of U to date. We compare to lower-redshift samples from SDSS, LEGA-C, and KBSS to probe the evolution of U from 0<z<9. We find that U increases with redshift and specific star formation rate (sSFR), and decreases with stellar mass. We combine the predictive power with multivariate relations to estimate U from redshift, stellar mass, and sSFR for use in cases where O32 is not available from spectroscopy, and show that U increases with redshift even at fixed stellar mass and sSFR by a factor of ~4 from z=2 to z=6. Crucially, and in contrast to previous linear best-fit calibrations, our inference results in a systematic uncertainty in log U of ~0.3 dex at zero measurement uncertainty due to the wide range of photoionization models that predict the same O32 ratio without informative priors. Finally, we discuss future modeling frameworks to accept many observed emission lines to simultaneously constrain gas-phase abundances, densities, ionizing sources, and ionization parameters to high accuracies for individual galaxies.