The exploitation of ultrafast electron dynamics in quantum cascade lasers (QCLs) holds enormous potential for intense, compact mode-locked terahertz (THz) sources, squeezed THz light, frequency mixers, and comb-based metrology systems. Yet the important sub-cycle dynamics have been notoriously difficult to access in operational THz QCLs. Here, we employ high-field THz pulses to perform the first ultrafast two-dimensional spectroscopy of a free-running THz QCL. Strong incoherent and coherent nonlinearities up to eight-wave mixing are detected below and above the laser threshold. These data not only reveal extremely short gain recovery times of 2 ps at the laser threshold, they also reflect the nonlinear polarization dynamics of the QCL laser transition for the first time, where we quantify the corresponding dephasing times between 0.9 and 1.5 ps with increasing bias currents. A density-matrix approach reproducing the emergence of all nonlinearities and their ultrafast evolution, simultaneously, allows us to map the coherently induced trajectory of the Bloch vector. The observed high-order multi-wave mixing nonlinearities benefit from resonant enhancement in the absence of absorption losses and bear potential for a number of future applications, ranging from efficient intracavity frequency conversion, mode proliferation to passive mode locking.