Pulsatile fluid flows through straight pipes undergo a sudden transition to turbulence that is extremely difficult to predict. The difficulty stems here from the linear Floquet stability of the laminar flow up to large Reynolds numbers, well above experimental observations of turbulent flow. This makes the instability problem fully nonlinear and thus dependent on the shape and amplitude of the flow perturbation, in addition to the Reynolds and Womersley numbers and the pulsation amplitude.In our paper, we present an adjoint optimization code, based on a GPU implementation of the pseudo-spectral Navier--Stokes solver nspipe, which incorporates an automatic, optimal check-pointing strategy. We leverage this code to show that the flow is susceptible to two distinct instability routes: One in the deceleration phase, where the flow is prone to oblique instabilities, and another during the acceleration phase with similar mechanisms as in steady pipe flow. Instability is energetically more likely in the deceleration phase. This dataset includes optimal perturbations, their corresponding energy growth over time and post-processing scripts.