We prove global Lipschitz stability for an inverse source problem of a system of wave equations on a Lorentzian manifold. The method used in this paper is widely known as the Bukhgeim–Klibanov method. However, the conventional method is not sufficient for the application to the hyperbolic partial differential equation with time-dependent coefficients to obtain the Lipschitz stability, and further innovations are needed. In this paper, we present an improved global Carleman estimate and an energy estimate to obtain the Lipschitz stability.