  Low-dimensional quantum magnets, particularly those with strong spin
frustration, are characterized by their notable spin fluctuations. Nuclear
magnetic resonance (NMR) serves as a sensitive probe of low-energy fluctuations
that offers valuable insight into rich magnetic phases and emergent phenomena
in quantum magnets. Although experimentally accessible, the numerical
simulation of NMR relaxation rates, specifically the spin-lattice relaxation
rate $1/T_1$, remains a significant challenge. Analytical continuation based on
Monte Carlo calculations are hampered by the notorious negative sign for
frustrated systems, and the real-time simulations incur significant costs to
capture low-energy fluctuations. Here we propose computing the relaxation rate
using thermal tensor networks (TTNs), which provides a streamlined approach by
calculating its imaginary-time proxy. We showcase the accuracy and versatility
of our methodology by applying it to one-dimensional spin chains and
two-dimensional lattices, where we find that the critical exponents $\eta$ and
$z\nu$ can be extracted from the low-temperature scalings of the simulated
$1/T_1$ near quantum critical points. Our results also provide insights into
the low-dimensional and frustrated magnetic materials, elucidating universal
scaling behaviors in the Ising chain compound CoNb$_2$O$_6$ and revealing the
renormalized classical behaviors in the triangular-lattice antiferromagnet
Ba$_8$CoNb$_6$O$_{24}$. We apply the approach to effective model of the family
of frustrated magnets AYbCh$_2$ (A = Na, K, Cs, and Ch = O, S, Se), and find
dramatic changes from spin ordered to the proposed quantum spin liquid phase.
Overall, with high reliability and accuracy, the TTN methodology offers a
systematic strategy for studying the intricate dynamics observed across a broad
spectrum of quantum magnets and related fields.
