  Recently, a family of models that couple multifluid systems to the full
Maxwell equations draw a lot of attention in laboratory, space, and
astrophysical plasma modeling. These models are more complete descriptions of
the plasma than reduced models like magnetohydrodynamic (MHD) since they
naturally retain non-ideal effects like electron inertia, Hall term, pressure
anisotropy/nongyrotropy, etc. One obstacle to broader application of these
model is that an explicit treatment of their source terms leads to the need to
resolve rapid kinetic processes like plasma oscillation and electron cyclotron
motion, even when they are not important. In this paper, we suggest two ways to
address this issue. First, we derive the analytic forms solutions to the source
update equations, which can be implemented as a practical, but less generic
solver. We then develop a time-centered, locally implicit algorithm to update
the source terms, allowing stepping over the fast kinetic time-scales. For a
plasma with $S$ species, the locally implict algorithm involves inverting a
local $3S+3$ matrix only, thus is very efficient. The performance can be
further elevated by using the direct update formulas to skip null calculations.
Benchmarks illustrated the exact energy-conservation of the locally implicit
solver, as well as its efficiency and robustness for both small-scale,
idealized problems and large-scale, complex systems. The locally implicit
algorithm can be also easily extended to include other local sources, like
collisions and ionization, which are difficult to solve analytically.
