Identifying quantum spin liquids, magnon breakdown, or fractionalized excitations in quantum magnets is an ongoing challenge due to the ambiguity in interpreting excitation continua occurring in linear response probes. Recently, it was proposed that techniques measuring higher-order response, such as two-dimensional coherent spectroscopy (2DCS), could resolve such ambiguities. Numerically simulating nonlinear response functions can, however, be computationally very demanding. We present an efficient Lanczos-based method to compute second-order susceptibilities $\chi^{2}(\omega_t,\omega_\tau)$ directly in the frequency domain. Applying this to extended Kitaev models describing α-RuCl3, we find qualitatively different nonlinear responses between intermediate magnetic field strengths and the high-field regime. To put these results into context, we derive the general 2DCS response of partially-polarized magnets within the linear spin-wave approximation, establishing that \chi^2(\omega_t,\omega_\tau)$ is restricted to a distinct universal form if the excitations are conventional magnons. Deviations from this form, as predicted in our Lanczos-based simulations for α-RuCl3, can hence serve in 2DCS experiments as direct criteria to determine whether an observed excitation continuum is of conventional two-magnon type or of different nature.

Numerical results that are plotted in the manuscript's figures. Also attached are python plot scripts that plot the respective data.