  We present a novel numerical scheme for the efficient and accurate solution
of the isothermal two-fluid (electron and ion) equations coupled to Poisson's
equation for low-temperature plasmas. The model considers electrons and ions as
separate fluids, comprising the electron inertia and charge separation. The
discretization of this system with standard explicit schemes is constrained by
very restrictive time steps and cell sizes related to the resolution of the
Debye length, electron plasma frequency, and electron sound waves. Both sheath
and electron inertia are fundamental to fully explain the physics in
low-pressure and low-temperature plasmas. However, most of the phenomena of
interest for fluid models occur at speeds much slower than the electron thermal
speed and are quasi-neutral, except in small charged regions. In this work, we
present a scheme based on the Lagrange-projection operator splitting that
preserves the asymptotic regime where the plasma is quasi-neutral with massless
electrons. As a result, the quasi-neutral regime is treated without the need of
an implicit solver nor the resolution of the Debye length and electron plasma
frequency. Additionally, the scheme proves to accurately represent the dynamics
of the electrons both at low speeds and when the electron speed is comparable
to the thermal speed. In addition, a well-balanced treatment of the ion source
terms is proposed in order to tackle problems where the ion temperature is very
low compared to the electron temperature. The scheme significantly improves the
accuracy both in the quasi-neutral limit and in the presence of plasma sheaths
when the Debye length is resolved. The scheme is assessed in two specifically
designed test-cases: a quasi-neutral two-stream periodic perturbation with
analytical solution and a low-temperature discharge that includes sheaths.
