The use of atom-centered basis functions (as opposed to real-space grids), and a commensurate atom-centered CAP that requires only standard, atom-centered DFT quadrature grids, makes the present methodology extensible to larger molecular systems, perhaps including models of condensed-phase systems. The cost of a single TDKS time step is comparable to the cost of a ground-state SCF cycle and the memory footprint is only twice that of a ground-state DFT calculation, due to the use of complex-valued matrices. The Gaussian orbital representation is thus essential for extending strongfield calculations to polyatomic molecules, and others have already reported applications including HHG spectra of CH 2 Cl 2 , 68 as well as strong-field ionization rates of C 6 H 8 , 57 both at the TD-CIS level. The present formalism suggests an efficient means to add correlation to those calculations, as an alternative to time-dependent implementations of correlated wave function models. In moving to larger systems, one may need to deal with a linear dependency problem arising when very diffuse basis sets are employed (such as 5-aug-cc-pVTZ that is used here), which we circumvent in the present work using very tight thresholds. Judicious use of floating centers ("ghost atoms") may be an alternative means to circumvent this problem. This extensibility may open the door to simulations of strong-field phenomena in a wide variety of complex systems.