*2.2. Direct Numerical Simulation*

In order to study the departure from the laminar Sexl–Womersley flow and the dynamics of intermittent localised turbulence, we perform direct numerical simulations (DNS) using our open-source pseudo-spectral simulation code nsPipe [22]. In nsPipe, the governing Equation (2) are treated in cylindrical coordinates (*r*, θ, *z*) and discretised using a Fourier–Galerkin ansatz in θ and *z* and high-order finite differences in *r*. No-slip boundary conditions are imposed at the solid pipe wall and periodic boundary conditions in θ and *z*. The discretised NSEs are integrated forward in time using a second-order predictor–corrector method with variable time-step size; details are given in López et al. [22] and the references therein. We have modified nsPipe to account for a time-dependent driving force that maintains a pulsating flow rate according to Equation (1) and an additional volume force *Fp* to perturb the flow locally.

We use a computational domain of 100*D* in length, and the number of radial grid points and Fourier modes used in our DNS is (*Nr* × *N*<sup>θ</sup> × *Nz*) = (96 × 192 × 2400). After dealiasing, this results in a spatial resolution of Δθ*R*<sup>+</sup> = 3.1 and Δ*z*<sup>+</sup> = 3.8, whereas radial grid points are clustered towards the pipe wall such that 0.06 <sup>≤</sup> <sup>Δ</sup>*r*<sup>+</sup> <sup>≤</sup> 1.4, and 14 points lie within the buffer layer based on the shear Reynolds number for the statistically steady case (*A* = 0). By comparison to the resolution used in other contemporary DNS studies in the literature and the fact that we consider only moderate *A* (the instantaneous Reynolds number is never > (1 + *A*)*Re*), we expect our choice of resolution to be sufficient for all set-ups considered here. The adaptive time-step size is roughly <sup>Δ</sup>*<sup>t</sup>* = <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>3</sup> *<sup>D</sup> us* .
