  To investigate the novel quantum dynamic behaviors of magnetic materials that
arise from complex spin-spin interactions, it is necessary to probe the
magnetic response at a speed greater than the spin-relaxation and dephasing
processes. Recently developed two-dimensional (2D) terahertz magnetic resonance
(THz-MR) spectroscopy techniques use the magnetic components of laser pulses,
and this allows investigation of the details of the ultrafast dynamics of spin
systems. For such investigations, quantum treatment -- not only of the spin
system itself but also of the environment surrounding the spin system -- is
important. In our method, based on the theory of multidimensional optical
spectroscopy, we formulate nonlinear THz-MR spectra using an approach based on
the numerically rigorous hierarchical equations of motion. We conduct numerical
calculations of both linear (1D) and 2D THz-MR spectra for a linear chiral spin
chain. The pitch and direction of chirality (clockwise or anticlockwise) are
determined by the strength and sign of the Dzyaloshinskii-Moriya interaction
(DMI). We show that not only the strength but also the sign of the DMI can be
evaluated through the use of 2D THz-MR spectroscopic measurements, while 1D
measurements allow us to determine only the strength.
