A nonlinear terahertz (THz) spectroscopy technique uses a sample illuminated by two THz pulses separately. The illumination generates two signals BA and BB, corresponding to the first and second THz pulse, respectively, after interaction with the sample. The interaction includes excitation of at least one ESR transition in the sample. The sample is also illuminated by the two THz pulses together, with an inter-pulse delay τ, generating a third signal BAB. A nonlinear signal BNL is then derived via BNL=BAB−BA−BB. This nonlinear signal BNL can be then processed (e.g., Fourier transform) to study the properties of the sample.