  The propagation of ultrafast pulses in dispersion-engineered waveguides,
exhibiting strong field confinement in both space and time, is a promising
avenue towards single-photon nonlinearities in an all-optical platform.
However, quantum engineering in such systems requires new numerical tools and
physical insights to harness their complicated multimode and nonlinear quantum
dynamics. In this work, we use a self-consistent, multimode Gaussian-state
model to capture the nonlinear dynamics of broadband quantum fluctuations and
correlations, including entanglement. Notably, despite its parametrization by
Gaussian states, our model exhibits nonlinear dynamics in both the mean field
and the quantum correlations, giving it a marked advantage over conventional
linearized treatments of quantum noise, especially for systems exhibiting gain
saturation and strong nonlinearities. Numerically, our approach takes the form
of a Gaussian split-step Fourier (GSSF) method, naturally generalizing highly
efficient SSF methods used in classical ultrafast nonlinear optics; the
equations for GSSF evaluate in $O(M^2\log M)$ time for an $M$-mode system with
$O(M^2)$ quantum correlations. To demonstrate the broad applicability of GSSF,
we numerically study quantum noise dynamics and multimode entanglement in
several ultrafast systems, from canonical soliton propagation in third-order
($\chi^{(3)}$) waveguides to saturated $\chi^{(2)}$ broadband parametric
generation and supercontinuum generation, e.g., as recently demonstrated in
thin-film lithium niobate nanophotonics.
