The chemical Langevin equation (CLE) is a popular simulation method to probe the stochastic dynamics of chemical systems. The CLE’s main disadvantage is its break down in finite time due to the problem of evaluating square roots of negative quantities whenever the molecule numbers become sufficiently small. We show that this issue is not a numerical integration problem, rather in many systems it is intrinsic to all representations of the CLE. Various methods of correcting the CLE have been proposed which avoid its break down. We show that these methods introduce undesirable artefacts in the CLE’s predictions. In particular, for unimolecular systems, these correction methods lead to CLE predictions for the mean concentrations and variance of fluctuations which disagree with those of the chemical master equation. We show that, by extending the domain of the CLE to complex space, break down is eliminated, and the CLE’s accuracy for unimolecular systems is restored. Although the molecule numbers are generally complex, we show that the “complex CLE” predicts real-valued quantities for the mean concentrations, the moments of intrinsic noise, power spectra, and first passage times, hence admitting a physical interpretation. It is also shown to provide a more accurate approximation of the chemical master equation of simple biochemical circuits involving bimolecular reactions than the various corrected forms of the real-valued CLE, the linear-noise approximation and a commonly used two moment-closure approximation.