  We consider the numerical solution of the inertial version of
Landau-Lifshitz-Gilbert equation (iLLG), which describes high-frequency
nutation on top of magnetization precession due to angular momentum relaxation.
The iLLG equation defines a higher-order nonlinear dynamical system with very
different nature compared to the classical LLG equation, requiring twice as
many degrees of freedom for space-time discretization. It exhibits essential
conservation properties, namely magnetization amplitude preservation,
magnetization projection conservation, and a balance equation for generalized
free energy, leading to a Lyapunov structure (i.e. the free energy is a
decreasing function of time) when the external magnetic field is constant in
time. We propose two second-order numerical schemes for integrating the iLLG
dynamics over time, both based on implicit midpoint rule. The first scheme
unconditionally preserves all the conservation properties, making it the
preferred choice for simulating inertial magnetization dynamics. However, it
implies doubling the number of unknowns, necessitating significant changes in
numerical micromagnetic codes and increasing computational costs especially for
spatially inhomogeneous dynamics simulations. To address this issue, we present
a second time-stepping method that retains the same computational cost as the
implicit midpoint rule for classical LLG dynamics while unconditionally
preserving magnetization amplitude and projection. Special quasi-Newton
techniques are developed for solving the nonlinear system of equations required
at each time step due to the implicit nature of both time-steppings. The
numerical schemes are validated on analytical solution for macrospin terahertz
frequency response and the effectiveness of the second scheme is demonstrated
with full micromagnetic simulation of inertial spin waves propagation in a
magnetic thin-film.
