Massive black hole (MBH) binary inspiral time scales are uncertain, and their spins are even more poorly constrained. Spin misalignment, along with unequal mass ratios and spin magnitudes, introduces asymmetry in the gravitational radiation, which imparts a recoil kick to the merged MBH. Understanding how MBH binary spins evolve is crucial for determining their recoil velocities, their gravitational wave (GW) waveforms detectable with LISA, as well as their post-merger retention rate in galaxies and thus their subsequent merger rate. Here we present a novel study that introduces a sub-resolution model for gas- and GW-driven MBH binary spin evolution using a population of accreting MBHs from the Illustris cosmological hydrodynamics simulations. We also model sub-resolution binary inspiral via dynamical friction, stellar scattering, viscous gas drag, and GW emission. Our model assumes differential accretion, which causes greater alignment of the secondary MBH spin in unequal-mass mergers. We find that 47% of the MBHs in our population merge by $z=0$. Of these, 19% have misaligned primaries and 10% have misaligned secondaries at the time of merger in our (conservative) fiducial model. The MBH misalignment fraction depends strongly on the accretion disc parameters, however. Reducing accretion rates by a factor of 100, in a thicker disc, yields 79% and 42% misalignment for primaries and secondaries, respectively. Even in the fiducial model, more than 12% of binaries experience recoils of $>500$ km s$^{-1}$, which could displace them at least temporarily from galactic nuclei. We additionally find that a significant number of systems experience strong precession.

We numerically determine the cosmological branch of the free function in a nonlocal metric-based modification of gravity which provides a relativistic generalization of Milgrom's Modified Newtonian Dynamics. We are able to reproduce the $\Lambda$CDM expansion history over virtually all of cosmic history, including the era of radiation domination during Big Bang Nucleosynthesis, the era of matter domination during Recombination, and most of the era of vacuum energy domination. The very late period of $0 \leq z < 0.0880$, during which the model deviates from the $\Lambda$CDM expansion history, is interesting because it causes the current value of the Hubble parameter to be about 4.5\% larger than it would be for the $\Lambda$CDM model. This may resolve the tension between inferences of $H_0$ which are based on data from large redshift and inferences based on Hubble plots.