Second-order perturbations of the Schwarzschild spacetime: Practical, covariant, and gauge-invariant formalisms

High-accuracy gravitational-wave modeling demands going beyond linear, first-order perturbation theory. Particularly motivated by the need for second-order perturbative models of extreme-mass-ratio inspirals and black hole ringdowns, we present practical spherical-harmonic decompositions of the Einstein equation, Regge-Wheeler-Zerilli equations, and Teukolsky equation at second perturbative order in a Schwarzschild background. Our formulations are covariant on the $t\text{\ensuremath{-}}r$ plane and on the two-sphere, and we express the field equations in terms of gauge-invariant metric perturbations. In a companion Mathematica package, perturbationequations, we provide these invariant formulas as well as the analogous formulas in terms of raw, gauge-dependent metric perturbations. Our decomposition of the second-order Einstein equation, when specialized to the Lorenz gauge, was a key ingredient in recent second-order self-force calculations [Phys. Rev. Lett. 124, 021101 (2020); ibid. 127, 151102 (2021); ibid. 130, 241402 (2023)].