  The complex Langevin (CL) method is a classical numerical strategy to
alleviate the numerical sign problem in the computation of lattice field
theories. Mathematically, it is a simple numerical tool to compute a wide class
of high-dimensional and oscillatory integrals. However, it is often observed
that the CL method converges but the limiting result is incorrect. The
literature has several unclear or even conflicting statements, making the
method look mysterious. By an in-depth analysis of a model problem, we reveal
the mechanism of how the CL result turns biased as the parameter changes, and
it is demonstrated that such a transition is difficult to capture. Our analysis
also shows that the method works for any observables only if the probability
density function generated by the CL process is localized. To generalize such
observations to lattice field theories, we formulate the CL method on general
groups using rigorous mathematical languages for the first time, and we
demonstrate that such localized probability density function does not exist in
the simulation of lattice field theories for general compact groups, which
explains the unstable behavior of the CL method. Fortunately, we also find that
the gauge cooling technique creates additional velocity that helps confine the
samples, so that we can still see localized probability density functions in
certain cases, as significantly broadens the application of the CL method. The
limitations of gauge cooling are also discussed. In particular, we prove that
gauge cooling has no effect for Abelian groups, and we provide an example
showing that biased results still exist when gauge cooling is insufficient to
confine the probability density function.
