Recently, there has been remarkable progress in the complex Langevin method, which aims to solve the complex action problem by complexifying the dynamical variables in the original path integral. In particular, a new technique, called gauge cooling, has been introduced and the full QCD simulation at finite density has been made possible in the high-temperature (deconfined) phase or with heavy quarks. Here we provide an explicit justification of the complex Langevin method including the gauge cooling procedure. We first show that the gauge cooling can be formulated in the form of a modified complex Langevin equation involving a complexified gauge transformation, which is chosen appropriately as a function of the configuration before cooling. The probability distribution of the complexified dynamical variables is modified accordingly. However, this modification is shown not to affect the Fokker–Planck equation for the corresponding complex weight as long as observables are restricted to gauge-invariant ones. Thus we demonstrate explicitly that gauge cooling can be used as a viable technique to satisfy the convergence conditions for the complex Langevin method. We also discuss “gauge cooling” in 0D systems such as vector models or matrix models.