  First-principles calculation of nonlinear magneto-optical effects has become
an indispensable tool to reveal the geometric and topological nature of
electronic states and to understand light-matter interactions. While
intriguingly rich physics could emerge in magnetic materials, further
methodological developments are required to deal with time-reversal symmetry
breaking, due to the degeneracy and gauge problems caused by symmetry and the
low-frequency divergence problem in the existing calculation formalism. Here we
present a gauge-covariant and low-frequency convergent formalism for the
first-principles computation. Remarkably, this formalism generally works for
both non-magnetic and magnetic materials with or without band degeneracy.
Reliability and capability of our method are demonstrated by studying example
materials (i.e., bilayers of MnBi$_2$Te$_4$ and CrI$_3$) and comparing with
published results. Moreover, an importance correction term that ensures gauge
covariance of degenerate states is derived, whose influence on physical
responses is systematically checked. Our method enables computation of
nonlinear magneto-optical effects in magnetic materials and paves the way for
exploring rich physics created by the interplay of light and magnetism.
