The numerical simulation of moving boundaries or fluid-solid interfaces in incompressible viscous flow continues to remain a challenging task. Higher order methods offer the potential to achieve accurate results on computationally feasible grids with a minimum of numerical costs. However, constructing higher order numerical methods maintaining stability and inheriting most of the rich structure of the continuous problem becomes increasingly difficult. For discretizing the incompressible Navier--Stokes equations on time-dependent domains, we use arbitrary order discontinuous Galerkin methods in time along with inf-sup stable finite element pairs with discontinuous pressure in space. For the accurate representation of dynamically changing geometries, we employ a stabilized Nitsche fictitious domain method and cut finite element techniques [M. Anselmann, M. Bause, CutFEM and ghost stabilization techniques for higher order space-time discretizations of the Navier--Stokes equations, Int. J. Numer. Meth. Fluids, 94 (2022), pp. 775-802]. The geometry is immersed into an underlying computational background grid, which is not fitted to the geometric problem structure and kept fixed over the whole simulation time. Boundary and interface conditions are imposed in weak form. The stability, accuracy and performance properties of the approach are investigated carefully. In Newton iterations for the discrete problems, large linear systems with complex block structure are built, in particular if higher order discretizations and high spatial resolution are involved. Their iterative solution requires an efficient preconditioner that is tailored to the heterogeneous structure of discretization, involving non-cut, cut and ghost penalty elements. We propose using GMRES iterations that are preconditioned by a single step of the V-cycle geometric multigrid algorithm [M. Anselmann, M. Bause, A geometric multigrid method for space-time finite element discretizations of the Navier-Stokes equations and its application to 3d flow simulation, ACM Trans. Math. Softw., 49 (2023), Article No.: 5, pp. 1-25]. For the smoother, we use the local Vanka operator that is built on the tensor product of either a single finite element for non-cut elements or a patch of adjacent elements in the region of cut elements and a time subinterval. The performance properties of the solver are studied carefully and illustrated for a sequence of test problems.