The erratic path to coalescence of LISA massive black hole binaries in sub-pc resolution simulations of smooth circumnuclear gas disks

We report on high-resolution simulations that explore the orbital decay of massive black hole (MBH) pairs with masses between $10^5$ and $10^7 M_{\odot}$ embedded in a circumnuclear gas disk (CND). An adiabatic equation of state is adopted, with a range of adiabatic indices, which maintains a smooth flow. Mergers between MBHs in this mass range would be detectable by the upcoming Laser Inteferometer Space Antenna (LISA). The orbital evolution is followed from the CND scale ($100$~pc) down to separations of $0.1$--$0.01$~pc at which a circumbinary disk (CBD) could form. The decay is erratic and strongly dependent on the gas flow within the disk, that ultimately determines the net torques experienced by the sinking MBH. Overall, we can identify three different evolutionary stages: (i) an initially slow decay that leads to no significant change in the orbital angular momentum, resulting in some circularization; (ii) a fast migration phase in which the orbital angular momentum decreases rapidly; and (iii) a final, very slow decay phase, in which orbital angular momentum can even increase, and a CBD can form. The fast migration phase owes to disk-driven torques originating primarily from the co-orbital region of the secondary MBH, at a distance of 1--3 Hill radii. We find strong analogies with fast Type III migration for massive planets in protoplanetary disks. The CBD forms only when the decay rate becomes small enough to allow it enough time to carve a cavity around the primary MBH, at scales $\lesssim 1$~pc; when this happens, the MBH separation nearly stalls in our higher-resolution run. We suggest an empirically modified gap-opening criterion that takes into account such timescale effects as well as other deviations from standard assumptions made in the literature. [Abriged]